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 * API functions to support the grdredpol application.
19 *
20 * Brief synopsis: grdredpol reads a grid file with the magnetic field anomaly and
21 * computes the continuous reduction to the pole (RTP) anomaly by calculating the
22 * filter coefficients in the frequency, inverse FT and convolve in space domain.
23 * For details on the method see, Luis, J.F and J.M. Miranda (2008),
24 * "Reevaluation of magnetic chrons in the North Atlantic between 35N and 47N:
25 * Implications for the formation of the Azores Triple Junction and associated plateau,
26 * J. Geophys. Res., 113, B10105, doi:10.1029/2007JB005573
27 *
28 * Author: Joaquim Luis / Miguel Miranda
29 * Date: 18-Feb-2003 (original GMT4 version)
30 * Version: 6 API
31 */
32
33 #include "gmt_dev.h"
34
35 #define THIS_MODULE_CLASSIC_NAME "grdredpol"
36 #define THIS_MODULE_MODERN_NAME "grdredpol"
37 #define THIS_MODULE_LIB "potential"
38 #define THIS_MODULE_PURPOSE "Compute the Continuous Reduction To the Pole, AKA differential RTP"
39 #define THIS_MODULE_KEYS "<G{,EG(,GG},ZG)"
40 #define THIS_MODULE_NEEDS "g"
41 #define THIS_MODULE_OPTIONS "-RVn"
42
43 struct GRDREDPOL_CTRL {
44 struct GRDREDPOL_In {
45 bool active;
46 char *file;
47 } In;
48 struct GRDREDPOL_C { /* -C */
49 bool active;
50 bool use_igrf;
51 bool const_f;
52 double dec;
53 double dip;
54 } C;
55 struct GRDREDPOL_E { /* -E */
56 bool active;
57 bool dip_grd_only;
58 bool dip_dec_grd;
59 char *decfile;
60 char *dipfile;
61 } E;
62 struct GRDREDPOL_F { /* -F */
63 bool active;
64 unsigned int ncoef_row;
65 unsigned int ncoef_col;
66 unsigned int compute_n; /* Compute ncoef_col */
67 double width;
68 } F;
69 struct GRDREDPOL_G { /* -G<file> */
70 bool active;
71 char *file;
72 } G;
73 struct GRDREDPOL_M { /* -M */
74 bool active;
75 bool pad_zero;
76 bool mirror;
77 } M;
78 struct GRDREDPOL_N { /* -N */
79 bool active;
80 } N;
81 struct GRDREDPOL_S { /* -S, size of working grid */
82 bool active;
83 unsigned int n_columns;
84 unsigned int n_rows;
85 } S;
86 struct GRDREDPOL_T { /* -T */
87 bool active;
88 double year;
89 } T;
90 struct GRDREDPOL_W { /* -W */
91 double wid;
92 bool active;
93 } W;
94 struct GRDREDPOL_Z { /* -Z */
95 bool active;
96 char *file;
97 } Z;
98 };
99
100
101 #define ij0_data(Ctrl,i,j) ((Ctrl->S.n_columns+Ctrl->F.ncoef_col-1)*(i)+(j))
102 #define ij_mn(Ctrl,i,j) (Ctrl->F.ncoef_row*(j)+(i))
103
New_Ctrl(struct GMT_CTRL * GMT)104 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
105 struct GRDREDPOL_CTRL *C;
106
107 C = gmt_M_memory (GMT, NULL, 1, struct GRDREDPOL_CTRL);
108
109 /* Initialize values whose defaults are not 0/false/NULL */
110 C->C.use_igrf = true;
111 C->M.pad_zero = true;
112 C->N.active = true;
113 C->F.ncoef_row = 25;
114 C->F.ncoef_col = 25;
115 C->T.year = 2000;
116 C->W.wid = 5;
117 return (C);
118 }
119
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDREDPOL_CTRL * C)120 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDREDPOL_CTRL *C) { /* Deallocate control structure */
121 if (!C) return;
122 gmt_M_str_free (C->In.file);
123 gmt_M_str_free (C->G.file);
124 gmt_M_str_free (C->E.dipfile);
125 gmt_M_str_free (C->E.decfile);
126 gmt_M_str_free (C->Z.file);
127 gmt_M_free (GMT, C);
128 }
129
grdgravmag3d_rtp_filt_colinear(int i,int j,int n21,double * gxr,double * gxi,double * gxar,double * gxai,double * gxbr,double * gxbi,double * gxgr,double * gxgi,double u,double v,double alfa,double beta,double gama,struct GRDREDPOL_CTRL * Ctrl)130 GMT_LOCAL void grdgravmag3d_rtp_filt_colinear (int i, int j, int n21, double *gxr,double *gxi, double *gxar,
131 double *gxai, double *gxbr, double *gxbi, double *gxgr, double *gxgi, double u,
132 double v, double alfa, double beta, double gama, struct GRDREDPOL_CTRL *Ctrl) {
133
134 uint64_t ij = ij_mn(Ctrl,i,j-n21+1);
135 double ro, ro2, ro3, ro4, ro5, alfa_u, beta_v, gama_ro, gama_ro_2;
136 double alfa_u_beta_v, alfa_u_beta_v_2, rnr, rni, t2, t3;
137 ro2 = u * u + v * v; ro = sqrt(ro2); ro3 = ro2 * ro;
138 ro4 = ro2 * ro2; ro5 = ro3 * ro2;
139 alfa_u = alfa * u;
140 beta_v = beta * v;
141 gama_ro = gama * ro; gama_ro_2 = gama_ro * gama_ro;
142 alfa_u_beta_v = alfa_u + beta_v;
143 alfa_u_beta_v_2 = alfa_u_beta_v * alfa_u_beta_v;
144
145 t2 = 1 / ((alfa_u_beta_v_2 + gama_ro_2) * (alfa_u_beta_v_2 + gama_ro_2));
146 t3 = t2 / (alfa_u_beta_v_2 + gama_ro_2);
147
148 rnr = (gama_ro_2 - alfa_u_beta_v_2) * ro2;
149 rni = 2 * gama_ro * alfa_u_beta_v * ro2;
150
151 gxr[ij] = rnr * t2;
152 gxi[ij] = rni * t2;
153
154 if (Ctrl->N.active) { /* Use the Taylor expansion */
155 gxar[ij] = -2*alfa_u_beta_v*u*ro2 * t2 - 4*(gama_ro_2-alfa_u_beta_v_2)*ro2*alfa_u_beta_v*u * t3;
156 gxai[ij] = 2*gama*ro3*u * t2 - 8*gama*ro3*alfa_u_beta_v_2*u * t3;
157 gxbr[ij] = -2*alfa_u_beta_v*v*ro2 * t2 - 4*(gama_ro_2-alfa_u_beta_v_2)*ro2*alfa_u_beta_v*v * t3;
158 gxbi[ij] = 2*gama*ro3*v * t2 - 8*gama*ro3*alfa_u_beta_v_2*v * t3;
159 gxgr[ij] = 2*gama*ro4 * t2 - 4*(gama_ro_2-alfa_u_beta_v_2)*ro4*gama * t3;
160 gxgi[ij] = 2*ro3*alfa_u_beta_v * t2 - 8*gama*gama*ro5*alfa_u_beta_v * t3;
161 }
162 }
163
164
grdgravmag3d_rtp_filt_not_colinear(int i,int j,int n21,double * gxr,double * gxi,double * gxar,double * gxai,double * gxbr,double * gxbi,double * gxgr,double * gxgi,double * gxtr,double * gxti,double * gxmr,double * gxmi,double * gxnr,double * gxni,double u,double v,double alfa,double beta,double gama,double tau,double mu,double nu,struct GRDREDPOL_CTRL * Ctrl)165 GMT_LOCAL void grdgravmag3d_rtp_filt_not_colinear (int i, int j, int n21, double *gxr, double *gxi, double *gxar,
166 double *gxai, double *gxbr, double *gxbi, double *gxgr, double *gxgi, double *gxtr,
167 double *gxti, double *gxmr, double *gxmi, double *gxnr, double *gxni, double u, double v, double alfa,
168 double beta, double gama, double tau, double mu, double nu, struct GRDREDPOL_CTRL *Ctrl) {
169
170 uint64_t ij = ij_mn(Ctrl,i,j-n21+1);
171 double ro, ro2, ro3, ro4, ro5, alfa_u, beta_v, gama_ro, gama_ro_2;
172 double alfa_u_beta_v, alfa_u_beta_v_2, rnr, rni, den_r, den_i1, den_i2;
173 double tau_u, mu_v, nu_ro, nu_ro_2, tau_u_mu_v, tau_u_mu_v_2, tau_u_mu_v_gama;
174 double alfa_u_beta_v_nu, alfa_u_beta_v_tau_u_mu_v;
175
176 ro2 = u * u + v * v; ro = sqrt(ro2); ro3 = ro2 * ro;
177 ro4 = ro2 * ro2; ro5 = ro3 * ro2;
178 alfa_u = alfa * u;
179 beta_v = beta * v;
180 gama_ro = gama * ro; gama_ro_2 = gama_ro * gama_ro;
181 alfa_u_beta_v = alfa_u + beta_v;
182 alfa_u_beta_v_2 = alfa_u_beta_v * alfa_u_beta_v;
183
184 tau_u = tau * u; mu_v = mu * v;
185 nu_ro = nu * ro; nu_ro_2 = nu_ro * nu_ro;
186 tau_u_mu_v = tau_u + mu_v;
187 tau_u_mu_v_2 = tau_u_mu_v * tau_u_mu_v;
188
189 tau_u_mu_v_gama = tau_u_mu_v * gama;
190 alfa_u_beta_v_nu = alfa_u_beta_v * nu;
191 alfa_u_beta_v_tau_u_mu_v = alfa_u_beta_v * tau_u_mu_v;
192
193 rnr = (-alfa_u_beta_v_tau_u_mu_v + gama_ro * nu_ro) * ro2;
194 rni = (alfa_u_beta_v_nu + tau_u_mu_v_gama) * ro3;
195
196 den_r = (alfa_u_beta_v_2 + gama_ro_2)*(tau_u_mu_v_2 + nu_ro_2);
197 den_i1 = den_r * (alfa_u_beta_v_2+gama_ro_2);
198 den_i2 = den_r * (tau_u_mu_v_2 + nu_ro_2);
199
200 gxr[ij] = rnr / den_r;
201 gxi[ij] = rni / den_r;
202
203 if (Ctrl->N.active) { /* Use the Taylor expansion */
204 gxar[ij] = -u*tau_u_mu_v*ro2/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro2*alfa_u_beta_v*u/den_i1;
205 gxai[ij] = u*nu*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro3*alfa_u_beta_v*u/den_i1;
206 gxbr[ij] = -v*tau_u_mu_v*ro2/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro2*alfa_u_beta_v*v/den_i1;
207 gxbi[ij] = v*nu*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro3*alfa_u_beta_v*v/den_i1;
208 gxgr[ij] = nu*ro4/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro4*gama/den_i1;
209 gxgi[ij] = tau_u_mu_v*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro5*gama/den_i1;
210 gxtr[ij] = -alfa_u_beta_v*u*ro2/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro2*u*tau_u_mu_v/den_i2;
211 gxti[ij] = u*gama*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro3*u*tau_u_mu_v/den_i2;
212 gxmr[ij] = -alfa_u_beta_v*v*ro2/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro2*v*tau_u_mu_v/den_i2;
213 gxmi[ij] = v*gama*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro3*v*tau_u_mu_v/den_i2;
214 gxnr[ij] = gama*ro4/den_r -2*(-alfa_u_beta_v_tau_u_mu_v+gama*ro2*nu)*ro4*nu/den_i2;
215 gxni[ij] = alfa_u_beta_v*ro3/den_r -2*(alfa_u_beta_v_nu+tau_u_mu_v_gama)*ro5*nu/den_i2;
216 }
217 }
218
grdgravmag3d_mirror_edges(gmt_grdfloat * grid,int nc,int i_data_start,int j_data_start,struct GRDREDPOL_CTRL * Ctrl)219 GMT_LOCAL void grdgravmag3d_mirror_edges (gmt_grdfloat *grid, int nc, int i_data_start, int j_data_start, struct GRDREDPOL_CTRL *Ctrl) {
220 /* This routine mirrors or replicates the West and East borders j_data_start times
221 and the South and North borders by i_data_start times.
222 nc is the total number of columns by which the grid is extended
223 Ctrl->S.n_columns & Ctrl->S.n_rows are the grid's original number of column/rows before extension */
224 int i, j, ins, isn, iss, jww, jwe, jee, jew, upper_nx, upper_ny;
225
226 /* First reflect about xmin and xmax, point symmetric about edge point */
227
228 upper_ny = Ctrl->S.n_rows+i_data_start;
229 for (j = 1; j <= j_data_start; j++) { /* COLUMNS */
230 jww = j_data_start-j; /* Minimum Outside xmin and approaching West border */
231 jee = Ctrl->S.n_columns + j_data_start + j-1; /* Minimum Outside xmax and approaching East border */
232 if (Ctrl->M.mirror) {
233 jwe = j_data_start+j; /* Minimum Inside xmin and approaching center */
234 jew = Ctrl->S.n_columns + j_data_start - j-1; /* Minimum Inside xmax and approaching center */
235 }
236 else {
237 jwe = j_data_start; /* West border */
238 jew = Ctrl->S.n_columns + j_data_start - 1; /* East border */
239 }
240 for (i = i_data_start; i < upper_ny; i++) { /* ROWS */
241 grid[ij0_data(Ctrl,i,jww)] = grid[ij0_data(Ctrl,i,jwe)]; /* West border */
242 grid[ij0_data(Ctrl,i,jee)] = grid[ij0_data(Ctrl,i,jew)]; /* East border */
243 }
244 }
245
246 /* Next, reflect about ymin and ymax. At the same time, since x has been reflected, we can use these vals */
247
248 upper_nx = Ctrl->S.n_columns + nc;
249 for (i = 0; i < i_data_start; i++) { /* ROWS */
250 iss = Ctrl->S.n_rows+i_data_start+i; /* Minimum Outside ymin and approaching South border */
251 if (Ctrl->M.mirror) {
252 ins = 2*i_data_start - i; /* Maximum Inside ymax and approaching North border */
253 isn = Ctrl->S.n_rows+i_data_start-2-i; /* Minimum Inside ymin and approaching center */
254 }
255 else {
256 ins = i_data_start; /* North border */
257 isn = Ctrl->S.n_rows+i_data_start-1; /* South border */
258 }
259 for (j = 0; j < upper_nx; j++) { /* COLUMNS */
260 grid[ij0_data(Ctrl,i,j)] = grid[ij0_data(Ctrl,ins,j)]; /* North border */
261 grid[ij0_data(Ctrl,iss,j)] = grid[ij0_data(Ctrl,isn,j)]; /* South border */
262 }
263 }
264 }
265
grdgravmag3d_tfpoeq(double * w,int m,int n,double * greel,double * gim,double * cosphi,double * sinphi,double * cospsi,double * sinpsi)266 GMT_LOCAL void grdgravmag3d_tfpoeq(double *w, int m, int n, double *greel, double *gim,
267 double *cosphi, double *sinphi, double *cospsi, double *sinpsi) {
268 /* Initialized data */
269
270 static int mkeep = -9999;
271 static int nkeep = -9999;
272
273 /* System generated locals */
274 int w_offset, greel_offset, gim_offset;
275 int i, k, l, k1, k2, m1, n1, l1, l2, ir, is, lr, ks, ky, lx, ir1, ir2, lx1, lrm, ksn;
276 static double co1, co2, si2, si1, c1c2, c1s2, arg, c2s1, s1s2, xmn, arg1, somi, somr;
277
278 /* THIS SUBROUTINE COMPUTES THE INVERSE FOURIER TRANSFORM OF A FILTER */
279 /* WHOSE RESULTANT INDEX IS 2 */
280
281 /* INPUT ARGUMENTS DESCRIPTION: */
282 /* - GREL => REAL PART OF THE FOURIER TRANSFORM OF THE FILTER */
283 /* - GIM => IMAGINARY PART OF THE FOURIER TRANSFORM */
284 /* - M/N => NUMBER OF LINES/COLUMN OF THE FILTER */
285
286 /* OUTPUT ARGUMENTS DESCRIPTION: */
287 /* - W => AN ARRAY WHICH CONTAINS THE FILTER COEFFICIENTS STORED */
288 /* COLUMNWISE */
289
290 /* REMARKS: */
291 /* M AND N MUST BE ODD */
292 /* GREL AND GIM MUST CONTAIN AT LEAST M*(N+1)/2 REAL*8 LOCATIONS */
293
294 /* Translated to C by f2c (and further massaged) from a routine of A. Galdeano */
295
296 /* Parameter adjustments */
297 gim_offset = 1 + m;
298 gim -= gim_offset;
299 greel_offset = 1 + m;
300 greel -= greel_offset;
301 w_offset = 1 + m;
302 w -= w_offset;
303
304 /* Function Body */
305 ky = (n + 1) / 2;
306 lx = (m + 1) / 2;
307 /* CALCUL DES COEF. DU FILTRE */
308 xmn = (double)(m * n);
309 if (m == mkeep) goto L2;
310 mkeep = m;
311 arg = TWO_PI / m;
312 for (i = 0; i < m; i++) {
313 arg1 = arg * i;
314 sinphi[i] = sin(arg1);
315 cosphi[i] = cos(arg1);
316 }
317 L2:
318 if (n == nkeep) goto L4;
319 nkeep = n;
320 arg = TWO_PI / n;
321 for (i = 0; i < n; i++) {
322 arg1 = arg * i;
323 sinpsi[i] = sin(arg1);
324 cospsi[i] = cos(arg1);
325 }
326 L4:
327 lx1 = lx + 1;
328 m1 = m + 1;
329 n1 = n + 1;
330 for (k1 = 1; k1 <= n; ++k1) {
331 k = k1 - ky;
332 k2 = n1 - k1;
333 for (l1 = lx; l1 <= m; ++l1) {
334 l = l1 - lx;
335 l2 = m1 - l1;
336 somr = 0.;
337 somi = 0.;
338 for (ir = lx1; ir <= m; ++ir) {
339 lr = ir - lx;
340 lrm = lr * l % m + 1;
341 somr += greel[ir + m] * cosphi[lrm - 1];
342 somi += gim[ir + m] * sinphi[lrm - 1];
343 }
344 for (is = 2; is <= ky; ++is) {
345 ks = is - 1;
346 ksn = ks * k % n + 1;
347 if (ksn <= 0) ksn += n;
348 co2 = cospsi[ksn - 1];
349 si2 = sinpsi[ksn - 1];
350 somr += greel[lx + is * m] * co2;
351 somi += gim[lx + is * m] * si2;
352 for (ir1 = lx1; ir1 <= m; ++ir1) {
353 lr = ir1 - lx;
354 ir2 = lx - lr;
355 lrm = lr * l % m + 1;
356 co1 = cosphi[lrm - 1];
357 si1 = sinphi[lrm - 1];
358 c1c2 = co1 * co2;
359 s1s2 = si1 * si2;
360 c1s2 = co1 * si2;
361 c2s1 = co2 * si1;
362 somr = somr + greel[ir1 + is * m] * (c1c2 - s1s2) + greel[ir2 + is * m] * (c1c2 + s1s2);
363 somi = somi + gim[ir1 + is * m] * (c1s2 + c2s1) + gim[ir2 + is * m] * (c1s2 - c2s1);
364 }
365 }
366 somr = somr + somr + greel[lx + m];
367 somi += somi;
368 w[l1 + k1 * m] = (somr + somi) / xmn;
369 if (l1 != lx)
370 w[l2 + k2 * m] = (somr - somi) / xmn;
371 }
372 }
373 }
374
375
grdgravmag3d_igrf10syn(struct GMT_CTRL * C,int isv,double date,int itype,double alt,double elong,double lat,double * out)376 GMT_LOCAL int grdgravmag3d_igrf10syn (struct GMT_CTRL *C, int isv, double date, int itype, double alt, double elong, double lat, double *out) {
377 /* This is a synthesis routine for the 10th generation IGRF as agreed
378 * in December 2004 by IAGA Working Group V-MOD. It is valid 1900.0 to
379 * 2010.0 inclusive. Values for dates from 1945.0 to 2000.0 inclusive are
380 * definitive, otherwise they are non-definitive.
381 * INPUT
382 * isv = 0 if main-field values are required
383 * isv = 1 if secular variation values are required
384 * date = year A.D. Must be greater than or equal to 1900.0 and
385 * less than or equal to 2015.0. Warning message is given
386 * for dates greater than 2010.0. Must be double precision.
387 * itype = 1 if geodetic (spheroid)
388 * itype = 2 if geocentric (sphere)
389 * alt = height in km above sea level if itype = 1
390 * = distance from centre of Earth in km if itype = 2 (>3485 km)
391 * lat = latitude (90-90)
392 * elong = east-longitude (0-360) -- it works also in [-180;+180]
393 * OUTPUT
394 * out[0] F = total intensity (nT) if isv = 0, rubbish if isv = 1
395 * out[1] H = horizontal intensity (nT)
396 * out[2] X = north component (nT) if isv = 0, nT/year if isv = 1
397 * out[3] Y = east component (nT) if isv = 0, nT/year if isv = 1
398 * out[4] Z = vertical component (nT) if isv = 0, nT/year if isv = 1
399 * out[5] D = declination
400 * out[6] I = inclination
401 *
402 * To get the other geomagnetic elements (D, I, H and secular
403 * variations dD, dH, dI and dF) use routines ptoc and ptocsv.
404 *
405 * Adapted from 8th generation version to include new maximum degree for
406 * main-field models for 2000.0 and onwards and use WGS84 spheroid instead
407 * of International Astronomical Union 1966 spheroid as recommended by IAGA
408 * in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
409 * radius (= 6371.0 km) but 6371.2 km is what is used in determining the
410 * coefficients. Adaptation by Susan Macmillan, August 2003 (for
411 * 9th generation) and December 2004.
412 *
413 * Joaquim Luis 1-MARS-2005
414 * Converted to C (with help of f2c, which explains the ugliness)
415 * 1995.0 coefficients as published in igrf9coeffs.xls and igrf10coeffs.xls
416 * used - (Kimmo Korhonen spotted 1 nT difference in 11 coefficients)
417 * Susan Macmillan July 2005 (PW update Oct 2006)
418 *
419 * Joaquim Luis 21-JAN-2010
420 * Updated for IGRF 11th generation
421 */
422
423 struct GRDREDPOL_IGRF {
424 double e_1[3450];
425 };
426 /* Initialized data */
427 static struct GRDREDPOL_IGRF equiv_22 = {
428 {-31543.,-2298., 5922., -677., 2905.,-1061., 924., 1121., /* g0 (1900) */
429 1022.,-1469., -330., 1256., 3., 572., 523., 876.,
430 628., 195., 660., -69., -361., -210., 134., -75.,
431 -184., 328., -210., 264., 53., 5., -33., -86.,
432 -124., -16., 3., 63., 61., -9., -11., 83.,
433 -217., 2., -58., -35., 59., 36., -90., -69.,
434 70., -55., -45., 0., -13., 34., -10., -41.,
435 -1., -21., 28., 18., -12., 6., -22., 11.,
436 8., 8., -4., -14., -9., 7., 1., -13.,
437 2., 5., -9., 16., 5., -5., 8., -18.,
438 8., 10., -20., 1., 14., -11., 5., 12.,
439 -3., 1., -2., -2., 8., 2., 10., -1.,
440 -2., -1., 2., -3., -4., 2., 2., 1.,
441 -5., 2., -2., 6., 6., -4., 4., 0.,
442 0., -2., 2., 4., 2., 0., 0., -6.,
443 -31464.,-2298., 5909., -728., 2928.,-1086., 1041., 1065., /* g1 (1905) */
444 1037.,-1494., -357., 1239., 34., 635., 480., 880.,
445 643., 203., 653., -77., -380., -201., 146., -65.,
446 -192., 328., -193., 259., 56., -1., -32., -93.,
447 -125., -26., 11., 62., 60., -7., -11., 86.,
448 -221., 4., -57., -32., 57., 32., -92., -67.,
449 70., -54., -46., 0., -14., 33., -11., -41.,
450 0., -20., 28., 18., -12., 6., -22., 11.,
451 8., 8., -4., -15., -9., 7., 1., -13.,
452 2., 5., -8., 16., 5., -5., 8., -18.,
453 8., 10., -20., 1., 14., -11., 5., 12.,
454 -3., 1., -2., -2., 8., 2., 10., 0.,
455 -2., -1., 2., -3., -4., 2., 2., 1.,
456 -5., 2., -2., 6., 6., -4., 4., 0.,
457 0., -2., 2., 4., 2., 0., 0., -6.,
458 -31354.,-2297., 5898., -769., 2948.,-1128., 1176., 1000., /* g2 (1910) */
459 1058.,-1524., -389., 1223., 62., 705., 425., 884.,
460 660., 211., 644., -90., -400., -189., 160., -55.,
461 -201., 327., -172., 253., 57., -9., -33., -102.,
462 -126., -38., 21., 62., 58., -5., -11., 89.,
463 -224., 5., -54., -29., 54., 28., -95., -65.,
464 71., -54., -47., 1., -14., 32., -12., -40.,
465 1., -19., 28., 18., -13., 6., -22., 11.,
466 8., 8., -4., -15., -9., 6., 1., -13.,
467 2., 5., -8., 16., 5., -5., 8., -18.,
468 8., 10., -20., 1., 14., -11., 5., 12.,
469 -3., 1., -2., -2., 8., 2., 10., 0.,
470 -2., -1., 2., -3., -4., 2., 2., 1.,
471 -5., 2., -2., 6., 6., -4., 4., 0.,
472 0., -2., 2., 4., 2., 0., 0., -6.,
473 -31212.,-2306., 5875., -802., 2956.,-1191., 1309., 917., /* g3 (1915) */
474 1084.,-1559., -421., 1212., 84., 778., 360., 887.,
475 678., 218., 631., -109., -416., -173., 178., -51.,
476 -211., 327., -148., 245., 58., -16., -34., -111.,
477 -126., -51., 32., 61., 57., -2., -10., 93.,
478 -228., 8., -51., -26., 49., 23., -98., -62.,
479 72., -54., -48., 2., -14., 31., -12., -38.,
480 2., -18., 28., 19., -15., 6., -22., 11.,
481 8., 8., -4., -15., -9., 6., 2., -13.,
482 3., 5., -8., 16., 6., -5., 8., -18.,
483 8., 10., -20., 1., 14., -11., 5., 12.,
484 -3., 1., -2., -2., 8., 2., 10., 0.,
485 -2., -1., 2., -3., -4., 2., 2., 1.,
486 -5., 2., -2., 6., 6., -4., 4., 0.,
487 0., -2., 1., 4., 2., 0., 0., -6.,
488 -31060.,-2317., 5845., -839., 2959.,-1259., 1407., 823., /* g4 (1920) */
489 1111.,-1600., -445., 1205., 103., 839., 293., 889.,
490 695., 220., 616., -134., -424., -153., 199., -57.,
491 -221., 326., -122., 236., 58., -23., -38., -119.,
492 -125., -62., 43., 61., 55., 0., -10., 96.,
493 -233., 11., -46., -22., 44., 18., -101., -57.,
494 73., -54., -49., 2., -14., 29., -13., -37.,
495 4., -16., 28., 19., -16., 6., -22., 11.,
496 7., 8., -3., -15., -9., 6., 2., -14.,
497 4., 5., -7., 17., 6., -5., 8., -19.,
498 8., 10., -20., 1., 14., -11., 5., 12.,
499 -3., 1., -2., -2., 9., 2., 10., 0.,
500 -2., -1., 2., -3., -4., 2., 2., 1.,
501 -5., 2., -2., 6., 6., -4., 4., 0.,
502 0., -2., 1., 4., 3., 0., 0., -6.,
503 -30926.,-2318., 5817., -893., 2969.,-1334., 1471., 728., /* g5 (1925) */
504 1140.,-1645., -462., 1202., 119., 881., 229., 891.,
505 711., 216., 601., -163., -426., -130., 217., -70.,
506 -230., 326., -96., 226., 58., -28., -44., -125.,
507 -122., -69., 51., 61., 54., 3., -9., 99.,
508 -238., 14., -40., -18., 39., 13., -103., -52.,
509 73., -54., -50., 3., -14., 27., -14., -35.,
510 5., -14., 29., 19., -17., 6., -21., 11.,
511 7., 8., -3., -15., -9., 6., 2., -14.,
512 4., 5., -7., 17., 7., -5., 8., -19.,
513 8., 10., -20., 1., 14., -11., 5., 12.,
514 -3., 1., -2., -2., 9., 2., 10., 0.,
515 -2., -1., 2., -3., -4., 2., 2., 1.,
516 -5., 2., -2., 6., 6., -4., 4., 0.,
517 0., -2., 1., 4., 3., 0., 0., -6.,
518 -30805.,-2316., 5808., -951., 2980.,-1424., 1517., 644., /* g6 (1930) */
519 1172.,-1692., -480., 1205., 133., 907., 166., 896.,
520 727., 205., 584., -195., -422., -109., 234., -90.,
521 -237., 327., -72., 218., 60., -32., -53., -131.,
522 -118., -74., 58., 60., 53., 4., -9., 102.,
523 -242., 19., -32., -16., 32., 8., -104., -46.,
524 74., -54., -51., 4., -15., 25., -14., -34.,
525 6., -12., 29., 18., -18., 6., -20., 11.,
526 7., 8., -3., -15., -9., 5., 2., -14.,
527 5., 5., -6., 18., 8., -5., 8., -19.,
528 8., 10., -20., 1., 14., -12., 5., 12.,
529 -3., 1., -2., -2., 9., 3., 10., 0.,
530 -2., -2., 2., -3., -4., 2., 2., 1.,
531 -5., 2., -2., 6., 6., -4., 4., 0.,
532 0., -2., 1., 4., 3., 0., 0., -6.,
533 -30715.,-2306., 5812.,-1018., 2984.,-1520., 1550., 586., /* g7 (1935) */
534 1206.,-1740., -494., 1215., 146., 918., 101., 903.,
535 744., 188., 565., -226., -415., -90., 249., -114.,
536 -241., 329., -51., 211., 64., -33., -64., -136.,
537 -115., -76., 64., 59., 53., 4., -8., 104.,
538 -246., 25., -25., -15., 25., 4., -106., -40.,
539 74., -53., -52., 4., -17., 23., -14., -33.,
540 7., -11., 29., 18., -19., 6., -19., 11.,
541 7., 8., -3., -15., -9., 5., 1., -15.,
542 6., 5., -6., 18., 8., -5., 7., -19.,
543 8., 10., -20., 1., 15., -12., 5., 11.,
544 -3., 1., -3., -2., 9., 3., 11., 0.,
545 -2., -2., 2., -3., -4., 2., 2., 1.,
546 -5., 2., -2., 6., 6., -4., 4., 0.,
547 0., -1., 2., 4., 3., 0., 0., -6.,
548 -30654.,-2292., 5821.,-1106., 2981.,-1614., 1566., 528., /* g8 (1940) */
549 1240.,-1790., -499., 1232., 163., 916., 43., 914.,
550 762., 169., 550., -252., -405., -72., 265., -141.,
551 -241., 334., -33., 208., 71., -33., -75., -141.,
552 -113., -76., 69., 57., 54., 4., -7., 105.,
553 -249., 33., -18., -15., 18., 0., -107., -33.,
554 74., -53., -52., 4., -18., 20., -14., -31.,
555 7., -9., 29., 17., -20., 5., -19., 11.,
556 7., 8., -3., -14., -10., 5., 1., -15.,
557 6., 5., -5., 19., 9., -5., 7., -19.,
558 8., 10., -21., 1., 15., -12., 5., 11.,
559 -3., 1., -3., -2., 9., 3., 11., 1.,
560 -2., -2., 2., -3., -4., 2., 2., 1.,
561 -5., 2., -2., 6., 6., -4., 4., 0.,
562 0., -1., 2., 4., 3., 0., 0., -6.,
563 -30594.,-2285., 5810.,-1244., 2990.,-1702., 1578., 477., /* g9 (1945) */
564 1282.,-1834., -499., 1255., 186., 913., -11., 944.,
565 776., 144., 544., -276., -421., -55., 304., -178.,
566 -253., 346., -12., 194., 95., -20., -67., -142.,
567 -119., -82., 82., 59., 57., 6., 6., 100.,
568 -246., 16., -25., -9., 21., -16., -104., -39.,
569 70., -40., -45., 0., -18., 0., 2., -29.,
570 6., -10., 28., 15., -17., 29., -22., 13.,
571 7., 12., -8., -21., -5., -12., 9., -7.,
572 7., 2., -10., 18., 7., 3., 2., -11.,
573 5., -21., -27., 1., 17., -11., 29., 3.,
574 -9., 16., 4., -3., 9., -4., 6., -3.,
575 1., -4., 8., -3., 11., 5., 1., 1.,
576 2., -20., -5., -1., -1., -6., 8., 6.,
577 -1., -4., -3., -2., 5., 0., -2., -2.,
578 -30554.,-2250., 5815.,-1341., 2998.,-1810., 1576., 381., /* ga (1950) */
579 1297.,-1889., -476., 1274., 206., 896., -46., 954.,
580 792., 136., 528., -278., -408., -37., 303., -210.,
581 -240., 349., 3., 211., 103., -20., -87., -147.,
582 -122., -76., 80., 54., 57., -1., 4., 99.,
583 -247., 33., -16., -12., 12., -12., -105., -30.,
584 65., -55., -35., 2., -17., 1., 0., -40.,
585 10., -7., 36., 5., -18., 19., -16., 22.,
586 15., 5., -4., -22., -1., 0., 11., -21.,
587 15., -8., -13., 17., 5., -4., -1., -17.,
588 3., -7., -24., -1., 19., -25., 12., 10.,
589 2., 5., 2., -5., 8., -2., 8., 3.,
590 -11., 8., -7., -8., 4., 13., -1., -2.,
591 13., -10., -4., 2., 4., -3., 12., 6.,
592 3., -3., 2., 6., 10., 11., 3., 8.,
593 -30500.,-2215., 5820.,-1440., 3003.,-1898., 1581., 291., /* gb (1955) */
594 1302.,-1944., -462., 1288., 216., 882., -83., 958.,
595 796., 133., 510., -274., -397., -23., 290., -230.,
596 -229., 360., 15., 230., 110., -23., -98., -152.,
597 -121., -69., 78., 47., 57., -9., 3., 96.,
598 -247., 48., -8., -16., 7., -12., -107., -24.,
599 65., -56., -50., 2., -24., 10., -4., -32.,
600 8., -11., 28., 9., -20., 18., -18., 11.,
601 9., 10., -6., -15., -14., 5., 6., -23.,
602 10., 3., -7., 23., 6., -4., 9., -13.,
603 4., 9., -11., -4., 12., -5., 7., 2.,
604 6., 4., -2., 1., 10., 2., 7., 2.,
605 -6., 5., 5., -3., -5., -4., -1., 0.,
606 2., -8., -3., -2., 7., -4., 4., 1.,
607 -2., -3., 6., 7., -2., -1., 0., -3.,
608 -30421.,-2169., 5791.,-1555., 3002.,-1967., 1590., 206., /* gc (1960) */
609 1302.,-1992., -414., 1289., 224., 878., -130., 957.,
610 800., 135., 504., -278., -394., 3., 269., -255.,
611 -222., 362., 16., 242., 125., -26., -117., -156.,
612 -114., -63., 81., 46., 58., -10., 1., 99.,
613 -237., 60., -1., -20., -2., -11., -113., -17.,
614 67., -56., -55., 5., -28., 15., -6., -32.,
615 7., -7., 23., 17., -18., 8., -17., 15.,
616 6., 11., -4., -14., -11., 7., 2., -18.,
617 10., 4., -5., 23., 10., 1., 8., -20.,
618 4., 6., -18., 0., 12., -9., 2., 1.,
619 0., 4., -3., -1., 9., -2., 8., 3.,
620 0., -1., 5., 1., -3., 4., 4., 1.,
621 0., 0., -1., 2., 4., -5., 6., 1.,
622 1., -1., -1., 6., 2., 0., 0., -7.,
623 -30334.,-2119., 5776.,-1662., 2997.,-2016., 1594., 114., /* gd (1965) */
624 1297.,-2038., -404., 1292., 240., 856., -165., 957.,
625 804., 148., 479., -269., -390., 13., 252., -269.,
626 -219., 358., 19., 254., 128., -31., -126., -157.,
627 -97., -62., 81., 45., 61., -11., 8., 100.,
628 -228., 68., 4., -32., 1., -8., -111., -7.,
629 75., -57., -61., 4., -27., 13., -2., -26.,
630 6., -6., 26., 13., -23., 1., -12., 13.,
631 5., 7., -4., -12., -14., 9., 0., -16.,
632 8., 4., -1., 24., 11., -3., 4., -17.,
633 8., 10., -22., 2., 15., -13., 7., 10.,
634 -4., -1., -5., -1., 10., 5., 10., 1.,
635 -4., -2., 1., -2., -3., 2., 2., 1.,
636 -5., 2., -2., 6., 4., -4., 4., 0.,
637 0., -2., 2., 3., 2., 0., 0., -6.,
638 -30220.,-2068., 5737.,-1781., 3000.,-2047., 1611., 25., /* ge (1970) */
639 1287.,-2091., -366., 1278., 251., 838., -196., 952.,
640 800., 167., 461., -266., -395., 26., 234., -279.,
641 -216., 359., 26., 262., 139., -42., -139., -160.,
642 -91., -56., 83., 43., 64., -12., 15., 100.,
643 -212., 72., 2., -37., 3., -6., -112., 1.,
644 72., -57., -70., 1., -27., 14., -4., -22.,
645 8., -2., 23., 13., -23., -2., -11., 14.,
646 6., 7., -2., -15., -13., 6., -3., -17.,
647 5., 6., 0., 21., 11., -6., 3., -16.,
648 8., 10., -21., 2., 16., -12., 6., 10.,
649 -4., -1., -5., 0., 10., 3., 11., 1.,
650 -2., -1., 1., -3., -3., 1., 2., 1.,
651 -5., 3., -1., 4., 6., -4., 4., 0.,
652 1., -1., 0., 3., 3., 1., -1., -4.,
653 -30100.,-2013., 5675.,-1902., 3010.,-2067., 1632., -68., /* gf (1975) */
654 1276.,-2144., -333., 1260., 262., 830., -223., 946.,
655 791., 191., 438., -265., -405., 39., 216., -288.,
656 -218., 356., 31., 264., 148., -59., -152., -159.,
657 -83., -49., 88., 45., 66., -13., 28., 99.,
658 -198., 75., 1., -41., 6., -4., -111., 11.,
659 71., -56., -77., 1., -26., 16., -5., -14.,
660 10., 0., 22., 12., -23., -5., -12., 14.,
661 6., 6., -1., -16., -12., 4., -8., -19.,
662 4., 6., 0., 18., 10., -10., 1., -17.,
663 7., 10., -21., 2., 16., -12., 7., 10.,
664 -4., -1., -5., -1., 10., 4., 11., 1.,
665 -3., -2., 1., -3., -3., 1., 2., 1.,
666 -5., 3., -2., 4., 5., -4., 4., -1.,
667 1., -1., 0., 3., 3., 1., -1., -5.,
668 -29992.,-1956., 5604.,-1997., 3027.,-2129., 1663., -200., /* gg (1980) */
669 1281.,-2180., -336., 1251., 271., 833., -252., 938.,
670 782., 212., 398., -257., -419., 53., 199., -297.,
671 -218., 357., 46., 261., 150., -74., -151., -162.,
672 -78., -48., 92., 48., 66., -15., 42., 93.,
673 -192., 71., 4., -43., 14., -2., -108., 17.,
674 72., -59., -82., 2., -27., 21., -5., -12.,
675 16., 1., 18., 11., -23., -2., -10., 18.,
676 6., 7., 0., -18., -11., 4., -7., -22.,
677 4., 9., 3., 16., 6., -13., -1., -15.,
678 5., 10., -21., 1., 16., -12., 9., 9.,
679 -5., -3., -6., -1., 9., 7., 10., 2.,
680 -6., -5., 2., -4., -4., 1., 2., 0.,
681 -5., 3., -2., 6., 5., -4., 3., 0.,
682 1., -1., 2., 4., 3., 0., 0., -6.,
683 -29873.,-1905., 5500.,-2072., 3044.,-2197., 1687., -306., /* gi (1985) */
684 1296.,-2208., -310., 1247., 284., 829., -297., 936.,
685 780., 232., 361., -249., -424., 69., 170., -297.,
686 -214., 355., 47., 253., 150., -93., -154., -164.,
687 -75., -46., 95., 53., 65., -16., 51., 88.,
688 -185., 69., 4., -48., 16., -1., -102., 21.,
689 74., -62., -83., 3., -27., 24., -2., -6.,
690 20., 4., 17., 10., -23., 0., -7., 21.,
691 6., 8., 0., -19., -11., 5., -9., -23.,
692 4., 11., 4., 14., 4., -15., -4., -11.,
693 5., 10., -21., 1., 15., -12., 9., 9.,
694 -6., -3., -6., -1., 9., 7., 9., 1.,
695 -7., -5., 2., -4., -4., 1., 3., 0.,
696 -5., 3., -2., 6., 5., -4., 3., 0.,
697 1., -1., 2., 4., 3., 0., 0., -6.,
698 -29775.,-1848., 5406.,-2131., 3059.,-2279., 1686., -373., /* gj (1990) */
699 1314.,-2239., -284., 1248., 293., 802., -352., 939.,
700 780., 247., 325., -240., -423., 84., 141., -299.,
701 -214., 353., 46., 245., 154., -109., -153., -165.,
702 -69., -36., 97., 61., 65., -16., 59., 82.,
703 -178., 69., 3., -52., 18., 1., -96., 24.,
704 77., -64., -80., 2., -26., 26., 0., -1.,
705 21., 5., 17., 9., -23., 0., -4., 23.,
706 5., 10., -1., -19., -10., 6., -12., -22.,
707 3., 12., 4., 12., 2., -16., -6., -10.,
708 4., 9., -20., 1., 15., -12., 11., 9.,
709 -7., -4., -7., -2., 9., 7., 8., 1.,
710 -7., -6., 2., -3., -4., 2., 2., 1.,
711 -5., 3., -2., 6., 4., -4., 3., 0.,
712 1., -2., 3., 3., 3., -1., 0., -6.,
713 -29692.,-1784., 5306.,-2200., 3070.,-2366., 1681., -413., /* gk (1995) */
714 1335.,-2267., -262., 1249., 302., 759., -427., 940.,
715 780., 262., 290., -236., -418., 97., 122., -306.,
716 -214., 352., 46., 235., 165., -118., -143., -166.,
717 -55., -17., 107., 68., 67., -17., 68., 72.,
718 -170., 67., -1., -58., 19., 1., -93., 36.,
719 77., -72., -69., 1., -25., 28., 4., 5.,
720 24., 4., 17., 8., -24., -2., -6., 25.,
721 6., 11., -6., -21., -9., 8., -14., -23.,
722 9., 15., 6., 11., -5., -16., -7., -4.,
723 4., 9., -20., 3., 15., -10., 12., 8.,
724 -6., -8., -8., -1., 8., 10., 5., -2.,
725 -8., -8., 3., -3., -6., 1., 2., 0.,
726 -4., 4., -1., 5., 4., -5., 2., -1.,
727 2., -2., 5., 1., 1., -2., 0., -7.,
728 0., 0., 0., 0., 0., 0., 0., 0.,
729 0., 0., 0., 0., 0., 0., 0., 0.,
730 0., 0., 0., 0., 0., 0., 0., 0.,
731 0., 0., 0., 0., 0., 0., 0., 0.,
732 0., 0., 0., 0., 0., 0., 0., 0.,
733 0., 0., 0., 0., 0., 0., 0., 0.,
734 0., 0., 0., 0., 0., 0., 0., 0.,
735 0., 0., 0., 0., 0., 0., 0., 0.,
736 0., 0., 0., 0., 0., 0., 0., 0.,
737 0., 0., 0.,
738 -29619.4,-1728.2, 5186.1,-2267.7, 3068.4,-2481.6, 1670.9, /* gl (2000) */
739 -458.0, 1339.6,-2288.0, -227.6, 1252.1, 293.4, 714.5,
740 -491.1, 932.3, 786.8, 272.6, 250.0, -231.9, -403.0,
741 119.8, 111.3, -303.8, -218.8, 351.4, 43.8, 222.3,
742 171.9, -130.4, -133.1, -168.6, -39.3, -12.9, 106.3,
743 72.3, 68.2, -17.4, 74.2, 63.7, -160.9, 65.1,
744 -5.9, -61.2, 16.9, 0.7, -90.4, 43.8, 79.0,
745 -74.0, -64.6, 0.0, -24.2, 33.3, 6.2, 9.1,
746 24.0, 6.9, 14.8, 7.3, -25.4, -1.2, -5.8,
747 24.4, 6.6, 11.9, -9.2, -21.5, -7.9, 8.5,
748 -16.6, -21.5, 9.1, 15.5, 7.0, 8.9, -7.9,
749 -14.9, -7.0, -2.1, 5.0, 9.4, -19.7, 3.0,
750 13.4, -8.4, 12.5, 6.3, -6.2, -8.9, -8.4,
751 -1.5, 8.4, 9.3, 3.8, -4.3, -8.2, -8.2,
752 4.8, -2.6, -6.0, 1.7, 1.7, 0.0, -3.1,
753 4.0, -0.5, 4.9, 3.7, -5.9, 1.0, -1.2,
754 2.0, -2.9, 4.2, 0.2, 0.3, -2.2, -1.1,
755 -7.4, 2.7, -1.7, 0.1, -1.9, 1.3, 1.5,
756 -0.9, -0.1, -2.6, 0.1, 0.9, -0.7, -0.7,
757 0.7, -2.8, 1.7, -0.9, 0.1, -1.2, 1.2,
758 -1.9, 4.0, -0.9, -2.2, -0.3, -0.4, 0.2,
759 0.3, 0.9, 2.5, -0.2, -2.6, 0.9, 0.7,
760 -0.5, 0.3, 0.3, 0.0, -0.3, 0.0, -0.4,
761 0.3, -0.1, -0.9, -0.2, -0.4, -0.4, 0.8,
762 -0.2, -0.9, -0.9, 0.3, 0.2, 0.1, 1.8,
763 -0.4, -0.4, 1.3, -1.0, -0.4, -0.1, 0.7,
764 0.7, -0.4, 0.3, 0.3, 0.6, -0.1, 0.3,
765 0.4, -0.2, 0.0, -0.5, 0.1, -0.9,
766 -29554.63, -1669.05, 5077.99, -2337.24, 3047.69, -2594.50, 1657.76, /* gm (2005) */
767 -515.43, 1336.30, -2305.83, -198.86, 1246.39, 269.72, 672.51,
768 -524.72, 920.55, 797.96, 282.07, 210.65, -225.23, -379.86,
769 145.15, 100.00, -305.36, -227.00, 354.41, 42.72, 208.95,
770 180.25, -136.54, -123.45, -168.05, -19.57, -13.55, 103.85,
771 73.60, 69.56, -20.33, 76.74, 54.75, -151.34, 63.63,
772 -14.58, -63.53, 14.58, 0.24, -86.36, 50.94, 79.88,
773 -74.46, -61.14, -1.65, -22.57, 38.73, 6.82, 12.30,
774 25.35, 9.37, 10.93, 5.42, -26.32, 1.94, -4.64,
775 24.80, 7.62, 11.20, -11.73, -20.88, -6.88, 9.83,
776 -18.11, -19.71, 10.17, 16.22, 9.36, 7.61, -11.25,
777 -12.76, -4.87, -0.06, 5.58, 9.76, -20.11, 3.58,
778 12.69, -6.94, 12.67, 5.01, -6.72, -10.76, -8.16,
779 -1.25, 8.10, 8.76, 2.92, -6.66, -7.73, -9.22,
780 6.01, -2.17, -6.12, 2.19, 1.42, 0.10, -2.35,
781 4.46, -0.15, 4.76, 3.06, -6.58, 0.29, -1.01,
782 2.06, -3.47, 3.77, -0.86, -0.21, -2.31, -2.09,
783 -7.93, 2.95, -1.60, 0.26, -1.88, 1.44, 1.44,
784 -0.77, -0.31, -2.27, 0.29, 0.90, -0.79, -0.58,
785 0.53, -2.69, 1.80, -1.08, 0.16, -1.58, 0.96,
786 -1.90, 3.99, -1.39, -2.15, -0.29, -0.55, 0.21,
787 0.23, 0.89, 2.38, -0.38, -2.63, 0.96, 0.61,
788 -0.30, 0.40, 0.46, 0.01, -0.35, 0.02, -0.36,
789 0.28, 0.08, -0.87, -0.49, -0.34, -0.08, 0.88,
790 -0.16, -0.88, -0.76, 0.30, 0.33, 0.28, 1.72,
791 -0.43, -0.54, 1.18, -1.07, -0.37, -0.04, 0.75,
792 0.63, -0.26, 0.21, 0.35, 0.53, -0.05, 0.38,
793 0.41, -0.22, -0.10, -0.57, -0.18, -0.82,
794 -29496.57, -1586.42, 4944.26, -2396.06, 3026.34, -2708.54, /* gp (2010) */
795 1668.17, -575.73, 1339.85, -2326.54, -160.40, 1232.10,
796 251.75, 633.73, -537.03, 912.66, 808.97, 286.48,
797 166.58, -211.03, -356.83, 164.46, 89.40, -309.72,
798 -230.87, 357.29, 44.58, 200.26, 189.01, -141.05,
799 -118.06, -163.17, -0.01, -8.03, 101.04, 72.78,
800 68.69, -20.90, 75.92, 44.18, -141.40, 61.54,
801 -22.83, -66.26, 13.10, 3.02, -78.09, 55.40,
802 80.44, -75.00, -57.80, -4.55, -21.20, 45.24,
803 6.54, 14.00, 24.96, 10.46, 7.03, 1.64,
804 -27.61, 4.92, -3.28, 24.41, 8.21, 10.84,
805 -14.50, -20.03, -5.59, 11.83, -19.34, -17.41,
806 11.61, 16.71, 10.85, 6.96, -14.05, -10.74,
807 -3.54, 1.64, 5.50, 9.45, -20.54, 3.45,
808 11.51, -5.27, 12.75, 3.13, -7.14, -12.38,
809 -7.42, -0.76, 7.97, 8.43, 2.14, -8.42,
810 -6.08, -10.08, 7.01, -1.94, -6.24, 2.73,
811 0.89, -0.10, -1.07, 4.71, -0.16, 4.44,
812 2.45, -7.22, -0.33, -0.96, 2.13, -3.95,
813 3.09, -1.99, -1.03, -1.97, -2.80, -8.31,
814 3.05, -1.48, 0.13, -2.03, 1.67, 1.65,
815 -0.66, -0.51, -1.76, 0.54, 0.85, -0.79,
816 -0.39, 0.37, -2.51, 1.79, -1.27, 0.12,
817 -2.11, 0.75, -1.94, 3.75, -1.86, -2.12,
818 -0.21, -0.87, 0.30, 0.27, 1.04, 2.13,
819 -0.63, -2.49, 0.95, 0.49, -0.11, 0.59,
820 0.52, 0.00, -0.39, 0.13, -0.37, 0.27,
821 0.21, -0.86, -0.77, -0.23, 0.04, 0.87,
822 -0.09, -0.89, -0.87, 0.31, 0.30, 0.42,
823 1.66, -0.45, -0.59, 1.08, -1.14, -0.31,
824 -0.07, 0.78, 0.54, -0.18, 0.10, 0.38,
825 0.49, 0.02, 0.44, 0.42, -0.25, -0.26,
826 -0.53, -0.26, -0.79,
827 -29442.0, -1501.0, 4797.1,-2445.1, 3012.9,-2845.6, 1676.7, /* gq (2015) */
828 -641.9, 1350.7, -2352.3, -115.3, 1225.6, 244.9, 582.0,
829 -538.4, 907.6, 813.7, 283.3, 120.4, -188.7, -334.9,
830 180.9, 70.4, -329.5, -232.6, 360.1, 47.3, 192.4,
831 197.0, -140.9, -119.3, -157.5, 16.0, 4.1, 100.2,
832 70.0, 67.7, -20.8, 72.7, 33.2, -129.9, 58.9,
833 -28.9, -66.7, 13.2, 7.3, -70.9, 62.6, 81.6,
834 -76.1, -54.1, -6.8, -19.5, 51.8, 5.7, 15.0,
835 24.4, 9.4, 3.4, -2.8, -27.4, 6.8, -2.2,
836 24.2, 8.8, 10.1, -16.9, -18.3, -3.2, 13.3,
837 -20.6, -14.6, 13.4, 16.2, 11.7, 5.7, -15.9,
838 -9.1, -2.0, 2.1, 5.4, 8.8, -21.6, 3.1,
839 10.8, -3.3, 11.8, 0.7, -6.8, -13.3, -6.9,
840 -0.1, 7.8, 8.7, 1.0, -9.1, -4.0, -10.5,
841 8.4, -1.9, -6.3, 3.2, 0.1, -0.4, 0.5,
842 4.6, -0.5, 4.4, 1.8, -7.9, -0.7, -0.6,
843 2.1, -4.2, 2.4, -2.8, -1.8, -1.2, -3.6,
844 -8.7, 3.1, -1.5, -0.1, -2.3, 2.0, 2.0,
845 -0.7, -0.8, -1.1, 0.6, 0.8, -0.7, -0.2,
846 0.2, -2.2, 1.7, -1.4, -0.2, -2.5, 0.4,
847 -2.0, 3.5, -2.4, -1.9, -0.2, -1.1, 0.4,
848 0.4, 1.2, 1.9, -0.8, -2.2, 0.9, 0.3,
849 0.1, 0.7, 0.5, -0.1, -0.3, 0.3, -0.4,
850 0.2, 0.2, -0.9, -0.9, -0.1, 0.0, 0.7,
851 0.0, -0.9, -0.9, 0.4, 0.4, 0.5, 1.6,
852 -0.5, -0.5, 1.0, -1.2, -0.2, -0.1, 0.8,
853 0.4, -0.1, -0.1, 0.3, 0.4, 0.1, 0.5,
854 0.5, -0.3, -0.4, -0.4, -0.3, -0.8,
855 10.3, 18.1, -26.6, -8.7, -3.3, -27.4, 2.1, /* sv (2015) */
856 -14.1, 3.4, -5.5, 8.2, -0.7, -0.4, -10.1,
857 1.8, -0.7, 0.2, -1.3, -9.1, 5.3, 4.1,
858 2.9, -4.3, -5.2, -0.2, 0.5, 0.6, -1.3,
859 1.7, -0.1, -1.2, 1.4, 3.4, 3.9, 0.0,
860 -0.3, -0.1, 0.0, -0.7, -2.1, 2.1, -0.7,
861 -1.2, 0.2, 0.3, 0.9, 1.6, 1.0, 0.3,
862 -0.2, 0.8, -0.5, 0.4, 1.3, -0.2, 0.1,
863 -0.3, -0.6, -0.6, -0.8, 0.1, 0.2, -0.2,
864 0.2, 0.0, -0.3, -0.6, 0.3, 0.5, 0.1,
865 -0.2, 0.5, 0.4, -0.2, 0.1, -0.3, -0.4,
866 0.3, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0,
867 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
868 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
869 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
870 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
871 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
872 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
873 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
874 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
875 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
876 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
877 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
878 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
879 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
880 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
881 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
882 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}
883 };
884 #define gh ((double *)&equiv_22)
885
886 int i, j, k, l, m, n, ll, lm, kmx, nmx, nc;
887 double cd, cl[13], tc, ct, sd, fn = 0.0, gn = 0.0, fm, sl[13];
888 double rr, st, one, gmm, rho, two, three, ratio;
889 double p[105], q[105], r, t, a2, b2;
890 double H, F, X = 0, Y = 0, Z = 0, dec, dip;
891
892 if (date < 1900.0 || date > 2020.0) {
893 GMT_Report (C->parent, GMT_MSG_ERROR, "Your date (%g) is outside valid extrapolated range for IGRF (1900-2020)\n", date);
894 return (true);
895 }
896
897 if (date < 2015.) {
898 t = 0.2 * (date - 1900.);
899 ll = (int) t;
900 one = (double) ll;
901 t -= one;
902 if (date < 1995.) {
903 nmx = 10;
904 nc = nmx * (nmx + 2);
905 ll = nc * ll;
906 kmx = (nmx + 1) * (nmx + 2) / 2;
907 } else {
908 nmx = 13;
909 nc = nmx * (nmx + 2);
910 ll = (int) ((date - 1995.) * .2);
911 ll = nc * ll + 2280 /* 2280, position of first coeff of 1995 */;
912 kmx = (nmx + 1) * (nmx + 2) / 2;
913 }
914 tc = 1. - t;
915 if (isv == 1) {
916 tc = -.2;
917 t = .2;
918 }
919 }
920 else {
921 t = date - 2015.;
922 tc = 1.;
923 if (isv == 1) {
924 t = 1.;
925 tc = 0.;
926 }
927 ll = 3060; /* nth position corresponding to first coeff of 2015 */
928 nmx = 13;
929 nc = nmx * (nmx + 2);
930 kmx = (nmx + 1) * (nmx + 2) / 2;
931 }
932 r = alt;
933 sincosd (90.0 - lat, &st, &ct);
934 sincosd (elong, &(sl[0]), &(cl[0]));
935 cd = 1.;
936 sd = 0.;
937 l = 1;
938 m = 1;
939 n = 0;
940 if (itype == 1) { /* conversion from geodetic to geocentric coordinates (using the WGS84 spheroid) */
941 a2 = 40680631.6;
942 b2 = 40408296.0;
943 one = a2 * st * st;
944 two = b2 * ct * ct;
945 three = one + two;
946 rho = sqrt(three);
947 r = sqrt(alt * (alt + rho * 2.) + (a2 * one + b2 * two) / three);
948 cd = (alt + rho) / r;
949 sd = (a2 - b2) / rho * ct * st / r;
950 one = ct;
951 ct = ct * cd - st * sd;
952 st = st * cd + one * sd;
953 }
954 ratio = 6371.2 / r;
955 rr = ratio * ratio;
956
957 /* computation of Schmidt quasi-normal coefficients p and x(=q) */
958
959 p[0] = 1.;
960 p[2] = st;
961 q[0] = 0.;
962 q[2] = ct;
963 for (k = 2; k <= kmx; ++k) {
964 if (n < m) {
965 m = 0;
966 n++;
967 rr *= ratio;
968 fn = (double) n;
969 gn = (double) (n - 1);
970 }
971 fm = (double) m;
972 if (k != 3) {
973 if (m == n) {
974 one = sqrt(1. - .5 / fm);
975 j = k - n - 1;
976 p[k-1] = one * st * p[j-1];
977 q[k-1] = one * (st * q[j-1] + ct * p[j-1]);
978 cl[m-1] = cl[m-2] * cl[0] - sl[m-2] * sl[0];
979 sl[m-1] = sl[m-2] * cl[0] + cl[m-2] * sl[0];
980 }
981 else {
982 gmm = (double) (m * m);
983 one = sqrt(fn * fn - gmm);
984 two = sqrt(gn * gn - gmm) / one;
985 three = (fn + gn) / one;
986 i = k - n;
987 j = i - n + 1;
988 p[k-1] = three * ct * p[i-1] - two * p[j-1];
989 q[k-1] = three * (ct * q[i-1] - st * p[i-1]) - two * q[j-1];
990 }
991 }
992
993 /* synthesis of x, y and z in geocentric coordinates */
994
995 lm = ll + l;
996 one = (tc * gh[lm-1] + t * gh[lm+nc-1]) * rr;
997 if (m == 0) {
998 X += one * q[k-1];
999 Z -= (fn + 1.) * one * p[k-1];
1000 l++;
1001 }
1002 else {
1003 two = (tc * gh[lm] + t * gh[lm+nc]) * rr;
1004 three = one * cl[m-1] + two * sl[m - 1];
1005 X += three * q[k-1];
1006 Z -= (fn + 1.) * three * p[k-1];
1007 if (st != 0.)
1008 Y += (one * sl[m-1] - two * cl[m-1]) * fm * p[k-1] / st;
1009 else
1010 Y += (one * sl[m-1] - two * cl[m-1]) * q[k-1] * ct;
1011 l += 2;
1012 }
1013 m++;
1014 }
1015
1016 /* conversion to coordinate system specified by itype */
1017 one = X;
1018 X = X * cd + Z * sd;
1019 Z = Z * cd - one * sd;
1020 H = sqrt(X*X + Y*Y);
1021 F = sqrt(H*H + Z*Z);
1022 dec = atan2d(Y,X); dip = atan2d(Z,H);
1023 out[0] = F; out[1] = H;
1024 out[2] = X; out[3] = Y;
1025 out[4] = Z;
1026 out[5] = dec; out[6] = dip;
1027
1028 return (GMT_NOERROR);
1029 }
1030
usage(struct GMTAPI_CTRL * API,int level)1031 static int usage (struct GMTAPI_CTRL *API, int level) {
1032 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
1033 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
1034 GMT_Usage (API, 0, "usage: %s %s -G<rtp_grdfile> [-C<dec>/<dip>] [-Ei|d<grid>] [-F<m>/<n>] "
1035 "[-Mm|r] [-N] [-T<year>] [%s] [%s] [-W<win_width>] [-Z<filterfile>] [%s] [%s]\n",
1036 name, GMT_INGRID, GMT_Rgeo_OPT, GMT_V_OPT, GMT_n_OPT, GMT_PAR_OPT);
1037
1038 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
1039
1040 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
1041 gmt_ingrid_syntax (API, 0, "Name of grid with the magnetic anomaly");
1042 GMT_Usage (API, 1, "\n-G<rtp_grdfile>");
1043 GMT_Usage (API, -2, "Set filename for output grid with the RTP solution.");
1044 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
1045 GMT_Usage (API, 1, "\n-C<dec>/<dip>");
1046 GMT_Usage (API, -2, "Set <dec>/<dip> and uses this constant values in the RTP procedure.");
1047 GMT_Usage (API, 1, "\n-Ei|d<grid>");
1048 GMT_Usage (API, -2, "Specify input grids for inclination and/or declination:");
1049 GMT_Usage (API, 3, "i: Append file name for magnetization inclination [default uses IGRF].");
1050 GMT_Usage (API, 3, "d: Append file name for the magnetization declination [default uses IGRF].");
1051 GMT_Usage (API, 1, "\n-F<m>/<n>");
1052 GMT_Usage (API, -2, "Set <m>/<n> filter widths [25/25].");
1053 GMT_Usage (API, 1, "\n-Mm|r");
1054 GMT_Usage (API, -2, "Specify optional boundary conditions [Default is zero padding]:");
1055 GMT_Usage (API, 3, "m: mirror boundary condition.");
1056 GMT_Usage (API, 3, "r: Replicate edges.");
1057 GMT_Usage (API, 1, "\n-N Do NOT use Taylor expansion.");
1058 GMT_Option (API, "R");
1059 GMT_Usage (API, 1, "\n-T<year>");
1060 GMT_Usage (API, -2, "Set year used by the IGRF routine to compute the various DECs & DIPs [default is 2000].");
1061 GMT_Option (API, "V");
1062 GMT_Usage (API, 1, "\n-W<win_width>");
1063 GMT_Usage (API, -2, "Set window width in degrees [5].");
1064 GMT_Usage (API, 1, "\n-Z<filterfile>");
1065 GMT_Usage (API, -2, "Write filter file <filterfile> to disk.");
1066 GMT_Option (API, "n,.");
1067
1068 return (GMT_MODULE_USAGE);
1069 }
1070
parse(struct GMT_CTRL * GMT,struct GRDREDPOL_CTRL * Ctrl,struct GMT_OPTION * options)1071 static int parse (struct GMT_CTRL *GMT, struct GRDREDPOL_CTRL *Ctrl, struct GMT_OPTION *options) {
1072 /* This parses the options provided to grdredpol and sets parameters in Ctrl.
1073 * Note Ctrl has already been initialized and non-zero default values set.
1074 * Any GMT common options will override values set previously by other commands.
1075 * It also replaces any file names specified as input or output with the data ID
1076 * returned when registering these sources/destinations with the API.
1077 */
1078
1079 unsigned int n_errors = 0, n_files = 0, pos = 0;
1080 int j;
1081 char p[GMT_LEN256] = {""};
1082 struct GMT_OPTION *opt = NULL;
1083 struct GMTAPI_CTRL *API = GMT->parent;
1084
1085 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
1086
1087 switch (opt->option) {
1088 case '<': /* Input file (only one is accepted) */
1089 Ctrl->In.active = true;
1090 if (n_files++ == 0) Ctrl->In.file = strdup (opt->arg);
1091 break;
1092
1093 /* Processes program-specific parameters */
1094
1095 case 'C':
1096 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
1097 Ctrl->C.active = true;
1098 sscanf (opt->arg, "%lf/%lf", &Ctrl->C.dec, &Ctrl->C.dip);
1099 Ctrl->C.dec *= D2R;
1100 Ctrl->C.dip *= D2R;
1101 Ctrl->C.const_f = true;
1102 Ctrl->C.use_igrf = false;
1103 break;
1104 case 'E': /* -Ei<dip_grid> -Ee<dec_grid> */
1105 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
1106 Ctrl->E.active = true;
1107 Ctrl->C.use_igrf = false;
1108 Ctrl->E.dip_grd_only = true;
1109
1110 if (opt->arg[0] == 'i') /* Will fail if old syntax is used and grid name starts with an 'i' */
1111 Ctrl->E.dipfile = strdup (&opt->arg[1]);
1112 else if (opt->arg[0] == 'd') { /* Will fail if old syntax is used and grid name starts with an 'd' */
1113 Ctrl->E.decfile = strdup (&opt->arg[1]);
1114 Ctrl->E.dip_dec_grd = true;
1115 Ctrl->E.dip_grd_only = false;
1116 }
1117 else { /* Old syntax -E<dip_grd>[/<dec_grd>] */
1118 j = 0;
1119 while (gmt_strtok (opt->arg, "/", &pos, p)) {
1120 switch (j) {
1121 case 0:
1122 Ctrl->E.dipfile = strdup (p);
1123 break;
1124 case 1:
1125 Ctrl->E.decfile = strdup (p);
1126 Ctrl->E.dip_grd_only = false;
1127 Ctrl->E.dip_dec_grd = true;
1128 break;
1129 default:
1130 GMT_Report (API, GMT_MSG_ERROR, "Syntax error using option -E\n");
1131 n_errors++;
1132 break;
1133 }
1134 j++;
1135 }
1136 }
1137 break;
1138 case 'F':
1139 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
1140 Ctrl->F.active = true;
1141 j = sscanf (opt->arg, "%d/%d", &Ctrl->F.ncoef_row, &Ctrl->F.ncoef_col);
1142 if (j == 1) Ctrl->F.compute_n = true; /* Case of only one filter dimension was given */
1143 if (Ctrl->F.ncoef_row %2 != 1 || Ctrl->F.ncoef_col %2 != 1) {
1144 GMT_Report (API, GMT_MSG_ERROR, "Number of filter coefficients must be odd\n");
1145 n_errors++;
1146 }
1147 if (Ctrl->F.ncoef_row < 5 || Ctrl->F.ncoef_col < 5) {
1148 GMT_Report (API, GMT_MSG_ERROR, "That was a ridiculous number of filter coefficients\n");
1149 n_errors++;
1150 }
1151 break;
1152 case 'G':
1153 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
1154 Ctrl->G.active = true;
1155 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
1156 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
1157 break;
1158 case 'M':
1159 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
1160 Ctrl->M.pad_zero = false;
1161 for (j = 0; opt->arg[j]; j++) {
1162 if (opt->arg[j] == 'm')
1163 Ctrl->M.mirror = true;
1164 else if (opt->arg[j] == 'r')
1165 Ctrl->M.mirror = false;
1166 else {
1167 GMT_Report (API, GMT_MSG_ERROR, "Failure while using option -M (option ignored)\n");
1168 Ctrl->M.pad_zero = true;
1169 }
1170 }
1171 break;
1172 case 'N':
1173 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
1174 Ctrl->N.active = false;
1175 break;
1176 case 'T':
1177 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
1178 Ctrl->T.active = false;
1179 sscanf (opt->arg, "%lf", &Ctrl->T.year);
1180 break;
1181 case 'W':
1182 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
1183 Ctrl->W.active = false;
1184 sscanf (opt->arg, "%lf", &Ctrl->W.wid);
1185 break;
1186 case 'Z':
1187 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
1188 Ctrl->Z.active = true;
1189 Ctrl->Z.file = strdup (opt->arg);
1190 break;
1191 default: /* Report bad options */
1192 n_errors += gmt_default_error (GMT, opt->option);
1193 break;
1194 }
1195 }
1196
1197 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input file\n");
1198 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify output file\n");
1199
1200 if (Ctrl->C.const_f && Ctrl->C.use_igrf) {
1201 GMT_Report (API, GMT_MSG_ERROR, "-E option overrides -C\n");
1202 Ctrl->C.const_f = false;
1203 }
1204
1205 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
1206 }
1207
1208 #define bailout(code) {gmt_M_free_options (mode); return (code);}
1209 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
1210
GMT_grdredpol(void * V_API,int mode,void * args)1211 EXTERN_MSC int GMT_grdredpol (void *V_API, int mode, void *args) {
1212
1213 bool wrote_one = false;
1214 unsigned int i, j, row, col, nx_new, ny_new, one_or_zero, m21, n21, i2, j2;
1215 unsigned int k, l, i3, n_jlon, n_jlat, n_coef;
1216 int error = 0;
1217 uint64_t ij, jj;
1218 double tmp_d, sloni, slati, slonf, slatf, slonm, slatm;
1219 double *ftlon = NULL, *ftlat = NULL, *gxr = NULL, *gxi = NULL, *fxr = NULL;
1220 double *gxar = NULL, *gxai = NULL, *gxbr = NULL, *gxbi = NULL, *gxgr = NULL;
1221 double *gxgi = NULL, *fxar = NULL, *fxbr = NULL, *fxgr = NULL, *fix = NULL;
1222 double *gxtr = NULL, *gxti = NULL, *gxmr = NULL, *gxmi = NULL, *gxnr = NULL;
1223 double *gxni = NULL, *fxtr = NULL, *fxmr = NULL, *fxnr = NULL;
1224 double *cosphi = NULL, *sinphi = NULL, *cospsi = NULL, *sinpsi = NULL;
1225 double fi, psi, alfa = 0, beta = 0, gama = 0, r, s, u, v;
1226 double alfa1, beta1, gama1, da = 0, db = 0, dg = 0, aniso;
1227 double dec_m, dip_m, tau1, mu1, nu1, dt = 0, dm = 0, dn = 0, tau = 0, mu = 0, nu = 0;
1228 double wesn_new[4], out_igrf[7];
1229
1230 struct GRDREDPOL_CTRL *Ctrl = NULL;
1231 struct GMT_GRID *Gin = NULL, *Gout = NULL, *Gdip = NULL, *Gdec = NULL, *Gfilt = NULL;
1232 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
1233 struct GMT_OPTION *options = NULL;
1234 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1235
1236 /*----------------------- Standard module initialization and parsing ----------------------*/
1237
1238 if (API == NULL) return (GMT_NOT_A_SESSION);
1239 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
1240 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
1241
1242 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
1243
1244 /* Parse the command-line arguments */
1245
1246 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 */
1247 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
1248 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
1249 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
1250
1251 /*--------------------------- This is the grdredpol main code --------------------------*/
1252
1253 if ((Gin = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) /* Get header only */
1254 Return (API->error);
1255
1256 if (Ctrl->F.compute_n) {
1257 aniso = Gin->header->inc[GMT_X] / Gin->header->inc[GMT_Y] * cos(Gin->header->wesn[YHI]*D2R);
1258 Ctrl->F.ncoef_col = (int) ((double)Ctrl->F.ncoef_row / aniso);
1259 if (Ctrl->F.ncoef_col %2 != 1) Ctrl->F.ncoef_col++;
1260 }
1261
1262 m21 = (Ctrl->F.ncoef_row+1) / 2; n21 = (Ctrl->F.ncoef_col+1) / 2;
1263 GMT->current.io.pad[XLO] = GMT->current.io.pad[XHI] = n21-1;
1264 GMT->current.io.pad[YLO] = GMT->current.io.pad[YHI] = m21-1;
1265
1266 if (!GMT->common.R.active[RSET])
1267 gmt_M_memcpy (wesn_new, Gin->header->wesn, 4, double);
1268 else
1269 gmt_M_memcpy (wesn_new, GMT->common.R.wesn, 4, double);
1270
1271 one_or_zero = (Gin->header->registration == GMT_GRID_PIXEL_REG) ? 0 : 1;
1272 nx_new = urint ((wesn_new[XHI] - wesn_new[XLO]) / Gin->header->inc[GMT_X]) + one_or_zero;
1273 ny_new = urint ((wesn_new[YHI] - wesn_new[YLO]) / Gin->header->inc[GMT_Y]) + one_or_zero;
1274
1275 Ctrl->S.n_columns = nx_new; Ctrl->S.n_rows = ny_new;
1276
1277 gmt_grd_init (GMT, Gin->header, options, true);
1278
1279 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn_new, Ctrl->In.file, Gin) == NULL) { /* Get subset */
1280 Return (API->error);
1281 }
1282
1283 if (GMT->common.R.active[RSET]) {
1284 if (wesn_new[XLO] < Gin->header->wesn[XLO] || wesn_new[XHI] > Gin->header->wesn[XHI]) {
1285 GMT_Report (API, GMT_MSG_ERROR, " Selected region exceeds the X-boundaries of the grid file!\n");
1286 return (GMT_RUNTIME_ERROR);
1287 }
1288 else if (wesn_new[YLO] < Gin->header->wesn[YLO] || wesn_new[YHI] > Gin->header->wesn[YHI]) {
1289 GMT_Report (API, GMT_MSG_ERROR, " Selected region exceeds the Y-boundaries of the grid file!\n");
1290 return (GMT_RUNTIME_ERROR);
1291 }
1292 }
1293
1294 if (!Ctrl->M.pad_zero) /* That is, if we want edges reflected or replicated */
1295 grdgravmag3d_mirror_edges (Gin->data, Ctrl->F.ncoef_col-1, m21-1, n21-1, Ctrl);
1296
1297 /* Section to deal with possible external grids with dip and dec for interpolation */
1298
1299 if (Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd) {
1300 if ((Gdip = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->E.dipfile, NULL)) == NULL) /* Get header only */
1301 Return (API->error);
1302
1303 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn_new, Ctrl->E.dipfile, Gdip) == NULL)
1304 Return (API->error);
1305 }
1306 if (Ctrl->E.dip_dec_grd) {
1307 if ((Gdec = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->E.decfile, NULL)) == NULL)
1308 Return (API->error);
1309
1310 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, wesn_new, Ctrl->E.decfile, Gdec) == NULL)
1311 Return (API->error);
1312 }
1313
1314 n_coef = Ctrl->F.ncoef_row * Ctrl->F.ncoef_col;
1315 cosphi = gmt_M_memory (GMT, NULL, n_coef, double);
1316 sinphi = gmt_M_memory (GMT, NULL, n_coef, double);
1317 cospsi = gmt_M_memory (GMT, NULL, n_coef, double);
1318 sinpsi = gmt_M_memory (GMT, NULL, n_coef, double);
1319 gxr = gmt_M_memory (GMT, NULL, n_coef, double);
1320 gxi = gmt_M_memory (GMT, NULL, n_coef, double);
1321 gxar = gmt_M_memory (GMT, NULL, n_coef, double);
1322 gxai = gmt_M_memory (GMT, NULL, n_coef, double);
1323 gxbr = gmt_M_memory (GMT, NULL, n_coef, double);
1324 gxbi = gmt_M_memory (GMT, NULL, n_coef, double);
1325 gxgr = gmt_M_memory (GMT, NULL, n_coef, double);
1326 gxgi = gmt_M_memory (GMT, NULL, n_coef, double);
1327 fxr = gmt_M_memory (GMT, NULL, n_coef, double);
1328 fix = gmt_M_memory (GMT, NULL, n_coef, double);
1329 fxar = gmt_M_memory (GMT, NULL, n_coef, double);
1330 fxbr = gmt_M_memory (GMT, NULL, n_coef, double);
1331 fxgr = gmt_M_memory (GMT, NULL, n_coef, double);
1332 ftlon = gmt_M_memory (GMT, NULL, Gin->header->n_columns, double);
1333 ftlat = gmt_M_memory (GMT, NULL, Gin->header->n_rows, double);
1334
1335 if ((Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd)) {
1336 gxtr = gmt_M_memory (GMT, NULL, n_coef, double);
1337 gxti = gmt_M_memory (GMT, NULL, n_coef, double);
1338 gxmr = gmt_M_memory (GMT, NULL, n_coef, double);
1339 gxmi = gmt_M_memory (GMT, NULL, n_coef, double);
1340 gxnr = gmt_M_memory (GMT, NULL, n_coef, double);
1341 gxni = gmt_M_memory (GMT, NULL, n_coef, double);
1342 fxtr = gmt_M_memory (GMT, NULL, n_coef, double);
1343 fxmr = gmt_M_memory (GMT, NULL, n_coef, double);
1344 fxnr = gmt_M_memory (GMT, NULL, n_coef, double);
1345 }
1346
1347 /* Generate vectors of lon & lats */
1348 for (col = 0; col < Gin->header->n_columns; col++) ftlon[col] = gmt_M_grd_col_to_x (GMT, col, Gin->header);
1349 for (row = 0; row < Gin->header->n_rows; row++) ftlat[row] = gmt_M_grd_row_to_y (GMT, row, Gin->header);
1350
1351 n_jlon = urint ((Gin->header->wesn[XHI] - Gin->header->wesn[XLO]) / Ctrl->W.wid) + 1;
1352 n_jlat = urint ((Gin->header->wesn[YHI] - Gin->header->wesn[YLO]) / Ctrl->W.wid) + 1;
1353
1354 if (Ctrl->C.const_f) {
1355 alfa = -cos(Ctrl->C.dip) * cos(Ctrl->C.dec);
1356 beta = cos(Ctrl->C.dip) * sin(Ctrl->C.dec);
1357 gama = -sin(Ctrl->C.dip);
1358 }
1359
1360 fi = TWO_PI / Ctrl->F.ncoef_row;
1361 psi = TWO_PI / Ctrl->F.ncoef_col;
1362
1363 if ((Gout = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn_new, Gin->header->inc, \
1364 Gin->header->registration, GMT_NOTSET, NULL)) == NULL) Return (API->error);
1365
1366 if (Ctrl->Z.active) { /* Create one grid to hold the filter coefficients */
1367 double wesn[4], inc[2];
1368 wesn[XLO] = 1; wesn[XHI] = (double)Ctrl->F.ncoef_col;
1369 wesn[YLO] = 1; wesn[YHI] = (double)Ctrl->F.ncoef_row;
1370 inc[GMT_X] = inc[GMT_Y] = 1;
1371
1372 if ((Gfilt = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, wesn, inc,
1373 GMT_GRID_PIXEL_REG, 0, NULL)) == NULL) Return (API->error);
1374 strcpy (Gfilt->header->title, "Reduction To the Pole filter");
1375 strcpy (Gfilt->header->x_units, "radians");
1376 strcpy (Gfilt->header->y_units, "radians");
1377 }
1378
1379 for (l = 0; l < n_jlat; l++) { /* Main loop over the moving windows */
1380 for (k = 0; k < n_jlon; k++) {
1381 sloni = Gin->header->wesn[XLO] + k * Ctrl->W.wid;
1382 slati = Gin->header->wesn[YHI] - (l+1) * Ctrl->W.wid;
1383 slonf = sloni + Ctrl->W.wid;
1384 slatf = slati + Ctrl->W.wid;
1385 slonm = (sloni + slonf) / 2;
1386 slatm = (slati + slatf) / 2;
1387 if (Ctrl->F.compute_n) {
1388 aniso = Gin->header->inc[GMT_X] / Gin->header->inc[GMT_Y] * cos(slatm*D2R);
1389 Ctrl->F.ncoef_col = (int) ((double)Ctrl->F.ncoef_row / aniso);
1390 if (Ctrl->F.ncoef_col %2 != 1) Ctrl->F.ncoef_row += 1;
1391 psi = TWO_PI / Ctrl->F.ncoef_row;
1392 n21 = (Ctrl->F.ncoef_row+1) / 2;
1393 }
1394 /* Compute dec and dip at the central point of the moving window */
1395 grdgravmag3d_igrf10syn(GMT, 0, Ctrl->T.year, 1, 0, slonm, slatm, out_igrf);
1396 if (!Ctrl->C.const_f) {
1397 Ctrl->C.dec = out_igrf[5] * D2R;
1398 Ctrl->C.dip = out_igrf[6] * D2R;
1399 /* Compute the direction cosines */
1400 alfa = -cos(Ctrl->C.dip) * cos(Ctrl->C.dec);
1401 beta = cos(Ctrl->C.dip) * sin(Ctrl->C.dec);
1402 gama = -sin(Ctrl->C.dip);
1403 }
1404 if ((Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd)) { /* */
1405 if (Ctrl->E.dip_grd_only) { /* Use mag DEC = 0 */
1406 dip_m = gmt_bcr_get_z(GMT, Gdip, slonm, slatm) * D2R;
1407 tau = -cos(dip_m);
1408 mu = 0;
1409 nu = -sin(dip_m);
1410 }
1411 else { /* Get central window mag DEC & DIP from grids */
1412 dip_m = gmt_bcr_get_z(GMT, Gdip, slonm, slatm) * D2R;
1413 dec_m = gmt_bcr_get_z(GMT, Gdec, slonm, slatm) * D2R;
1414 tau = -cos(dip_m) * cos(dec_m);
1415 mu = cos(dip_m) * sin(dec_m);
1416 nu = -sin(dip_m);
1417 }
1418 }
1419 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION))
1420 GMT_Report (API, GMT_MSG_INFORMATION, "Dec %5.1f Dip %5.1f Bin_lon %6.1f Bin_lat %5.1f\r",
1421 Ctrl->C.dec/D2R, Ctrl->C.dip/D2R, slonm, slatm);
1422
1423 /* Compute the filter coefficients in the frequency domain */
1424 for (i = 0; i < Ctrl->F.ncoef_row; i++) {
1425 for (j = n21-1; j < Ctrl->F.ncoef_col; j++) {
1426 r = (double)(i - m21 + 1); s = (double)(j - n21 + 1);
1427 if (r == 0 && s == 0) continue;
1428 u = r * fi; v = s * psi;
1429
1430 if (!(Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd)) {
1431 grdgravmag3d_rtp_filt_colinear(i,j,n21,gxr,gxi,gxar,gxai,gxbr,gxbi,gxgr,gxgi,
1432 u,v,alfa,beta,gama, Ctrl);
1433 }
1434 else { /* Case more general. alfa, beta, gama, tau, mu e nu */
1435 grdgravmag3d_rtp_filt_not_colinear(i,j,n21,gxr,gxi,gxar,gxai,gxbr,gxbi,gxgr,gxgi,gxtr,gxti,
1436 gxmr,gxmi,gxnr,gxni,u,v,alfa,beta,gama,tau,mu,nu, Ctrl);
1437 }
1438 }
1439 }
1440 gxr[ij_mn(Ctrl,m21-1,0)] = 1; gxi[ij_mn(Ctrl,m21-1,0)] = 0;
1441 gxar[ij_mn(Ctrl,m21-1,0)] = 0; gxai[ij_mn(Ctrl,m21-1,0)] = 0;
1442 gxbr[ij_mn(Ctrl,m21-1,0)] = 0; gxbi[ij_mn(Ctrl,m21-1,0)] = 0;
1443 gxgr[ij_mn(Ctrl,m21-1,0)] = 0; gxgi[ij_mn(Ctrl,m21-1,0)] = 0;
1444
1445 /* Compute iFT of the filter */
1446 grdgravmag3d_tfpoeq(fxr, Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxr,gxi, cosphi,sinphi,cospsi,sinpsi);
1447 grdgravmag3d_tfpoeq(fxar,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxar,gxai, cosphi,sinphi,cospsi,sinpsi);
1448 grdgravmag3d_tfpoeq(fxbr,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxbr,gxbi, cosphi,sinphi,cospsi,sinpsi);
1449 grdgravmag3d_tfpoeq(fxgr,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxgr,gxgi, cosphi,sinphi,cospsi,sinpsi);
1450
1451 if ((Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd)) {
1452 gxtr[ij_mn(Ctrl,m21-1,0)] = 0; gxti[ij_mn(Ctrl,m21-1,0)] = 0;
1453 gxmr[ij_mn(Ctrl,m21-1,0)] = 0; gxmi[ij_mn(Ctrl,m21-1,0)] = 0;
1454 gxnr[ij_mn(Ctrl,m21-1,0)] = 0; gxni[ij_mn(Ctrl,m21-1,0)] = 0;
1455 grdgravmag3d_tfpoeq(fxtr,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxtr,gxti, cosphi,sinphi,cospsi,sinpsi);
1456 grdgravmag3d_tfpoeq(fxmr,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxmr,gxmi, cosphi,sinphi,cospsi,sinpsi);
1457 grdgravmag3d_tfpoeq(fxnr,Ctrl->F.ncoef_row, Ctrl->F.ncoef_col,gxnr,gxni, cosphi,sinphi,cospsi,sinpsi);
1458 }
1459
1460 /* Convolve filter with input data that is inside current window (plus what filter width imposes) */
1461 gmt_M_row_loop (GMT, Gout,row) {
1462 if (ftlat[row] < slati || ftlat[row] > slatf) continue; /* Current point outside WOI */
1463 gmt_M_col_loop (GMT, Gout,row,col,ij) {
1464 if (ftlon[col] < sloni || ftlon[col] > slonf) continue; /* Current point outside WOI */
1465 /* Compute dec and dip at current point */
1466 if (!Ctrl->C.const_f) { /* It means we need to get F (& M) vector parameters */
1467 grdgravmag3d_igrf10syn(GMT, 0, Ctrl->T.year, 1, 0, ftlon[col], ftlat[row], out_igrf);
1468 Ctrl->C.dec = out_igrf[5] * D2R;
1469 Ctrl->C.dip = out_igrf[6] * D2R;
1470 if (Ctrl->E.dip_grd_only) {
1471 dip_m = gmt_bcr_get_z(GMT, Gdip, ftlon[row], ftlat[row]) * D2R;
1472 tau1 = -cos(dip_m);
1473 mu1 = 0;
1474 nu1 = -sin(dip_m);
1475 dt = tau1 - tau; dm = mu1 - mu; dn = nu1 - nu;
1476 }
1477 else if (Ctrl->E.dip_dec_grd) {
1478 dec_m = gmt_bcr_get_z(GMT, Gdec, ftlon[col], ftlat[row]) * D2R;
1479 dip_m = gmt_bcr_get_z(GMT, Gdip, ftlon[col], ftlat[row]) * D2R;
1480 tau1 = -cos(dip_m) * cos(dec_m);
1481 mu1 = cos(dip_m) * sin(dec_m);
1482 nu1 = -sin(dip_m);
1483 dt = tau1 - tau; dm = mu1 - mu; dn = nu1 - nu;
1484 }
1485 /* Compute director cosinus */
1486 alfa1 = -cos(Ctrl->C.dip) * cos(Ctrl->C.dec);
1487 beta1 = cos(Ctrl->C.dip) * sin(Ctrl->C.dec);
1488 gama1 = -sin(Ctrl->C.dip);
1489 da = alfa1 - alfa; db = beta1 - beta; dg = gama1 - gama;
1490 }
1491 if (!Ctrl->N.active) /* Do not use the Taylor expansion (What's the interest?) */
1492 da = db = dg = dt = dm = dn = 0;
1493
1494 /* Rebuild the filter */
1495 if (!Ctrl->C.const_f) {
1496 if (Ctrl->C.use_igrf) {
1497 for (i2 = 0; i2 < n_coef; i2++)
1498 fix[i2] = fxr[i2] + da * fxar[i2] + db * fxbr[i2] + dg * fxgr[i2];
1499 }
1500 else {
1501 for (i2 = 0; i2 < n_coef; i2++)
1502 fix[i2] = fxr[i2] + da * fxar[i2] + db * fxbr[i2] + dg * fxgr[i2] +
1503 dt * fxtr[i2] + dm * fxmr[i2] + dn * fxnr[i2];
1504 }
1505 }
1506 else
1507 gmt_M_memcpy (fix, fxr, n_coef, double);
1508
1509 if (Ctrl->Z.active && !wrote_one && l == 0 && k == 0) {
1510 for (jj = i2 = 0; i2 < Ctrl->F.ncoef_row; i2++) /* Remember, filter is columnwise */
1511 for (j2 = 0; j2 < Ctrl->F.ncoef_col; j2++, jj++)
1512 Gfilt->data[jj] = (gmt_grdfloat)fix[ij_mn(Ctrl,i2,j2)];
1513
1514 wrote_one = true;
1515 }
1516
1517 tmp_d = 0;
1518 for (i2 = 0; i2 < Ctrl->F.ncoef_row; i2++) {
1519 i3 = row + i2;
1520 for (j2 = 0; j2 < Ctrl->F.ncoef_col; j2++)
1521 tmp_d += fix[ij_mn(Ctrl,i2,j2)] * Gin->data[ij0_data(Ctrl,i3,(col + j2))];
1522 }
1523
1524 Gout->data[ij] = (gmt_grdfloat)tmp_d;
1525 }
1526 }
1527 }
1528 }
1529
1530 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) GMT_Report (API, GMT_MSG_LONG_VERBOSE, "\n");
1531
1532 gmt_M_free (GMT, cosphi); gmt_M_free (GMT, sinphi);
1533 gmt_M_free (GMT, cospsi); gmt_M_free (GMT, sinpsi);
1534 gmt_M_free (GMT, gxr); gmt_M_free (GMT, gxi);
1535 gmt_M_free (GMT, ftlat); gmt_M_free (GMT, ftlon);
1536 gmt_M_free (GMT, fxr);
1537
1538 gmt_M_free (GMT, fxar); gmt_M_free (GMT, fxbr);
1539 gmt_M_free (GMT, fxgr); gmt_M_free (GMT, fix);
1540
1541 gmt_M_free (GMT, gxar); gmt_M_free (GMT, gxai);
1542 gmt_M_free (GMT, gxbr); gmt_M_free (GMT, gxbi);
1543 gmt_M_free (GMT, gxgr); gmt_M_free (GMT, gxgi);
1544 if ((Ctrl->E.dip_grd_only || Ctrl->E.dip_dec_grd)) {
1545 gmt_M_free (GMT, gxtr); gmt_M_free (GMT, gxti);
1546 gmt_M_free (GMT, gxmr); gmt_M_free (GMT, gxmi);
1547 gmt_M_free (GMT, gxnr); gmt_M_free (GMT, gxni);
1548 gmt_M_free (GMT, fxtr); gmt_M_free (GMT, fxmr);
1549 gmt_M_free (GMT, fxnr);
1550 }
1551
1552 strcpy (Gout->header->title, "Anomaly reduced to the pole");
1553 strcpy (Gout->header->z_units, "nT");
1554
1555 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Gout)) Return (API->error);
1556 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Gout) != GMT_NOERROR) {
1557 Return (API->error);
1558 }
1559 if (Ctrl->Z.active) {
1560 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Gfilt)) Return (API->error);
1561 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Z.file, Gfilt) != GMT_NOERROR) {
1562 Return (API->error);
1563 }
1564 }
1565
1566 Return (GMT_NOERROR);
1567 }
1568