"Fossies" - the Fresh Open Source Software Archive 
Member "harminv-1.4.2/sines.c" (6 Feb 2023, 5767 Bytes) of package /linux/privat/harminv-1.4.2.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "sines.c" see the
Fossies "Dox" file reference documentation and the latest
Fossies "Diffs" side-by-side code changes report:
1.4.1_vs_1.4.2.
1 /* Copyright (C) 2017 Massachusetts Institute of Technology.
2 *
3 * This program is free software; you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation; either version 2 of the License, or
6 * (at your option) any later version.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program; if not, write to the Free Software
15 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 */
17
18 #include <stdlib.h>
19 #include <stdio.h>
20 #include <time.h>
21 #include <math.h>
22
23 #include "config.h"
24 #include "check.h"
25 #include "harminv-int.h"
26 #include "copyright.h"
27
28 #ifdef HAVE_UNISTD_H
29 #include <unistd.h>
30 #endif
31 #ifdef HAVE_GETOPT_H
32 #include <getopt.h>
33 #endif
34
35 int verbose = 0;
36
37 #define TWOPI 6.2831853071795864769252867665590057683943388
38 #define NPERIODS 10
39
40 void usage(FILE *f) {
41 fprintf(f,
42 "Usage: sines [options] <freq>...\n"
43 "Note that Im[freq] is the decay rate, or its inverse for -T.\n"
44 "Options: \n"
45 " -h : this help message\n"
46 " -V : print version number and copyright\n"
47 " -v : verbose output\n"
48 " -T : specify periods instead of frequencies\n"
49 " -R : real output\n"
50 " -r : random amplitudes\n"
51 " -N <a> : add white noise with amplitude <a>\n"
52 " -s <s> : use seed <s> for random #s (default <s>=time)\n"
53 " -n <n> : output <n> points (default %d * max period)\n"
54 " -t <dt> : time step <dt> (default 1.0)\n",
55 NPERIODS);
56 }
57
58 typedef struct {
59 double freq, decay, amplitude, phase;
60 } sinusoid;
61
62 #define MAX2(a, b) ((a) > (b) ? (a) : (b))
63
64 int main(int argc, char **argv) {
65 int c;
66 extern char *optarg;
67 extern int optind;
68 int iarg;
69 int specify_periods = 0;
70 int random_amplitudes = 0;
71 int real_output = 0;
72 sinusoid *sines = NULL;
73 int nsines = 0, nalloc = 0, n = 0;
74 double max_period = 0;
75 double dt = 1.0;
76 double noise = 0.0;
77 int i, is;
78
79 srand(time(NULL));
80
81 while ((c = getopt(argc, argv, "hVvTRrn:t:N:s:")) != -1)
82 switch (c) {
83 case 'h': usage(stdout); return EXIT_SUCCESS;
84 case 'V':
85 printf("sines " PACKAGE_VERSION " by Steven G. Johnson.\n"
86 "Test program for harminv.\n" COPYRIGHT);
87 return EXIT_SUCCESS;
88 case 'v': verbose = 1; break;
89 case 'T': specify_periods = 1; break;
90 case 'R': real_output = 1; break;
91 case 'r': random_amplitudes = 1; break;
92 case 'N': noise = atof(optarg); break;
93 case 's': srand(atoi(optarg)); break;
94 case 'n':
95 n = atoi(optarg);
96 if (n < 1) {
97 fprintf(stderr,
98 "sines: "
99 "invalid non-positive -n argument %d\n",
100 n);
101 return EXIT_FAILURE;
102 }
103 break;
104 case 't': dt = atof(optarg); break;
105 default:
106 fprintf(stderr, "sines: invalid argument -%c\n", c);
107 usage(stderr);
108 return EXIT_FAILURE;
109 }
110
111 if (optind == argc) { /* no parameters left */
112 fprintf(stderr, "sines: missing required frequency(ies)\n");
113 usage(stderr);
114 return EXIT_FAILURE;
115 }
116
117 for (iarg = optind; iarg < argc; ++iarg) {
118 sinusoid s = {0, 0, 0, 0};
119
120 if (sscanf(argv[iarg], "%lf+%lfi", &s.freq, &s.decay) < 1) {
121 fprintf(stderr, "sines: invalid argument \"%s\"\n", argv[iarg]);
122 return EXIT_FAILURE;
123 }
124 if (specify_periods) {
125 if (s.freq == 0) {
126 fprintf(stderr,
127 "sines: invalid argument \"%s\""
128 ": 0 not a valid period\n",
129 argv[iarg]);
130 return EXIT_FAILURE;
131 }
132 s.freq = 1 / s.freq;
133 if (s.decay != 0) s.decay = 1 / s.decay;
134 }
135
136 if (s.decay == 0 || fabs(1 / s.freq) > 1 / s.decay)
137 max_period = MAX2(max_period, fabs(1 / s.freq));
138 else
139 max_period = MAX2(max_period, 1 / s.decay);
140
141 if (random_amplitudes) {
142 s.amplitude = rand() * 1.0 / RAND_MAX;
143 s.phase = rand() * TWOPI / RAND_MAX - TWOPI / 2;
144 }
145 else {
146 s.amplitude = (iarg - optind + 1);
147 s.phase = (iarg - optind + 1) * TWOPI / (argc - optind) - TWOPI / 2;
148 }
149
150 if (verbose)
151 printf("# mode: frequency = %g (period %g), decay = %g (lifetime %g), amplitude = %g, phase "
152 "= %g\n",
153 s.freq, 1 / s.freq, s.decay, s.decay != 0 ? 1 / s.decay : 0, s.amplitude, s.phase);
154
155 if (nsines >= nalloc) {
156 nalloc = (nalloc + 1) * 2;
157 sines = (sinusoid *)realloc(sines, sizeof(sinusoid) * nalloc);
158 CHECK(sines, "out of memory");
159 }
160 sines[nsines++] = s;
161 }
162
163 if (n == 0 && max_period == 0) {
164 fprintf(stderr, "sines: must specify -n or non-zero period\n");
165 return EXIT_FAILURE;
166 }
167
168 if (n == 0) n = (int)(max_period / fabs(dt) + 0.5) * NPERIODS;
169
170 for (i = 0; i < n; ++i) {
171 cmplx output = 0;
172
173 for (is = 0; is < nsines; ++is)
174 output += sines[is].amplitude * cexp(I * (sines[is].phase - TWOPI * sines[is].freq * i * dt) -
175 sines[is].decay * i * dt);
176 output += noise * (rand() * 2.0 / RAND_MAX - 1);
177 if (real_output)
178 printf("%0.17e\n", creal(output));
179 else
180 printf("%0.17e%+0.17ei\n", creal(output), cimag(output));
181 }
182
183 free(sines);
184 return EXIT_SUCCESS;
185 }
186
187 #ifdef F77_DUMMY_MAIN
188 #ifdef __cplusplus
189 extern "C"
190 #endif
191 int
192 F77_DUMMY_MAIN() {
193 return 1;
194 }
195 #endif