1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: grdfft.c is a program to do various operations on
19 * grid files in the frequency domain.
20 *
21 * Author: W.H.F. Smith
22 * Date: 1-JAN-2010
23 * Version: 6 API
24 * Note: PW: As of 2/14/2013 the various setup and init functions for FFT use
25 * have been generalized and made available GMT-wide via new functions
26 * in gmt_fft.c, called GMT_fft_*.
27 */
28
29 #include "gmt_dev.h"
30
31 #define THIS_MODULE_CLASSIC_NAME "grdfft"
32 #define THIS_MODULE_MODERN_NAME "grdfft"
33 #define THIS_MODULE_LIB "core"
34 #define THIS_MODULE_PURPOSE "Mathematical operations on grids in the spectral domain"
35 #define THIS_MODULE_KEYS "<G{+,GG},GDE"
36 #define THIS_MODULE_NEEDS ""
37 #define THIS_MODULE_OPTIONS "-Vfh" GMT_OPT("T")
38
39 #define GMT_FFT_DIM 2 /* Dimension of FFT needed */
40
41 #ifdef DEBUG
42 /* For debugging -E; running this in debug and setting it to true will also output the number of estimates per radial k */
43 static bool show_n = false;
44 #endif
45
46 struct GRDFFT_CTRL {
47 unsigned int n_op_count, n_par;
48 int *operation;
49 double *par;
50
51 struct GRDFFT_In {
52 bool active;
53 unsigned int n_grids; /* 1 or 2 */
54 char *file[2];
55 } In;
56 struct GRDFFT_A { /* -A<azimuth> */
57 bool active;
58 double value;
59 } A;
60 struct GRDFFT_C { /* -C<zlevel> */
61 bool active;
62 double value;
63 } C;
64 struct GRDFFT_D { /* -D[<scale>|g] */
65 bool active;
66 double value;
67 } D;
68 struct GRDFFT_E { /* -E[r|x|y][w[k]][n] */
69 bool active;
70 bool give_wavelength;
71 bool normalize;
72 bool km;
73 int mode; /*-1/0/+1 */
74 } E;
75 struct GRDFFT_F { /* -F[r|x|y]<lc>/<lp>/<hp>/<hc> or -F[r|x|y]<lo>/<hi> */
76 bool active;
77 double lc, lp, hp, hc;
78 } F;
79 struct GRDFFT_G { /* -G<outfile> */
80 bool active;
81 char *file;
82 } G;
83 struct GRDFFT_I { /* -I[<scale>|g] */
84 bool active;
85 double value;
86 } I;
87 struct GRDFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
88 bool active;
89 struct GMT_FFT_INFO *info;
90 } N;
91 struct GRDFFT_Q { /* -Q */
92 bool active;
93 } Q;
94 struct GRDFFT_S { /* -S<scale> */
95 bool active;
96 double scale;
97 } S;
98 /* Now in gravfft in potential supplement; left for backwards compatibility */
99 struct GRDFFT_T { /* -T<te/rl/rm/rw/ri> */
100 bool active;
101 double te, rhol, rhom, rhow, rhoi;
102 } T;
103 };
104
105 enum Grdfft_operators {
106 GRDFFT_UP_DOWN_CONTINUE = 0,
107 GRDFFT_AZIMUTHAL_DERIVATIVE,
108 GRDFFT_DIFFERENTIATE ,
109 GRDFFT_INTEGRATE ,
110 GRDFFT_FILTER_EXP ,
111 GRDFFT_FILTER_BW ,
112 GRDFFT_FILTER_COS ,
113 GRDFFT_SPECTRUM ,
114 GRDFFT_ISOSTASY ,
115 GRDFFT_NOTHING };
116
117 #define MGAL_AT_45 980619.9203 /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
118
119 struct F_INFO {
120 double lc[3]; /* Low-cut frequency for r, x, and y */
121 double lp[3]; /* Low-pass frequency for r, x, and y */
122 double hp[3]; /* High-pass frequency for r, x, and y */
123 double hc[3]; /* High-cut frequency for r, x, and y */
124 double ltaper[3]; /* Low taper width for r, x, and y */
125 double htaper[3]; /* High taper width for r, x, and y */
126 double llambda[3]; /* Low full-wavelength where Gauss amp = 0.5 for r, x, and y */
127 double hlambda[3]; /* High full-wavelength where Gauss amp = 0.5 for r, x, and y */
128 double bw_order; /* Order, N, of Butterworth filter */
129 double (*filter) (struct F_INFO *, double, int); /* Points to the correct filter function */
130
131 bool set_already; /* true if we already filled in the structure */
132 unsigned int k_type; /* Which of the r, x, or y wavenumbers we need */
133 unsigned int kind; /* GRDFFT_FILTER_EXP, GRDFFT_FILTER_BW, GRDFFT_FILTER_COS */
134 unsigned int arg; /* 0 = Gaussian, 1 = Butterworth, 2 = cosine taper, */
135 };
136
New_Ctrl(struct GMT_CTRL * GMT)137 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
138 struct GRDFFT_CTRL *C = NULL;
139
140 C = gmt_M_memory (GMT, NULL, 1, struct GRDFFT_CTRL);
141
142 /* Initialize values whose defaults are not 0/false/NULL */
143
144 C->S.scale = 1.0;
145 return (C);
146 }
147
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDFFT_CTRL * C)148 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *C) { /* Deallocate control structure */
149 if (!C) return;
150 gmt_M_free (GMT, C->operation);
151 gmt_M_free (GMT, C->par);
152 gmt_M_str_free (C->In.file[GMT_IN]);
153 gmt_M_str_free (C->In.file[GMT_OUT]);
154 gmt_M_str_free (C->G.file);
155 gmt_M_free (GMT, C->N.info);
156 gmt_M_free (GMT, C);
157 }
158
grdfft_do_differentiate(struct GMT_GRID * Grid,double * par,struct GMT_FFT_WAVENUMBER * K)159 GMT_LOCAL unsigned int grdfft_do_differentiate (struct GMT_GRID *Grid, double *par, struct GMT_FFT_WAVENUMBER *K) {
160 uint64_t k;
161 double scale, fact;
162 gmt_grdfloat *datac = Grid->data; /* Shorthand */
163
164 /* Differentiate in frequency domain by multiplying by kr [scale optional] */
165
166 scale = (*par != 0.0) ? *par : 1.0;
167 datac[0] = datac[1] = 0.0f; /* Derivative of the mean is zero */
168 for (k = 2; k < Grid->header->size; k += 2) {
169 fact = scale * gmt_fft_get_wave (k, K);
170 datac[k] *= (gmt_grdfloat)fact;
171 datac[k+1] *= (gmt_grdfloat)fact;
172 }
173 return (1); /* Number of parameters used */
174 }
175
grdfft_do_integrate(struct GMT_GRID * Grid,double * par,struct GMT_FFT_WAVENUMBER * K)176 GMT_LOCAL unsigned int grdfft_do_integrate (struct GMT_GRID *Grid, double *par, struct GMT_FFT_WAVENUMBER *K) {
177 /* Integrate in frequency domain by dividing by kr [scale optional] */
178 uint64_t k;
179 double fact, scale;
180 gmt_grdfloat *datac = Grid->data; /* Shorthand */
181
182 scale = (*par != 0.0) ? *par : 1.0;
183 datac[0] = datac[1] = 0.0f;
184 for (k = 2; k < Grid->header->size; k += 2) {
185 fact = 1.0 / (scale * gmt_fft_get_wave (k, K));
186 datac[k] *= (gmt_grdfloat)fact;
187 datac[k+1] *= (gmt_grdfloat)fact;
188 }
189 return (1); /* Number of parameters used */
190 }
191
grdfft_do_continuation(struct GMT_GRID * Grid,double * zlevel,struct GMT_FFT_WAVENUMBER * K)192 GMT_LOCAL unsigned int grdfft_do_continuation (struct GMT_GRID *Grid, double *zlevel, struct GMT_FFT_WAVENUMBER *K) {
193 uint64_t k;
194 gmt_grdfloat tmp, *datac = Grid->data; /* Shorthand */
195
196 /* If z is positive, the field will be upward continued using exp[- k z]. */
197
198 for (k = 2; k < Grid->header->size; k += 2) {
199 tmp = (gmt_grdfloat)exp (-(*zlevel) * gmt_fft_get_wave (k, K));
200 datac[k] *= tmp;
201 datac[k+1] *= tmp;
202 }
203 return (1); /* Number of parameters used */
204 }
205
grdfft_do_azimuthal_derivative(struct GMT_GRID * Grid,double * azim,struct GMT_FFT_WAVENUMBER * K)206 GMT_LOCAL unsigned int grdfft_do_azimuthal_derivative (struct GMT_GRID *Grid, double *azim, struct GMT_FFT_WAVENUMBER *K) {
207 uint64_t k;
208 gmt_grdfloat tempr, tempi, fact, *datac = Grid->data; /* Shorthand */
209 double cos_azim, sin_azim;
210
211 sincosd (*azim, &sin_azim, &cos_azim);
212
213 datac[0] = datac[1] = 0.0f;
214 for (k = 2; k < Grid->header->size; k += 2) {
215 fact = (gmt_grdfloat)(sin_azim * gmt_fft_any_wave (k, GMT_FFT_K_IS_KX, K) + cos_azim * gmt_fft_any_wave (k, GMT_FFT_K_IS_KY, K));
216 tempr = -(datac[k+1] * fact);
217 tempi = (datac[k] * fact);
218 datac[k] = tempr;
219 datac[k+1] = tempi;
220 }
221 return (1); /* Number of parameters used */
222 }
223
224 /* Now obsolete but left for backwards compatibility. Users are encouraged to use gravfft */
225 #define POISSONS_RATIO 0.25
226 #define YOUNGS_MODULUS 1.0e11 /* Pascal = Nt/m**2 */
227 #define NORMAL_GRAVITY 9.806199203 /* m/s**2 */
228
grdfft_do_isostasy(struct GMT_GRID * Grid,struct GRDFFT_CTRL * Ctrl,double * par,struct GMT_FFT_WAVENUMBER * K)229 GMT_LOCAL unsigned int grdfft_do_isostasy (struct GMT_GRID *Grid, struct GRDFFT_CTRL *Ctrl, double *par, struct GMT_FFT_WAVENUMBER *K) {
230 /* Do the isostatic response function convolution in the Freq domain.
231 All units assumed to be in SI (that is kx, ky, modk wavenumbers in m**-1,
232 densities in kg/m**3, Te in m, etc.
233 rw, the water density, is used to set the Airy ratio and the restoring
234 force on the plate (rm - ri)*gravity if ri = rw; so use zero for topo in air. */
235 uint64_t k;
236 double airy_ratio, rigidity_d, d_over_restoring_force, mk, k2, k4, transfer_fn;
237
238 double te; /* Elastic thickness, SI units (m) */
239 double rl; /* Load density, SI units */
240 double rm; /* Mantle density, SI units */
241 double rw; /* Water density, SI units */
242 double ri; /* Infill density, SI units */
243
244 gmt_grdfloat *datac = Grid->data; /* Shorthand */
245
246 te = par[0]; rl = par[1]; rm = par[2]; rw = par[3]; ri = par[4];
247 airy_ratio = -(rl - rw)/(rm - ri);
248
249 if (te == 0.0) { /* Airy isostasy; scale variable Ctrl->S.scale and return */
250 Ctrl->S.scale *= airy_ratio;
251 return (5); /* Number of parameters used */
252 }
253
254 rigidity_d = (YOUNGS_MODULUS * pow (te, 3.0)) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
255 d_over_restoring_force = rigidity_d / ((rm - ri) * NORMAL_GRAVITY);
256
257 for (k = 0; k < Grid->header->size; k += 2) {
258 mk = gmt_fft_get_wave (k, K);
259 k2 = mk * mk;
260 k4 = k2 * k2;
261 transfer_fn = airy_ratio / ((d_over_restoring_force * k4) + 1.0);
262 datac[k] *= (gmt_grdfloat)transfer_fn;
263 datac[k+1] *= (gmt_grdfloat)transfer_fn;
264 }
265 return (5); /* Number of parameters used */
266 }
267
268 #ifndef M_LN2
269 #define M_LN2 0.69314718055994530942 /* log_e 2 */
270 #endif
271
grdfft_gauss_weight(struct F_INFO * f_info,double freq,int j)272 GMT_LOCAL double grdfft_gauss_weight (struct F_INFO *f_info, double freq, int j) {
273 double hi, lo;
274 lo = (f_info->llambda[j] == -1.0) ? 0.0 : exp (-M_LN2 * pow (freq * f_info->llambda[j], 2.0)); /* Low-pass part */
275 hi = (f_info->hlambda[j] == -1.0) ? 1.0 : exp (-M_LN2 * pow (freq * f_info->hlambda[j], 2.0)); /* Hi-pass given by its complementary low-pass */
276 return (hi - lo);
277 }
278
grdfft_bw_weight(struct F_INFO * f_info,double freq,int j)279 GMT_LOCAL double grdfft_bw_weight (struct F_INFO *f_info, double freq, int j) { /* Butterworth filter */
280 double hi, lo;
281 lo = (f_info->llambda[j] == -1.0) ? 0.0 : sqrt (1.0 / (1.0 + pow (freq * f_info->llambda[j], f_info->bw_order))); /* Low-pass part */
282 hi = (f_info->hlambda[j] == -1.0) ? 1.0 : sqrt (1.0 / (1.0 + pow (freq * f_info->hlambda[j], f_info->bw_order))); /* Hi-pass given by its complementary low-pass */
283 return (hi - lo);
284 }
285
grdfft_cosine_weight_grdfft(struct F_INFO * f_info,double freq,int j)286 GMT_LOCAL double grdfft_cosine_weight_grdfft (struct F_INFO *f_info, double freq, int j) {
287 if (freq <= f_info->lc[j] || freq >= f_info->hc[j]) return(0.0); /* In fully cut range. Weight is zero. */
288 if (freq > f_info->lc[j] && freq < f_info->lp[j]) return (0.5 * (1.0 + cos (M_PI * (freq - f_info->lp[j]) * f_info->ltaper[j])));
289 if (freq > f_info->hp[j] && freq < f_info->hc[j]) return (0.5 * (1.0 + cos (M_PI * (freq - f_info->hp[j]) * f_info->htaper[j])));
290 return (1.0); /* Freq is in the fully passed range, so weight is multiplied by 1.0 */
291 }
292
grdfft_get_filter_weight(uint64_t k,struct F_INFO * f_info,struct GMT_FFT_WAVENUMBER * K)293 GMT_LOCAL double grdfft_get_filter_weight (uint64_t k, struct F_INFO *f_info, struct GMT_FFT_WAVENUMBER *K) {
294 double freq, return_value;
295
296 freq = gmt_fft_any_wave (k, f_info->k_type, K);
297 return_value = f_info->filter (f_info, freq, f_info->k_type);
298
299 return (return_value);
300 }
301
grdfft_do_filter(struct GMT_GRID * Grid,struct F_INFO * f_info,struct GMT_FFT_WAVENUMBER * K)302 GMT_LOCAL void grdfft_do_filter (struct GMT_GRID *Grid, struct F_INFO *f_info, struct GMT_FFT_WAVENUMBER *K) {
303 uint64_t k;
304 gmt_grdfloat weight, *datac = Grid->data; /* Shorthand */
305
306 for (k = 0; k < Grid->header->size; k += 2) {
307 weight = (gmt_grdfloat) grdfft_get_filter_weight (k, f_info, K);
308 datac[k] *= weight;
309 datac[k+1] *= weight;
310 }
311 }
312
grdfft_do_spectrum(struct GMT_CTRL * GMT,struct GMT_GRID * GridX,struct GMT_GRID * GridY,double * par,bool give_wavelength,bool km,bool normalize,char * file,struct GMT_FFT_WAVENUMBER * K)313 GMT_LOCAL int grdfft_do_spectrum (struct GMT_CTRL *GMT, struct GMT_GRID *GridX, struct GMT_GRID *GridY, double *par, bool give_wavelength, bool km, bool normalize, char *file, struct GMT_FFT_WAVENUMBER *K) {
314 /* Compute [cross-]spectral estimates from the two grids X and Y and return frequency f and 8 quantities:
315 * Xpower[f], Ypower[f], coherent power[f], noise power[f], phase[f], admittance[f], gain[f], coherency[f].
316 * Each quantity comes with its own 1-std dev error estimate, hence output is 17 columns. If GridY == NULL
317 * then just XPower[f] and its 1-std dev error estimate are returned, hence just 3 columns.
318 * Equations based on spectrum1d.c */
319
320 uint64_t dim[GMT_DIM_SIZE] = {1, 1, 0, 0}; /* One table and one segment, with either 1 + 1*2 = 3 or 1 + 8*2 = 17 columns and yet unknown rows */
321 uint64_t k, nk, ifreq, *nused = NULL;
322 unsigned int col;
323 gmt_grdfloat *X = GridX->data, *Y = NULL; /* Short-hands */
324 double delta_k, r_delta_k, freq, coh_k, sq_norm, powfactor, k_pow_factor, tmp, eps_pow;
325 double *X_pow = NULL, *Y_pow = NULL, *co_spec = NULL, *quad_spec = NULL;
326 struct GMT_DATASET *D = NULL;
327 struct GMT_DATASEGMENT *S = NULL;
328
329 dim[GMT_COL] = (GridY) ? 17 : 3; /* Either 3 or 17 estimates */
330 #ifdef DEBUG
331 if (show_n) dim[GMT_COL]++; /* Also write out nused[k] so add one more column */
332 #endif
333 if (*par > 0.0) { /* X spectrum desired */
334 delta_k = K->delta_kx; nk = K->nx2 / 2;
335 gmt_fft_set_wave (GMT, GMT_FFT_K_IS_KX, K);
336 }
337 else if (*par < 0.0) { /* Y spectrum desired */
338 delta_k = K->delta_ky; nk = K->ny2 / 2;
339 gmt_fft_set_wave (GMT, GMT_FFT_K_IS_KY, K);
340 }
341 else { /* R spectrum desired */
342 if (K->delta_kx < K->delta_ky) {
343 delta_k = K->delta_kx; nk = K->nx2 / 2;
344 }
345 else {
346 delta_k = K->delta_ky; nk = K->ny2 / 2;
347 }
348 gmt_fft_set_wave (GMT, GMT_FFT_K_IS_KR, K);
349 }
350
351 /* Get arrays for summing stuff */
352 X_pow = gmt_M_memory (GMT, NULL, nk, double );
353 nused = gmt_M_memory (GMT, NULL, nk, uint64_t);
354 if (GridY) { /* For cross-spectral estimates */
355 Y = GridY->data; /* Shorthand for Y data */
356 Y_pow = gmt_M_memory (GMT, NULL, nk, double);
357 co_spec = gmt_M_memory (GMT, NULL, nk, double);
358 quad_spec = gmt_M_memory (GMT, NULL, nk, double);
359 }
360
361 /* Loop over it all, summing and storing, checking range for r */
362
363 r_delta_k = 1.0 / delta_k;
364
365 for (k = 2; k < GridX->header->size; k += 2) {
366 freq = gmt_fft_get_wave (k, K);
367 ifreq = lrint (fabs (freq) * r_delta_k); /* Smallest value returned might be 0 when doing r spectrum*/
368 if (ifreq > 0) --ifreq;
369 if (ifreq >= nk) continue; /* Might happen when doing r spectrum */
370 X_pow[ifreq] += (X[k] * X[k] + X[k+1] * X[k+1]); /* X x X* = Power of grid X */
371 nused[ifreq]++;
372 if (GridY) { /* For cross-spectral estimates */
373 Y_pow[ifreq] += (Y[k] * Y[k] + Y[k+1] * Y[k+1]); /* Y x Y* = Power of grid Y */
374 co_spec[ifreq] += (Y[k] * X[k] + Y[k+1] * X[k+1]); /* Real part of Y x X* */
375 quad_spec[ifreq] += (X[k+1] * Y[k] - Y[k+1] * X[k]); /* Imag part of Y x X* */
376 }
377 }
378
379 gmt_set_cartesian (GMT, GMT_OUT); /* To counter-act any -fg setting */
380
381 delta_k /= (2.0 * M_PI); /* Write out frequency, not wavenumber */
382 powfactor = 4.0 / pow ((double)GridX->header->size, 2.0); /* Squared normalization of FFT */
383 dim[GMT_ROW] = nk;
384 if ((D = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_NONE, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
385 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create a data set for spectral estimates\n");
386 gmt_M_free (GMT, X_pow);
387 gmt_M_free (GMT, Y_pow);
388 gmt_M_free (GMT, nused);
389 gmt_M_free (GMT, co_spec);
390 gmt_M_free (GMT, quad_spec);
391 return (GMT->parent->error);
392 }
393 S = D->table[0]->segment[0]; /* Only one table with one segment here, with 17 cols and nk rows */
394 S->n_rows = nk;
395 if (give_wavelength && km) delta_k *= 1000.0; /* Wanted distances measured in km */
396
397 k_pow_factor = powfactor; /* Standard normalization */
398 for (k = 0; k < nk; k++) {
399 eps_pow = 1.0 / sqrt ((double)nused[k]); /* Multiplicative error bars for power spectra */
400 freq = (k + 1) * delta_k;
401 if (give_wavelength) freq = 1.0 / freq;
402
403 if (normalize) k_pow_factor = powfactor / (double)nused[k];
404 X_pow[k] *= k_pow_factor;
405
406 col = 0;
407 /* Col 0 is the frequency (or wavelength) */
408 S->data[col++][k] = freq;
409 /* Cols 1-2 are xpower and std.err estimate */
410 S->data[col++][k] = X_pow[k];
411 S->data[col++][k] = X_pow[k] * eps_pow;
412
413 if (!GridY) { /* Nothing more to do (except add nused[k] if true and debug) */
414 #ifdef DEBUG
415 if (show_n) S->data[col][k] = (double)nused[k];
416 #endif
417 continue;
418 }
419
420 Y_pow[k] *= k_pow_factor;
421 co_spec[k] *= k_pow_factor;
422 quad_spec[k] *= k_pow_factor;
423 /* Compute coherence first since it is needed by many of the other estimates */
424 coh_k = (co_spec[k] * co_spec[k] + quad_spec[k] * quad_spec[k]) / (X_pow[k] * Y_pow[k]);
425 sq_norm = sqrt ((1.0 - coh_k) / (2.0 * coh_k)); /* Save repetitive expression further down */
426 /* Cols 3-4 are ypower and std.err estimate */
427 S->data[col++][k] = Y_pow[k];
428 S->data[col++][k] = Y_pow[k] * eps_pow;
429 /* Cols 5-6 are coherent power and std.err estimate */
430 S->data[col++][k] = tmp = Y_pow[k] * coh_k;
431 S->data[col++][k] = tmp * eps_pow * sqrt ((2.0 - coh_k) / coh_k);
432 /* Cols 7-8 are noise power and std.err estimate */
433 S->data[col++][k] = tmp = Y_pow[k] * (1.0 - coh_k);
434 S->data[col++][k] = tmp * eps_pow;
435 /* Cols 9-10 are phase and std.err estimate */
436 S->data[col++][k] = tmp = d_atan2 (quad_spec[k], co_spec[k]);
437 S->data[col++][k] = tmp * eps_pow * sq_norm;
438 /* Cols 11-12 are admittance and std.err estimate */
439 S->data[col++][k] = tmp = co_spec[k] / X_pow[k];
440 S->data[col++][k] = tmp * eps_pow * fabs (sq_norm);
441 /* Cols 13-14 are gain and std.err estimate */
442 S->data[col++][k] = tmp = sqrt (quad_spec[k]) / X_pow[k];
443 S->data[col++][k] = tmp * eps_pow * sq_norm;
444 /* Cols 15-16 are coherency and std.err estimate */
445 S->data[col++][k] = coh_k;
446 S->data[col++][k] = coh_k * eps_pow * (1.0 - coh_k) * sqrt (2.0 / coh_k);
447 #ifdef DEBUG
448 if (show_n) S->data[col][k] = (double)nused[k];
449 #endif
450 }
451
452 if (GMT->common.h.add_colnames) {
453 char header[GMT_BUFSIZ] = {""}, *name[2] = {"freq", "wlength"};
454 if (GridY) { /* Long header record - number in [] is GMT column; useful for -i option */
455 sprintf (header, "%s[0]\txpow[1]\tstd_xpow[2]\typow[3]\tstd_ypow[4]\tcpow[5]\tstd_cpow[6]\tnpow[7]\tstd_npow[8]\t" \
456 "phase[9]\tstd_phase[10]\tadm[11]\tstd_ad[12]\tgain[13]\tstd_gain[14]\tcoh[15]\tstd_coh[16]", name[give_wavelength]);
457 }
458 else
459 sprintf (header, "%s[0]\tpow[1]\tstd_pow[2]", name[give_wavelength]);
460 if (GMT_Set_Comment (GMT->parent, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, D)) return (GMT->parent->error);
461 }
462
463 if (GMT_Write_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, file, D) != GMT_NOERROR) {
464 return (GMT->parent->error);
465 }
466 gmt_M_free (GMT, X_pow);
467 gmt_M_free (GMT, nused);
468 if (GridY) {
469 gmt_M_free (GMT, Y_pow);
470 gmt_M_free (GMT, co_spec);
471 gmt_M_free (GMT, quad_spec);
472 }
473
474 return (1); /* Number of parameters used */
475 }
476
grdfft_parse_f_string(struct GMT_CTRL * GMT,struct F_INFO * f_info,char * c)477 GMT_LOCAL bool grdfft_parse_f_string (struct GMT_CTRL *GMT, struct F_INFO *f_info, char *c) {
478 unsigned int i, j, n_tokens, pos;
479 bool descending;
480 double fourvals[4];
481 char line[GMT_LEN256] = {""}, p[GMT_LEN256] = {""};
482
483 /* Syntax is either -F[r|x|y]lc/hc/lp/hp (Cosine taper), -F[r|x|y]lo/hi (Gaussian), or -F[r|x|y]lo/hi/order (Butterworth) */
484
485 strncpy (line, c, GMT_LEN256-1);
486 i = 0;
487 f_info->k_type = GMT_FFT_K_IS_KR; /* j is Filter type: r=2, x=0, y=1 [r] */
488
489 if (line[i] == 'r') {
490 i++; f_info->k_type = GMT_FFT_K_IS_KR;
491 }
492 else if (line[i] == 'x') {
493 i++; f_info->k_type = GMT_FFT_K_IS_KX;
494 }
495 else if (line[i] == 'y') {
496 i++; f_info->k_type = GMT_FFT_K_IS_KY;
497 }
498 fourvals[0] = fourvals[1] = fourvals[2] = fourvals[3] = GMT_NOTSET;
499
500 n_tokens = pos = 0;
501 while ((gmt_strtok (&line[i], "/", &pos, p))) {
502 if (n_tokens > 3) {
503 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Too many slashes in -F.\n");
504 return (true);
505 }
506 if(p[0] == '-')
507 fourvals[n_tokens] = GMT_NOTSET;
508 else {
509 if ((sscanf(p, "%lf", &fourvals[n_tokens])) != 1) {
510 GMT_Report (GMT->parent, GMT_MSG_ERROR, " Cannot read token %d.\n", n_tokens);
511 return (true);
512 }
513 }
514 n_tokens++;
515 }
516
517 if (!(n_tokens == 2 || n_tokens == 3 || n_tokens == 4)) {
518 GMT_Report (GMT->parent, GMT_MSG_ERROR, "-F Cannot find 2-4 tokens separated by slashes.\n");
519 return (true);
520 }
521 descending = true;
522 if (f_info->kind == GRDFFT_FILTER_BW && n_tokens == 3) n_tokens = 2; /* So we don't check the order as a wavelength */
523
524 for (i = 1; i < n_tokens; i++) {
525 if (fourvals[i] == GMT_NOTSET || fourvals[i-1] == GMT_NOTSET) continue;
526 if (fourvals[i] > fourvals[i-1]) descending = false;
527 }
528 if (!(descending)) {
529 GMT_Report (GMT->parent, GMT_MSG_ERROR, "-F Wavelengths are not in descending order.\n");
530 return (true);
531 }
532 j = f_info->k_type;
533 if (f_info->kind == GRDFFT_FILTER_COS) { /* Cosine band-pass specification */
534 if ((fourvals[0] * fourvals[1]) < 0.0 || (fourvals[2] * fourvals[3]) < 0.0) {
535 GMT_Report (GMT->parent, GMT_MSG_ERROR, "-F Pass/Cut specification error.\n");
536 return (true);
537 }
538
539 /* Now everything is OK */
540
541 if (fourvals[0] >= 0.0 || fourvals[1] >= 0.0) { /* Lower end values are set */
542 f_info->lc[j] = (2.0 * M_PI)/fourvals[0];
543 f_info->lp[j] = (2.0 * M_PI)/fourvals[1];
544 if (fourvals[0] != fourvals[1]) f_info->ltaper[j] = 1.0/(f_info->lc[j] - f_info->lp[j]);
545 }
546
547 if (fourvals[2] >= 0.0 || fourvals[3] >= 0.0) { /* Higher end values are set */
548 f_info->hp[j] = (2.0 * M_PI)/fourvals[2];
549 f_info->hc[j] = (2.0 * M_PI)/fourvals[3];
550 if (fourvals[2] != fourvals[3]) f_info->htaper[j] = 1.0/(f_info->hc[j] - f_info->hp[j]);
551 }
552 f_info->filter = &grdfft_cosine_weight_grdfft;
553 }
554 else if (f_info->kind == GRDFFT_FILTER_BW) { /* Butterworth specification */
555 f_info->llambda[j] = (fourvals[0] == GMT_NOTSET) ? -1.0 : fourvals[0] / TWO_PI; /* TWO_PI is used to counteract the 2*pi in the wavenumber */
556 f_info->hlambda[j] = (fourvals[1] == GMT_NOTSET) ? -1.0 : fourvals[1] / TWO_PI;
557 f_info->bw_order = 2.0 * fourvals[2];
558 f_info->filter = &grdfft_bw_weight;
559 }
560 else { /* Gaussian half-amp specifications */
561 f_info->llambda[j] = (fourvals[0] == GMT_NOTSET) ? -1.0 : fourvals[0] / TWO_PI; /* TWO_PI is used to counteract the 2*pi in the wavenumber */
562 f_info->hlambda[j] = (fourvals[1] == GMT_NOTSET) ? -1.0 : fourvals[1] / TWO_PI;
563 f_info->filter = &grdfft_gauss_weight;
564 }
565 f_info->arg = f_info->kind - GRDFFT_FILTER_EXP;
566 return (false);
567 }
568
usage(struct GMTAPI_CTRL * API,int level)569 static int usage (struct GMTAPI_CTRL *API, int level) {
570 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
571 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
572 GMT_Usage (API, 0, "usage: %s %s [<ingrid2>] [-A<azimuth>] [-C<zlevel>] [-D[<scale>|g]] "
573 "[-E[r|x|y][+w[k]][+n]] [-F[r|x|y]<parameters>] [-G<outgrid>|<table>] [-I[<scale>|g]] [-N%s] [-Q] "
574 "[-S<scale>|d] [%s] [-fg] [%s] [%s]\n", name, GMT_INGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_ho_OPT, GMT_PAR_OPT);
575
576 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
577
578 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
579 gmt_ingrid_syntax (API, 0, "Name of input grid. For cross-spectrum also supply <ingrid2> (same modifiers apply)");
580 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
581 GMT_Usage (API, 1, "\n-A<azimuth>");
582 GMT_Usage (API, -2, "Take azimuthal derivative along line <azimuth> degrees CW from North.");
583 GMT_Usage (API, 1, "\n-C<zlevel>");
584 GMT_Usage (API, -2, "Continue field upward (+) or downward (-) to <zlevel> (meters). ");
585 GMT_Usage (API, 1, "\n-D[<scale>|g]");
586 GMT_Usage (API, -2, "Differentiate, i.e., multiply by kr [ and optionally by <scale>]. Use -Dg to get mGal from m].");
587 GMT_Usage (API, 1, "\n-E[r|x|y][+w[k]][+n]");
588 GMT_Usage (API, -2, "Estimate spEctrum in the radial r [Default], x, or y-direction. "
589 "Given one grid X, write f, Xpower[f] to output file (see -G) or standard output. "
590 "Given two grids X and Y, write f, Xpower[f], Ypower[f], coherent power[f], "
591 "noise power[f], phase[f], admittance[f], gain[f], coherency[f]. "
592 "Each quantity is followed by a column of 1 std dev. error estimates.");
593 GMT_Usage (API, 3, "+w Write wavelength instead of frequency; append k to report "
594 "wavelength in km (geographic grids only) [m].");
595 GMT_Usage (API, 3, "+n Yield mean power instead of total power per frequency.");
596 GMT_Usage (API, 1, "\n-F[r|x|y]<parameters>");
597 GMT_Usage (API, -2, "Filter r [x] [y] frequencies according to one of three kinds of filter specifications:");
598 GMT_Usage (API, 3, "%s Cosine band-pass: Append four wavelengths <lc>/<lp>/<hp>/<hc>. "
599 "frequencies outside <lc>/<hc> are cut; inside <lp>/<hp> are passed, rest are tapered. "
600 "Replace wavelength by - to skip, e.g., -F-/-/500/100 is a low-pass filter.\n", GMT_LINE_BULLET);
601 GMT_Usage (API, 3, "%s Gaussian band-pass: Append two wavelengths <lo>/<hi> where filter amplitudes = 0.5. "
602 "Replace wavelength by - to skip, e.g., -F300/- is a high-pass Gaussian filter.", GMT_LINE_BULLET);
603 GMT_Usage (API, 3, "%s Butterworth band-pass: Append two cut-off wavelengths and order, i.e., <lo>/<hi>/<order> "
604 "where filter amplitudes are 0.707. Replace wavelength by - to skip, e.g., "
605 "try -F300/-/2 for a high-pass 2nd-order Butterworth filter.", GMT_LINE_BULLET);
606 GMT_Usage (API, 1, "\n-G<outgrid>|<table>");
607 GMT_Usage (API, -2, "Filename for output grid file OR 1-D data table (see -E). "
608 "Note: Optional for -E (spectrum written to standard output); required otherwise.");
609 GMT_Usage (API, 1, "\n-I[<scale>|g]");
610 GMT_Usage (API, -2, "Integrate, i.e., divide by kr [ and optionally by scale]. Use -Ig to get m from mGal].");
611 GMT_FFT_Option (API, 'N', GMT_FFT_DIM, "Choose or inquire about suitable grid dimensions for FFT, and set modifiers.");
612 GMT_Usage (API, 1, "\n-Q Perform no operations, just do forward FFF and write output if selected in -N.");
613 GMT_Usage (API, 1, "\n-S<scale>|d");
614 GMT_Usage (API, -2, "Multiply field by <scale> after inverse FFT [1.0]. "
615 "Alternatively, append d to convert deflection of vertical to micro-radians.");
616 GMT_Option (API, "V");
617 GMT_Usage (API, 1, "\n-fg Convert geographic grids to meters using a \"Flat Earth\" approximation.");
618 GMT_Usage (API, 1, "\n-ho Write header record for spectral estimates (requires -E) [no header].");
619 GMT_Option (API, ".");
620 GMT_Usage (API, 1, "\nNote: List operations in the order desired for execution.");
621
622 return (GMT_MODULE_USAGE);
623 }
624
grdfft_add_operation(struct GMT_CTRL * GMT,struct GRDFFT_CTRL * Ctrl,int operation,unsigned int n_par,double * par)625 GMT_LOCAL void grdfft_add_operation (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, int operation, unsigned int n_par, double *par) {
626 Ctrl->n_op_count++;
627 Ctrl->operation = gmt_M_memory (GMT, Ctrl->operation, Ctrl->n_op_count, int);
628 Ctrl->operation[Ctrl->n_op_count-1] = operation;
629 if (n_par) {
630 Ctrl->par = gmt_M_memory (GMT, Ctrl->par, Ctrl->n_par + n_par, double);
631 gmt_M_memcpy (&Ctrl->par[Ctrl->n_par], par, n_par, double);
632 Ctrl->n_par += n_par;
633 }
634 }
635
parse(struct GMT_CTRL * GMT,struct GRDFFT_CTRL * Ctrl,struct F_INFO * f_info,struct GMT_OPTION * options)636 static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO *f_info, struct GMT_OPTION *options) {
637 /* This parses the options provided to grdfft and sets parameters in Ctrl.
638 * Note Ctrl has already been initialized and non-zero default values set.
639 * Any GMT common options will override values set previously by other commands.
640 * It also replaces any file names specified as input or output with the data ID
641 * returned when registering these sources/destinations with the API.
642 */
643
644 unsigned int j, k, pos, n_errors = 0, filter_type = 0;
645 int n_scan;
646 double par[5];
647 char combined[GMT_BUFSIZ] = {""}, argument[GMT_LEN16] = {""}, p[GMT_LEN64] = {""}, *c = NULL;
648 struct GMT_OPTION *opt = NULL, *ptr = NULL;
649 struct GMTAPI_CTRL *API = GMT->parent;
650
651 if (gmt_M_compat_check (GMT, 4)) {
652 char *mod = NULL;
653 if ((ptr = GMT_Find_Option (API, 'L', options))) { /* Gave old -L */
654 mod = ptr->arg; /* Gave old -L option */
655 if (mod[0] == '\0') strcat (argument, "+l"); /* Leave trend alone -L */
656 else if (mod[0] == 'm') strcat (argument, "+a"); /* Remove mean -Lm */
657 else if (mod[0] == 'h') strcat (argument, "+h"); /* Remove mid-value -Lh */
658 }
659 }
660
661 gmt_M_memset (f_info, 1, struct F_INFO);
662 for (j = 0; j < 3; j++) {
663 f_info->lc[j] = f_info->lp[j] = -1.0; /* Set negative, below valid frequency range */
664 f_info->hp[j] = f_info->hc[j] = DBL_MAX; /* Set huge positive, above valid frequency range */
665 }
666
667 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
668 gmt_M_memset (par, 5, double);
669
670 switch (opt->option) {
671 case '<': /* Input file (only 1 or 2 are accepted) */
672 Ctrl->In.active = true;
673 if (Ctrl->In.n_grids >= 2) {
674 n_errors++;
675 GMT_Report (API, GMT_MSG_ERROR, "A maximum of two input grids may be processed\n");
676 }
677 else {
678 Ctrl->In.file[Ctrl->In.n_grids] = strdup (opt->arg);
679 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file[Ctrl->In.n_grids]))) n_errors++;;
680 Ctrl->In.n_grids++;
681 }
682 break;
683
684 /* Processes program-specific parameters */
685
686 case 'A': /* Directional derivative */
687 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
688 Ctrl->A.active = true;
689 n_errors += gmt_M_check_condition (GMT, sscanf(opt->arg, "%lf", &par[0]) != 1,
690 "Option -A: Cannot read azimuth\n");
691 grdfft_add_operation (GMT, Ctrl, GRDFFT_AZIMUTHAL_DERIVATIVE, 1, par);
692 break;
693 case 'C': /* Upward/downward continuation */
694 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
695 Ctrl->C.active = true;
696 n_errors += gmt_M_check_condition (GMT, sscanf(opt->arg, "%lf", &par[0]) != 1,
697 "Option -C: Cannot read zlevel\n");
698 grdfft_add_operation (GMT, Ctrl, GRDFFT_UP_DOWN_CONTINUE, 1, par);
699 break;
700 case 'D': /* d/dz */
701 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
702 Ctrl->D.active = true;
703 par[0] = (opt->arg[0]) ? ((opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg)) : 1.0;
704 n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -D: scale must be nonzero\n");
705 grdfft_add_operation (GMT, Ctrl, GRDFFT_DIFFERENTIATE, 1, par);
706 break;
707 case 'E': /* x,y,or radial spectrum, w for wavelength; k for km if geographical, n for normalize */
708 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
709 /* Old syntax: -E[x|y|r][w[k]][n] new syntax: -E[x|y|r][+w[k]][+n] */
710 Ctrl->E.active = true;
711 j = 0;
712 switch (opt->arg[0]) { /* Determine direction for spectrum or default to radial */
713 case 'r': Ctrl->E.mode = 0; j = 1; break;
714 case 'x': Ctrl->E.mode = +1; j = 1; break;
715 case 'y': Ctrl->E.mode = -1; j = 1; break;
716 default :Ctrl->E.mode = 0; break;
717 }
718 if ((c = gmt_first_modifier (GMT, opt->arg, "wn"))) { /* Process new modifiers [+w[k]][+n] */
719 pos = 0; /* Reset to start of new word */
720 while (gmt_getmodopt (GMT, 'E', c, "wn", &pos, p, &n_errors) && n_errors == 0) {
721 switch (p[0]) {
722 case 'w': Ctrl->E.give_wavelength = true; if (p[1] == 'k') Ctrl->E.km = true; break;
723 case 'n': Ctrl->E.normalize = true; break;
724 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
725 }
726 }
727 }
728 else { /* Old-style modifiers [w[k]][n] */
729 while (opt->arg[j]) {
730 switch (opt->arg[j]) {
731 case 'w': Ctrl->E.give_wavelength = true; break;
732 case 'k': if (Ctrl->E.give_wavelength) Ctrl->E.km = true; break;
733 case 'n': Ctrl->E.normalize = true; break;
734 }
735 j++;
736 }
737 }
738 par[0] = Ctrl->E.mode;
739 grdfft_add_operation (GMT, Ctrl, GRDFFT_SPECTRUM, 1, par);
740 break;
741 case 'F': /* Filter */
742 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
743 Ctrl->F.active = true;
744 if (!(f_info->set_already)) {
745 filter_type = gmt_count_char (GMT, opt->arg, '/');
746 f_info->kind = GRDFFT_FILTER_EXP + (filter_type - 1);
747 f_info->set_already = true;
748 grdfft_add_operation (GMT, Ctrl, f_info->kind, 0, NULL);
749 }
750 n_errors += gmt_M_check_condition (GMT, grdfft_parse_f_string (GMT, f_info, opt->arg), "Option -F");
751 break;
752 case 'G': /* Output file */
753 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
754 Ctrl->G.active = true;
755 Ctrl->G.file = strdup (opt->arg);
756 break;
757 case 'I': /* Integrate */
758 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
759 Ctrl->I.active = true;
760 par[0] = (opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg);
761 n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -I: scale must be nonzero\n");
762 grdfft_add_operation (GMT, Ctrl, GRDFFT_INTEGRATE, 1, par);
763 break;
764 case 'L': /* Leave trend alone */
765 if (gmt_M_compat_check (GMT, 4))
766 GMT_Report (API, GMT_MSG_COMPAT, "Option -L is deprecated; use -N modifiers in the future.\n");
767 else
768 n_errors += gmt_default_error (GMT, opt->option);
769 break;
770 case 'M': /* Geographic data */
771 if (gmt_M_compat_check (GMT, 4)) {
772 GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
773 if (gmt_M_is_cartesian (GMT, GMT_IN))
774 gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
775 }
776 else
777 n_errors += gmt_default_error (GMT, opt->option);
778 break;
779 case 'N': /* Grid dimension setting or inquiery */
780 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
781 Ctrl->N.active = true;
782 if (ptr && gmt_M_compat_check (GMT, 4)) { /* Got both old -L and -N; append */
783 sprintf (combined, "%s%s", opt->arg, argument);
784 Ctrl->N.info = GMT_FFT_Parse (API, 'N', GMT_FFT_DIM, combined);
785 }
786 else
787 Ctrl->N.info = GMT_FFT_Parse (API, 'N', GMT_FFT_DIM, opt->arg);
788 if (Ctrl->N.info == NULL) n_errors++;
789 break;
790 case 'Q': /* No operator */
791 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
792 Ctrl->Q.active = true;
793 break;
794 case 'S': /* Scale */
795 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
796 Ctrl->S.active = true;
797 Ctrl->S.scale = (opt->arg[0] == 'd' || opt->arg[0] == 'D') ? 1.0e6 : atof (opt->arg);
798 break;
799 case 'T': /* Flexural isostasy */
800 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
801 if (gmt_M_compat_check (GMT, 4)) {
802 GMT_Report (API, GMT_MSG_COMPAT,
803 "Option -T is deprecated; see gravfft for isostasy and gravity calculations.\n");
804 Ctrl->T.active = true;
805 n_scan = sscanf (opt->arg, "%lf/%lf/%lf/%lf/%lf", &par[0], &par[1], &par[2], &par[3], &par[4]);
806 for (j = 1, k = 0; j < 5; j++) if (par[j] < 0.0) k++;
807 n_errors += gmt_M_check_condition (GMT, n_scan != 5 || k > 0,
808 "Option -T: Correct syntax:\n\t-T<te>/<rhol>/<rhom>/<rhow>/<rhoi>, all densities >= 0\n");
809 grdfft_add_operation (GMT, Ctrl, GRDFFT_ISOSTASY, 5, par);
810 }
811 else
812 n_errors += gmt_default_error (GMT, opt->option);
813 break;
814 #ifdef DEBUG
815 case '=': /* Do nothing */
816 grdfft_add_operation (GMT, Ctrl, GRDFFT_NOTHING, 1, par);
817 if (opt->arg[0] == '+') show_n = true;
818 break;
819 #endif
820 default: /* Report bad options */
821 n_errors += gmt_default_error (GMT, opt->option);
822 break;
823 }
824 }
825 if (gmt_M_compat_check (GMT, 4) && !Ctrl->N.active && ptr) { /* User set -L but no -N so nothing got appended above... Sigh...*/
826 Ctrl->N.info = GMT_FFT_Parse (API, 'N', GMT_FFT_DIM, argument);
827 }
828 if (Ctrl->N.info && Ctrl->N.active && Ctrl->N.info->info_mode == GMT_FFT_LIST)
829 return (GMT_PARSE_ERROR); /* So that we exit the program */
830
831 if (Ctrl->N.info && gmt_M_compat_check (GMT, 4) && Ctrl->T.active)
832 Ctrl->N.info->trend_mode = GMT_FFT_REMOVE_NOTHING;
833
834 n_errors += gmt_M_check_condition (GMT, !(Ctrl->n_op_count) && !Ctrl->Q.active, "Must specify at least one operation (unless -Q is set)\n");
835 n_errors += gmt_M_check_condition (GMT, Ctrl->n_op_count && Ctrl->Q.active, "Cannot specify operations if -Q is set\n");
836 n_errors += gmt_M_check_condition (GMT, Ctrl->S.scale == 0.0, "Option -S: scale must be nonzero\n");
837 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file[0], "Must specify input file\n");
838 n_errors += gmt_M_check_condition (GMT, !Ctrl->E.active && !Ctrl->Q.active && !Ctrl->G.file, "Option -G: Must specify output grid file unless -E or -Q is set\n");
839
840 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
841 }
842
843 #define bailout(code) {gmt_M_free_options (mode); return (code);}
844 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
845
GMT_grdfft(void * V_API,int mode,void * args)846 EXTERN_MSC int GMT_grdfft (void *V_API, int mode, void *args) {
847 int error = 0, status;
848 unsigned int op_count = 0, par_count = 0, k;
849 char *spec_msg[2] = {"spectrum", "cross-spectrum"};
850 struct GMT_GRID *Grid[2] = {NULL, NULL}, *Orig[2] = {NULL, NULL};
851 struct F_INFO f_info;
852 struct GMT_FFT_WAVENUMBER *FFT_info[2] = {NULL, NULL}, *K = NULL;
853 struct GRDFFT_CTRL *Ctrl = NULL;
854 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
855 struct GMT_OPTION *options = NULL;
856 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
857
858 /*----------------------- Standard module initialization and parsing ----------------------*/
859
860 if (API == NULL) return (GMT_NOT_A_SESSION);
861 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
862 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
863
864 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
865
866 /* Parse the command-line arguments */
867
868 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
869 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
870 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
871 if ((error = parse (GMT, Ctrl, &f_info, options)) != 0) Return (error);
872
873 /*---------------------------- This is the grdfft main code ----------------------------*/
874
875 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid(s)\n");
876 for (k = 0; k < Ctrl->In.n_grids; k++) { /* First read the grid header(s) */
877 if ((Orig[k] = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file[k], NULL)) == NULL)
878 Return (API->error);
879 }
880
881 if (Ctrl->In.n_grids == 2) { /* If given 2 grids, make sure they are co-registered and has same size, registration, etc. */
882 if(Orig[0]->header->registration != Orig[1]->header->registration) {
883 GMT_Report (API, GMT_MSG_ERROR, "The two grids have different registrations!\n");
884 Return (GMT_RUNTIME_ERROR);
885 }
886 if (!gmt_M_grd_same_shape (GMT, Orig[0], Orig[1])) {
887 GMT_Report (API, GMT_MSG_ERROR, "The two grids have different dimensions\n");
888 Return (GMT_RUNTIME_ERROR);
889 }
890 if (!gmt_M_grd_same_region (GMT, Orig[0], Orig[1])) {
891 GMT_Report (API, GMT_MSG_ERROR, "The two grids have different regions\n");
892 Return (GMT_RUNTIME_ERROR);
893 }
894 if (!gmt_M_grd_same_inc (GMT, Orig[0], Orig[1])) {
895 GMT_Report (API, GMT_MSG_ERROR, "The two grids have different intervals\n");
896 Return (GMT_RUNTIME_ERROR);
897 }
898 }
899
900 /* Grids are compatible. Initialize FFT structs, grid headers, read data, and check for NaNs */
901
902 for (k = 0; k < Ctrl->In.n_grids; k++) { /* Read, and check that no NaNs are present in either grid */
903 /* Note: If input grid(s) are read-only then we must duplicate them; otherwise Grid[k] points to Orig[k]
904 * From here we address the first grid via Grid[0] and the 2nd grid (if given) as Grid[1];
905 * we are done with using the addresses Orig[k] directly. */
906 (void) gmt_set_outgrid (GMT, Ctrl->In.file[k], false, 0, Orig[k], &Grid[k]);
907 FFT_info[k] = GMT_FFT_Create (API, Grid[k], GMT_FFT_DIM, GMT_GRID_IS_COMPLEX_REAL, Ctrl->N.info);
908 }
909 K = FFT_info[0]; /* We only need one of these anyway; K is a shorthand */
910
911 #ifdef FTEST
912 /* PW: Used with -DFTEST to check that the radial filters compute correctly */
913 {
914 double f = 0.0;
915 while (f < 3.0) {
916 printf ("%g\t%g\n", f, f_info.filter (f, 0)); /* Radial filter */
917 f += 0.01;
918 }
919 Return (GMT_RUNTIME_ERROR);
920 }
921 #endif
922
923 for (k = 0; k < Ctrl->In.n_grids; k++) { /* Call the forward FFT, once per grid, optionally save raw FFT output */
924 GMT_Report (API, GMT_MSG_INFORMATION, "forward FFT...\n");
925 if (GMT_FFT (API, Grid[k], GMT_FFT_FWD, GMT_FFT_COMPLEX, FFT_info[k]))
926 Return (GMT_RUNTIME_ERROR);
927 }
928
929 if (Ctrl->Q.active) GMT_Report (API, GMT_MSG_INFORMATION, "No wavenumber operations selected\n");
930
931 for (op_count = par_count = 0; op_count < Ctrl->n_op_count; op_count++) {
932 switch (Ctrl->operation[op_count]) {
933 case GRDFFT_UP_DOWN_CONTINUE:
934 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING))
935 ((Ctrl->par[par_count] < 0.0) ? GMT_Report (API, GMT_MSG_INFORMATION, "downward continuation...\n") :
936 GMT_Report (API, GMT_MSG_INFORMATION, "upward continuation...\n"));
937 par_count += grdfft_do_continuation (Grid[0], &Ctrl->par[par_count], K);
938 break;
939 case GRDFFT_AZIMUTHAL_DERIVATIVE:
940 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "azimuthal derivative...\n");
941 par_count += grdfft_do_azimuthal_derivative (Grid[0], &Ctrl->par[par_count], K);
942 break;
943 case GRDFFT_DIFFERENTIATE:
944 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "differentiate...\n");
945 par_count += grdfft_do_differentiate (Grid[0], &Ctrl->par[par_count], K);
946 break;
947 case GRDFFT_INTEGRATE:
948 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "integrate...\n");
949 par_count += grdfft_do_integrate (Grid[0], &Ctrl->par[par_count], K);
950 break;
951 case GRDFFT_ISOSTASY:
952 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "isostasy...\n");
953 par_count += grdfft_do_isostasy (Grid[0], Ctrl, &Ctrl->par[par_count], K);
954 break;
955 case GRDFFT_FILTER_COS:
956 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "cosine filter...\n");
957 grdfft_do_filter (Grid[0], &f_info, K);
958 break;
959 case GRDFFT_FILTER_EXP:
960 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "Gaussian filter...\n");
961 grdfft_do_filter (Grid[0], &f_info, K);
962 break;
963 case GRDFFT_FILTER_BW:
964 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "Butterworth filter...\n");
965 grdfft_do_filter (Grid[0], &f_info, K);
966 break;
967 case GRDFFT_SPECTRUM: /* This operator writes a table to file (or stdout if -G is not used) */
968 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "%s...\n", spec_msg[Ctrl->In.n_grids-1]);
969 status = grdfft_do_spectrum (GMT, Grid[0], Grid[1], &Ctrl->par[par_count], Ctrl->E.give_wavelength, Ctrl->E.km, Ctrl->E.normalize, Ctrl->G.file, K);
970 if (status < 0) Return (status);
971 par_count += status;
972 break;
973 case GRDFFT_NOTHING:
974 break;
975 }
976 }
977
978 if (!Ctrl->E.active && !Ctrl->Q.active) { /* Since -E output is handled separately by grdfft_do_spectrum itself */
979 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) GMT_Report (API, GMT_MSG_INFORMATION, "inverse FFT...\n");
980
981 if (GMT_FFT (API, Grid[0], GMT_FFT_INV, GMT_FFT_COMPLEX, K))
982 Return (GMT_RUNTIME_ERROR);
983 #ifdef DEBUG
984 gmt_grd_dump (Grid[0]->header, Grid[0]->data, false, "After Inv FFT");
985 #endif
986
987 if (!doubleAlmostEqual (Ctrl->S.scale, 1.0)) gmt_scale_and_offset_f (GMT, Grid[0]->data, Grid[0]->header->size, Ctrl->S.scale, 0);
988
989 /* The data are in the middle of the padded array; only the interior (original dimensions) will be written to file */
990 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid[0])) Return (API->error);
991 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY |
992 GMT_GRID_IS_COMPLEX_REAL, NULL, Ctrl->G.file, Grid[0]) != GMT_NOERROR) {
993 Return (API->error);
994 }
995 }
996
997 for (k = 0; k < Ctrl->In.n_grids; k++)
998 GMT_FFT_Destroy (API, &(FFT_info[k]));
999
1000 Return (GMT_NOERROR);
1001 }
1002