"Fossies" - the Fresh Open Source Software Archive 
Member "harminv-1.4.2/harminv-main.c" (6 Feb 2023, 13012 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 "harminv-main.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 <ctype.h>
21 #include <math.h>
22
23 #include "harminv-int.h"
24 #include "check.h"
25 #include "copyright.h"
26
27 #ifdef HAVE_UNISTD_H
28 #include <unistd.h>
29 #endif
30 #ifdef HAVE_GETOPT_H
31 #include <getopt.h>
32 #endif
33
34 /* eat whitespace, including #... comments, from the file. Returns the
35 number of newlines read (so that a line count can be maintained). If
36 echo_comments != 0, then echo #... comments to stdout. Commas count
37 as whitespace, so that we can read comma-delimited text. */
38 static int eat_whitespace(FILE *f, int echo_comments) {
39 int c, newlines = 0;
40 do {
41 do {
42 c = getc(f);
43 newlines += c == '\n';
44 } while (isspace(c) || c == ',');
45
46 if (c == '#') { /* # begins comments that extend to the newline */
47 if (echo_comments) putc(c, stdout);
48 do {
49 c = getc(f);
50 if (echo_comments) {
51 if (c != EOF)
52 putc(c, stdout);
53 else /* terminate line if we hit EOF */
54 putc('\n', stdout);
55 }
56 newlines += c == '\n';
57 } while (c != EOF && c != '\n');
58 }
59 } while (isspace(c));
60
61 ungetc(c, f); /* put back the last character read */
62 newlines -= c == '\n';
63
64 return newlines;
65 }
66
67 static int eat_plus(FILE *f) {
68 int c = getc(f);
69 if (c != EOF && c != '+') ungetc(c, f);
70 return (c == '+' || c == '-');
71 }
72
73 static int eat_i(FILE *f) {
74 int c = getc(f);
75 if (c != EOF && tolower(c) != 'i') ungetc(c, f);
76 return (tolower(c) == 'i');
77 }
78
79 static cmplx *read_input_data(FILE *f, int *n, int verbose) {
80 cmplx *data = NULL;
81 int line = 1, n_alloc = 0;
82 *n = 0;
83
84 do {
85 double re = 0.0, im = 0.0;
86 int nread;
87 line += eat_whitespace(f, verbose);
88 nread = fscanf(f, "%lg", &re);
89 if (nread == 1 && eat_plus(f)) {
90 nread = fscanf(f, "%lg", &im);
91 if (nread == 1) nread = eat_i(f);
92 }
93 else
94 im = 0.0;
95 if (nread != EOF) {
96 if (nread < 1) {
97 fprintf(stderr, "harminv: invalid input on line %d.\n", line);
98 free(data);
99 *n = 0;
100 return NULL;
101 }
102 if (*n >= n_alloc) {
103 n_alloc = (n_alloc + 1) * 2;
104 data = (cmplx *)realloc(data, sizeof(cmplx) * n_alloc);
105 CHECK(data, "out of memory");
106 }
107 data[*n] = re + I * im;
108 ++*n;
109 }
110 } while (!feof(f));
111
112 data = (cmplx *)realloc(data, sizeof(cmplx) * *n);
113 return data;
114 }
115
116 static const double inf = HUGE_VAL;
117
118 #define DENSITY 0.0
119 #define NFMIN 100
120 #define NFMAX 300
121 #define ERR_THRESH 0.1
122 #define REL_ERR_THRESH inf
123 #define AMP_THRESH 0.0
124 #define REL_AMP_THRESH -1.0
125 #define Q_THRESH 10.0
126
127 static void usage(FILE *f) {
128 fprintf(f,
129 "Usage: harminv [options] <freq-min>-<freq-max>...\n"
130 "Options: \n"
131 " -h : this help message\n"
132 " -V : print version number and copyright\n"
133 " -v : verbose output\n"
134 " -T : specify periods instead of frequencies\n"
135 " -w : specify/output angular frequency, not frequency\n"
136 " -n : flip the sign of the frequency convention\n"
137 " -t <dt> : specify sampling interval dt [default: 1]\n"
138 " -d <d> : specify spectral density [default: %g]\n"
139 " -f <nf> : use at least <nf> basis functions [default: %d]\n"
140 " -s <sort> : sort by <sort> = freq/err/decay/amp [default: freq]\n"
141 " -F : discard frequencies outside of specified range\n"
142 " -a <a> : discard amplitudes < max * <a> [default: %e]\n"
143 " -A <A> : discard amplitudes < <A> [default: %g]\n"
144 " -e <e> : discard relative errors > min * <e> [default: %e]\n"
145 " -E <E> : discard relative errors > <E> [default: %e]\n"
146 " -Q <Q> : discard Q > <Q> [default: %g]\n",
147 DENSITY, NFMIN, AMP_THRESH, REL_AMP_THRESH, ERR_THRESH, REL_ERR_THRESH, Q_THRESH);
148 }
149
150 #define TWOPI 6.2831853071795864769252867665590057683943388
151
152 harminv_data hd;
153
154 enum { SORT_FREQUENCY, SORT_DECAY, SORT_ERROR, SORT_AMPLITUDE, SORT_Q } sortby = SORT_FREQUENCY;
155
156 static int cmp(double a, double b) { return a > b ? 1 : (a < b ? -1 : 0); }
157
158 static int compar(const void *a, const void *b) {
159 const int *ia = (const int *)a;
160 const int *ib = (const int *)b;
161 cmplx aa, ab;
162
163 switch (sortby) {
164 case SORT_FREQUENCY: return cmp(harminv_get_freq(hd, *ia), harminv_get_freq(hd, *ib));
165 case SORT_DECAY: return cmp(harminv_get_decay(hd, *ia), harminv_get_decay(hd, *ib));
166 case SORT_ERROR: return cmp(harminv_get_freq_error(hd, *ia), harminv_get_freq_error(hd, *ib));
167 case SORT_AMPLITUDE:
168 harminv_get_amplitude(&aa, hd, *ia);
169 harminv_get_amplitude(&ab, hd, *ia);
170 return cmp(cabs(aa), cabs(ab));
171 case SORT_Q:
172 return cmp(harminv_get_freq(hd, *ia) / harminv_get_decay(hd, *ia),
173 harminv_get_freq(hd, *ib) / harminv_get_decay(hd, *ib));
174 }
175 return 0;
176 }
177
178 typedef struct {
179 int verbose;
180 double fmin, fmax;
181 int only_f_inrange;
182 double err_thresh, rel_err_thresh, amp_thresh, rel_amp_thresh, Q_thresh;
183 double min_err, max_amp;
184 int num_ok;
185 } mode_ok_data;
186
187 static int mode_ok(harminv_data d, int k, void *ok_d_) {
188 mode_ok_data *ok_d = (mode_ok_data *)ok_d_;
189 double errk, ampk, f;
190 int ok;
191 cmplx aa;
192
193 if (k == -1) { /* initialize */
194 int i;
195 ok_d->num_ok = 0;
196 if (!harminv_get_num_freqs(d)) return 0;
197 ok_d->min_err = harminv_get_freq_error(d, 0);
198 ;
199 harminv_get_amplitude(&aa, d, 0);
200 ok_d->max_amp = cabs(aa);
201 for (i = 1; i < harminv_get_num_freqs(d); ++i) {
202 double err, amp;
203 if ((err = harminv_get_freq_error(d, i)) < ok_d->min_err) ok_d->min_err = err;
204 harminv_get_amplitude(&aa, d, i);
205 if ((amp = cabs(aa)) > ok_d->max_amp) ok_d->max_amp = amp;
206 }
207 return 0;
208 }
209 else if (k == -2) { /* finish */
210 if (ok_d->verbose && harminv_get_num_freqs(d))
211 printf("# harminv: %d/%d modes are ok: "
212 "errs <= %e and %e * %e\n, "
213 "amps >= %g, %e * %g, "
214 "|Q| >= %g\n",
215 ok_d->num_ok, harminv_get_num_freqs(d), ok_d->err_thresh, ok_d->rel_err_thresh,
216 ok_d->min_err, ok_d->amp_thresh, ok_d->rel_amp_thresh, ok_d->max_amp, ok_d->Q_thresh);
217 return 0;
218 }
219
220 f = fabs(harminv_get_freq(d, k));
221 errk = harminv_get_freq_error(d, k);
222 harminv_get_amplitude(&aa, d, k);
223 ampk = cabs(aa);
224
225 ok = ((!ok_d->only_f_inrange || (f >= ok_d->fmin && f <= ok_d->fmax)) &&
226 errk <= ok_d->err_thresh &&
227 (!isfinite(ok_d->rel_err_thresh) || errk <= ok_d->min_err * ok_d->rel_err_thresh) &&
228 ampk >= ok_d->amp_thresh && ampk >= ok_d->rel_amp_thresh * ok_d->max_amp &&
229 fabs(harminv_get_Q(d, k)) >= ok_d->Q_thresh);
230
231 ok_d->num_ok += ok;
232
233 return ok;
234 }
235
236 #define SOLVE_ONCE_ONLY 0 /* 1 to use harminv_solve_once */
237 #define SOLVE_OK_ONLY 0 /* 1 for experimental solver */
238
239 int main(int argc, char **argv) {
240 int verbose = 0;
241 int c;
242 extern char *optarg;
243 extern int optind;
244 int specify_periods = 0;
245 int specify_omega = 0;
246 int negate_omega = 0;
247 double dt = 1.0;
248 mode_ok_data ok_d;
249 int n, nf, nfmin = NFMIN;
250 double density = DENSITY;
251 int iarg;
252 cmplx *data;
253
254 ok_d.only_f_inrange = 0;
255 ok_d.err_thresh = ERR_THRESH;
256 ok_d.rel_err_thresh = REL_ERR_THRESH;
257 ok_d.amp_thresh = AMP_THRESH;
258 ok_d.rel_amp_thresh = REL_AMP_THRESH;
259 ok_d.Q_thresh = Q_THRESH;
260
261 while ((c = getopt(argc, argv, "hvVTFwnt:d:f:s:e:E:a:A:Q:")) != -1)
262 switch (c) {
263 case 'h': usage(stdout); return EXIT_SUCCESS;
264 case 'V':
265 printf("harminv " PACKAGE_VERSION " by Steven G. Johnson\n" COPYRIGHT);
266 return EXIT_SUCCESS;
267 case 'v': verbose = 1; break;
268 case 'T': specify_periods = 1; break;
269 case 'w': specify_omega = 1; break;
270 case 'n': negate_omega = 1; break;
271 case 'F': ok_d.only_f_inrange = 1; break;
272 case 'a': ok_d.rel_amp_thresh = atof(optarg); break;
273 case 'A': ok_d.amp_thresh = atof(optarg); break;
274 case 'E': ok_d.err_thresh = atof(optarg); break;
275 case 'e': ok_d.rel_err_thresh = atof(optarg); break;
276 case 'Q': ok_d.Q_thresh = atof(optarg); break;
277 case 't': dt = atof(optarg); break;
278 case 'f': nfmin = atoi(optarg); break;
279 case 'd':
280 density = atof(optarg);
281 if (density < 0) {
282 fprintf(stderr, "harminv: -d argument must be >= 0\n");
283 return EXIT_FAILURE;
284 }
285 break;
286 case 's':
287 switch (tolower(optarg[0])) {
288 case 'f': sortby = SORT_FREQUENCY; break;
289 case 'd': sortby = SORT_DECAY; break;
290 case 'e': sortby = SORT_ERROR; break;
291 case 'a': sortby = SORT_AMPLITUDE; break;
292 case 'q': sortby = SORT_Q; break;
293 default:
294 fprintf(stderr, "harminv: invalid sort type -s %c\n", tolower(optarg[0]));
295 usage(stderr);
296 return EXIT_FAILURE;
297 }
298 break;
299 default:
300 fprintf(stderr, "harminv: invalid argument -%c\n", c);
301 usage(stderr);
302 return EXIT_FAILURE;
303 }
304 if (optind == argc) { /* no parameters left */
305 fprintf(stderr, "harminv: missing required frequency range(s)\n");
306 usage(stderr);
307 return EXIT_FAILURE;
308 }
309
310 /* harminv requires nf > 1 */
311 if (nfmin < 2) nfmin = 2;
312
313 data = read_input_data(stdin, &n, verbose);
314
315 if (n < 1) {
316 fprintf(stderr, "harminv: no data read\n");
317 return EXIT_FAILURE;
318 }
319
320 if (verbose) printf("# harminv: %d inputs, dt = %g\n", n, dt);
321
322 printf("frequency, decay constant, Q, amplitude, phase, error\n");
323
324 ok_d.verbose = verbose;
325
326 for (iarg = optind; iarg < argc; ++iarg) {
327 double fmin, fmax;
328 int i;
329 int *isort = NULL;
330
331 if (sscanf(argv[iarg], "%lf-%lf", &fmin, &fmax) != 2) {
332 fprintf(stderr, "harminv: invalid argument \"%s\"\n", argv[iarg]);
333 return EXIT_FAILURE;
334 }
335 if (specify_periods) {
336 if (fmin == 0 || fmax == 0) {
337 fprintf(stderr,
338 "harminv: invalid argument \"%s\""
339 ": 0 not a valid period\n",
340 argv[iarg]);
341 return EXIT_FAILURE;
342 }
343 fmin = 1 / fmin;
344 fmax = 1 / fmax;
345 }
346 if (specify_omega) {
347 fmin /= TWOPI;
348 fmax /= TWOPI;
349 }
350 if (negate_omega) dt *= -1;
351 if ((fmin > fmax && dt > 0) || (fmin < fmax && dt < 0)) {
352 double dummy = fmin;
353 fmin = fmax;
354 fmax = dummy;
355 }
356 if (verbose) printf("# searching frequency range %g - %g\n", fmin, fmax);
357
358 ok_d.fmin = fabs(fmin * dt);
359 ok_d.fmax = fabs(fmax * dt);
360 if (ok_d.fmin > ok_d.fmax) {
361 double dummy = ok_d.fmin;
362 ok_d.fmin = ok_d.fmax;
363 ok_d.fmax = dummy;
364 }
365
366 nf = (fmax - fmin) * dt * n * density;
367 if (nf > NFMAX) nf = NFMAX;
368 if (nf < nfmin) nf = nfmin;
369 if (verbose)
370 printf("# using %d spectral basis functions, density %g\n", nf,
371 nf / ((fmax - fmin) * dt * n));
372
373 hd = harminv_data_create(n, data, fmin * dt, fmax * dt, nf);
374
375 #if SOLVE_OK_ONLY
376 harminv_solve_ok_modes(hd, mode_ok, &ok_d);
377 #elif SOLVE_ONCE_ONLY
378 harminv_solve_once(hd);
379 #else
380 harminv_solve(hd);
381 #endif
382
383 mode_ok(hd, -1, &ok_d); /* initialize ok_d */
384
385 CHK_MALLOC(isort, int, harminv_get_num_freqs(hd));
386 for (i = 0; i < harminv_get_num_freqs(hd); ++i)
387 isort[i] = i;
388 qsort(isort, harminv_get_num_freqs(hd), sizeof(int), compar);
389
390 for (i = 0; i < harminv_get_num_freqs(hd); ++i) {
391 double freq, decay, err;
392 cmplx amp;
393 int j = isort[i];
394 #if SOLVE_OK_ONLY
395 CHECK(mode_ok(hd, j, &ok_d), "bug: invalid mode");
396 #else
397 if (!mode_ok(hd, j, &ok_d)) continue;
398 #endif
399 freq = harminv_get_freq(hd, j) / dt;
400 decay = harminv_get_decay(hd, j) / fabs(dt);
401 harminv_get_amplitude(&, hd, j);
402 err = harminv_get_freq_error(hd, j);
403 printf("%g, %e, %g, %g, %g, %e\n", freq * (specify_omega ? TWOPI : 1.0), decay,
404 harminv_get_Q(hd, j), cabs(amp), carg(amp) * (dt < 0 ? -1 : 1), err);
405 }
406
407 #if !SOLVE_OK_ONLY
408 mode_ok(hd, -2, &ok_d);
409 #endif
410
411 harminv_data_destroy(hd);
412 }
413
414 free(data);
415 return EXIT_SUCCESS;
416 }
417
418 #ifdef F77_DUMMY_MAIN
419 #ifdef __cplusplus
420 extern "C"
421 #endif
422 int
423 F77_DUMMY_MAIN() {
424 return 1;
425 }
426 #endif