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 * Author: Paul Wessel
19 * Date: 1-JAN-2010
20 * Version: 6 API
21 *
22 * Brief synopsis: greenspline grids data using Green's functions for a selected spline.
23 * The data may be Cartesian or geographical, and gridding can be done
24 * in a Cartesian 1-D, 2-D, 3-D space or on a sphere. The spline may be evaluated on
25 * a grid, on parts of a grid, or at specified arbitrary locations.
26 * Five classes of splines (for a total of 10 splines) are implemented:
27 * 1. Minimum curvature Cartesian spline [Sandwell, Geophys. Res. Lett, 1987]
28 * 2. Minimum curvature Cartesian spline in tension [Wessel & Bercovici, Math., Geol., 1998]
29 * 3. Regularized Cartesian spline in tension [Mitasova & Mitas, Math. Geol., 1993]
30 * 4. Minimum curvature spherical spline [Parker, "Geophysical Inverse Theory", 1994]
31 * 5. Minimum curvature spherical spline in tension [Wessel & Becker, Geophys. J. Int, 2008]
32 *
33 * Originally published as:
34 * "Wessel, P., 2009. A general-purpose Green's function-based interpolator,
35 * Computers & Geosciences, 35: 1247-1254".
36 *
37 * PW Update June 2013. The numerical implementation of the Green's function found by Wessel & Becker [2008]
38 * was unstable (it required the difference between k*P_v and log, and P_v is difficult to compute accurately).
39 * Bob Parker (Scripps) helped develop a series solution for canceling out the two singular behaviors,
40 * resulting in a more stable expression that converges reasonably rapidly. We now use this new series
41 * solution for -Sq, combined with a (new) cubic spline interpolation. This replaces the old -SQ machinery
42 * with linear interpolation which is now deprecated.
43 *
44 * PW Update July 2015. With help from Dong Ju Choi, San Diego Supercomputing Center, we have added Open MP
45 * support in greenspline and the matrix solvers in gmt_vector.c. Requires open MP support and use of -x.
46 */
47
48 #include "gmt_dev.h"
49
50 #define THIS_MODULE_CLASSIC_NAME "greenspline"
51 #define THIS_MODULE_MODERN_NAME "greenspline"
52 #define THIS_MODULE_LIB "core"
53 #define THIS_MODULE_PURPOSE "Interpolate using Green's functions for splines in 1-3 dimensions"
54 #define THIS_MODULE_KEYS "<D{,AD(=,ED),ND(,TG(,CD)=f,G?},GDN"
55 #define THIS_MODULE_NEEDS "R"
56 #define THIS_MODULE_OPTIONS "-:>Vbdefghioqrsw" GMT_OPT("FH") GMT_ADD_x_OPT
57
58 EXTERN_MSC int gmtlib_cspline (struct GMT_CTRL *GMT, double *x, double *y, uint64_t n, double *c);
59
60 /* Control structure for greenspline */
61
62 struct GREENSPLINE_CTRL {
63 unsigned int dimension;
64 struct GREENSPLINE_A { /* -A<gradientfile> */
65 bool active;
66 unsigned int mode; /* 0 = azimuths, 1 = directions, 2 = dx,dy components, 3 = dx, dy, dz components */
67 char *file;
68 } A ;
69 struct GREENSPLINE_C { /* -C[[n|r|v]<cutoff>[%]][+c][+f<file>][+i][+n] */
70 bool active;
71 bool dryrun; /* Only report eigenvalues */
72 unsigned int history;
73 unsigned int mode;
74 double value;
75 char *file;
76 } C;
77 struct GREENSPLINE_D { /* -D[+x<xname>][+yyname>][+z<zname>][+v<zname>][+s<scale>][+o<offset>][+n<invalid>][+t<title>][+r<remark>] */
78 bool active;
79 char *information;
80 } D;
81 struct GREENSPLINE_E { /* -E[<file>] */
82 bool active;
83 unsigned int mode;
84 char *file;
85 } E;
86 struct GREENSPLINE_G { /* -G<output_grdfile> */
87 bool active;
88 char *file;
89 } G;
90 struct GREENSPLINE_I { /* -Idx[/dy[/dz]] */
91 bool active;
92 double inc[3];
93 } I;
94 struct GREENSPLINE_L { /* -L */
95 bool active;
96 } L;
97 struct GREENSPLINE_M { /* -M<gfuncfile> */
98 bool active;
99 unsigned int mode; /* GMT_IN or GMT_OUT */
100 char *file;
101 } M;
102 struct GREENSPLINE_N { /* -N<outputnode_file> */
103 bool active;
104 char *file;
105 } N;
106 struct GREENSPLINE_Q { /* -Qdaz */
107 bool active;
108 double az;
109 double dir[3];
110 } Q;
111 struct GREENSPLINE_R3 { /* -Rxmin/xmax[/ymin/ymax[/zmin/zmaz]] | -Ggridfile */
112 bool active;
113 bool mode; /* true if settings came from a grid file */
114 unsigned int dimension; /* 1, 2, or 3 */
115 unsigned int offset; /* 0 or 1 */
116 double range[6]; /* Min/max for each dimension */
117 double inc[2]; /* xinc/yinc when -Rgridfile was given*/
118 } R3;
119 struct GREENSPLINE_S { /* -S<mode>[<tension][+<mod>[args]] */
120 bool active;
121 unsigned int mode;
122 double value[4];
123 double rval[2];
124 char *arg;
125 } S;
126 struct GREENSPLINE_T { /* -T<mask_grdfile> */
127 bool active;
128 char *file;
129 } T;
130 struct GREENSPLINE_W { /* -W[w] */
131 bool active;
132 unsigned int mode; /* 0 = got sigmas, 1 = got weights */
133 } W;
134 struct GREENSPLINE_Z { /* -Z<distflag> */
135 bool active;
136 int mode; /* Can be negative */
137 } Z;
138 struct GREENSPLINE_DEBUG { /* -0 undocumented debugging option */
139 bool active;
140 } debug;
141 };
142
143 enum Greenspline_modes { /* Various integer mode flags */
144 SANDWELL_1987_1D = 0,
145 SANDWELL_1987_2D = 1,
146 SANDWELL_1987_3D = 2,
147 WESSEL_BERCOVICI_1998_1D = 3,
148 WESSEL_BERCOVICI_1998_2D = 4,
149 WESSEL_BERCOVICI_1998_3D = 5,
150 MITASOVA_MITAS_1993_2D = 6,
151 MITASOVA_MITAS_1993_3D = 7,
152 PARKER_1994 = 8,
153 WESSEL_BECKER_2008 = 9,
154 LINEAR_1D = 10,
155 LINEAR_2D = 11,
156 N_METHODS = 12,
157 N_PARAMS = 11,
158 GREENSPLINE_TREND = 1, /* Remove/Restore linear trend */
159 GREENSPLINE_NORM = 2, /* Normalize residual data to 0-1 range */
160 SQ_N_NODES = 10001, /* Default number of nodes in the precalculated -Sq spline */
161 GSP_GOT_SIG = 0,
162 GSP_GOT_WEIGHTS = 1
163 };
164
165 #ifndef M_LOG_2
166 #define M_LOG_2 0.69314718055994530942
167 #endif
168 #ifndef M_GAMMA
169 #define M_GAMMA 0.577215664901532860606512
170 #endif
171 #ifndef M_SQRT_PI
172 #define M_SQRT_PI 1.772453850905516027298167483341
173 #endif
174 #ifndef M_INV_SQRT_PI
175 #define M_INV_SQRT_PI (1.0 / M_SQRT_PI)
176 #endif
177
178 #define SQ_TRUNC_ERROR 1.0e-6 /* Max truncation error in Parker's simplified sum for WB'08 */
179
180 enum Greenspline_index { /* Indices for coeff array for normalization */
181 GSP_MEAN_X = 0,
182 GSP_MEAN_Y = 1,
183 GSP_MEAN_Z = 2,
184 GSP_SLP_X = 3,
185 GSP_SLP_Y = 4,
186 GSP_RANGE = 5,
187 GSP_LENGTH = 6};
188
189 struct GREENSPLINE_LOOKUP { /* Used to spline interpolation of precalculated function */
190 uint64_t n; /* Number of values in the spline setup */
191 double *y; /* Function values */
192 double *c; /* spline coefficients */
193 double *A, *B, *C; /* power/ratios of order l terms */
194 };
195
196 #ifdef DEBUG
197 static bool TEST = false; /* Global variable used for undocumented testing [under -DDEBUG only; see -+ hidden option] */
198 #endif
199
New_Ctrl(struct GMT_CTRL * GMT)200 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
201 struct GREENSPLINE_CTRL *C;
202
203 C = gmt_M_memory (GMT, NULL, 1, struct GREENSPLINE_CTRL);
204
205 /* Initialize values whose defaults are not 0/false/NULL */
206 C->C.mode = GMT_SVD_EIGEN_RATIO_CUTOFF;
207 C->S.mode = SANDWELL_1987_2D;
208 C->S.rval[0] = -1.0; C->S.rval[1] = 1.0;
209 C->S.value[3] = (double)SQ_N_NODES; /* Default number of spline nodes */
210 C->S.value[2] = SQ_TRUNC_ERROR; /* Default truncation error for Legendre sum in -Sq */
211 return (C);
212 }
213
Free_Ctrl(struct GMT_CTRL * GMT,struct GREENSPLINE_CTRL * C)214 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *C) { /* Deallocate control structure */
215 if (!C) return;
216 gmt_M_str_free (C->A.file);
217 gmt_M_str_free (C->C.file);
218 gmt_M_str_free (C->G.file);
219 gmt_M_str_free (C->N.file);
220 gmt_M_str_free (C->T.file);
221 gmt_M_str_free (C->S.arg);
222 gmt_M_free (GMT, C);
223 }
224
usage(struct GMTAPI_CTRL * API,int level)225 static int usage (struct GMTAPI_CTRL *API, int level) {
226 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
227 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
228 GMT_Usage (API, 0, "usage: %s [<table>] -G<outfile> [-A<gradientfile>+f<format>] [-C[[n|r|v]<val>[%%]][+c][+f<file>][+i][+n]] "
229 "[-D<information>] [-E[<misfitfile>]] [-I<dx>[/<dy>[/<dz>]]] [-L] [-N<nodefile>] [-Q<az>] "
230 "[-R<xmin>/<xmax[/<ymin>/<ymax>[/<zmin>/<zmax>]]] [-Sc|l|t|r|p|q[<pars>]] [-T<maskgrid>] "
231 "[%s] [-W[w]] [-Z<mode>] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s [%s]%s[%s] [%s]\n",
232 name, GMT_V_OPT,GMT_bi_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT,
233 GMT_o_OPT, GMT_q_OPT, GMT_r_OPT, GMT_s_OPT, GMT_w_OPT, GMT_x_OPT, GMT_colon_OPT, GMT_PAR_OPT);
234
235 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
236
237 GMT_Usage (API, -1, "Choose one of three ways to specify where to evaluate the spline:");
238 GMT_Usage (API, 2, "%s Specify an equidistant 1-, 2-, or 3-D domain with options -R, -I [and optionally -r].", GMT_LINE_BULLET);
239 GMT_Usage (API, 2, "%s Specify a set of 1- ,2-, or 3-D output locations via the -N option.", GMT_LINE_BULLET);
240 GMT_Usage (API, 2, "%s Supply a mask grid file via -T whose values are NaN or 0. The spline will then "
241 "only be evaluated at the nodes originally set to zero (2-D only).", GMT_LINE_BULLET);
242 GMT_Message (API, GMT_TIME_NONE, "\n REQUIRED ARGUMENTS:\n");
243 GMT_Option (API, "<");
244 GMT_Usage (API, 1, "\n-G<outfile>");
245 GMT_Usage (API, -2, "Output data. Append name of output file.");
246
247 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
248
249 GMT_Usage (API, 1, "\n-A<gradientfile>+f<format>");
250 GMT_Usage (API, -2, "ASCII file with surface gradients V to use in the modeling. Specify format:");
251 GMT_Usage (API, 3, "0: For 1-D: x, slope.");
252 GMT_Usage (API, 3, "1: X, Vmagnitude, Vazimuth(s).");
253 GMT_Usage (API, 3, "2: X, Vazimuth(s), Vmagnitude.");
254 GMT_Usage (API, 3, "3: X, Vmagnitude, Vangle(s).");
255 GMT_Usage (API, 3, "4: X, Vcomponents.");
256 GMT_Usage (API, 3, "5: X, Vunit-vector, Vmagnitude.");
257 GMT_Usage (API, -2, "Here, X = (x, y[, z]) is the position vector, V = (Vx, Vy[, Vz]) is the gradient vector.");
258 GMT_Usage (API, 1, "\n-C[[n|r|v]<val>[%%]][+c][+f<file>][+i][+n]");
259 GMT_Usage (API, -2, "Solve by SVD and control how many eigenvalues to use. Optionally append a directive and value:");
260 GMT_Usage (API, 3, "n: Only use the largest <val> eigenvalues [all].");
261 GMT_Usage (API, 3, "r: Eliminate eigenvalues whose ratio to largest eigenvalue is less than <val> [0].");
262 GMT_Usage (API, 3, "v: Include eigenvalues needed to reach a variance explained fraction of <val> [1].");
263 GMT_Usage (API, -2, "Note: For r|v you may append %% to give <val> as the percentage of total instead. Various optional modifiers are available:");
264 GMT_Usage (API, 3, "+c Valid for 2-D gridding only and will create a series of intermediate "
265 "grids for each eigenvalue holding the cumulative result. Requires -G with a valid filename "
266 "and extension and we will insert _cum_### before the extension.");
267 GMT_Usage (API, 3, "+f Save the eigenvalues to <filename>.");
268 GMT_Usage (API, 3, "+i As +c but save incremental results, inserting _inc_### before the extension.");
269 GMT_Usage (API, 3, "+n Stop execution after reporting the eigenvalues - no solution is computed.");
270 GMT_Usage (API, -2, "Note: Without -C we use Gauss-Jordan elimination to solve the linear system.");
271 gmt_grdcube_info_syntax (API->GMT, 'D');
272 GMT_Usage (API, 1, "\n-E[<misfitfile>]");
273 GMT_Usage (API, -2, "Evaluate solution at input locations and report misfit statistics. "
274 "Append a filename to save all data with two extra columns for model and misfit. "
275 "If -C+i|c are used then we instead report the history of model variance and rms misfit.");
276 GMT_Usage (API, 1, "\nI<dx>[/<dy>[/<dz>]]");
277 GMT_Usage (API, -2, "Specify a regular set of output locations. Give equidistant increment for each dimension. "
278 "Requires -R for specifying the output domain.");
279 GMT_Usage (API, 1, "\n-L Leave trend alone. Do not remove least squares plane from data before spline fit. "
280 "Only applies to -D0-2 [Default removes linear trend, fits residuals, and restores trend].");
281 GMT_Usage (API, 1, "\n-N<nodefile>");
282 GMT_Usage (API, -2, "ASCII file with desired output locations. "
283 "The resulting ASCII coordinates and interpolation are written to file given in -G "
284 "or standard output if no file specified (see -bo for binary output).");
285 GMT_Usage (API, 1, "\n-Q<az>");
286 GMT_Usage (API, -2, "Calculate the directional derivative in the <az> direction and return it instead of surface elevation.");
287 GMT_Option (API, "R");
288 if (gmt_M_showusage (API)) {
289 GMT_Usage (API, -2, "-R Specify a regular set of output locations. Give min and max coordinates for each dimension. "
290 "Requires -I for specifying equidistant increments. For 2D-gridding a gridfile may be given; "
291 "this then also sets -I (and perhaps -r); use those options to override the grid settings.");
292 }
293 GMT_Usage (API, 1, "\n-Sc|l|t|r|p|q[<pars>]");
294 GMT_Usage (API, -2, "Specify which spline to use; except for c|p, append normalized <tension> between 0 and 1:");
295 GMT_Usage (API, 3, "c: Minimum curvature spline (Sandwell, 1987) [Default].");
296 GMT_Usage (API, 3, "l: Linear (1-D) or bilinear (2-D) spline.");
297 GMT_Usage (API, 3, "t: Cartesian spline in tension (Wessel & Bercovici, 1998). Append <tension> and "
298 "optionally append /<scale> for length-scale [Default is the given output spacing].");
299 GMT_Usage (API, 3, "r: Regularized spline in tension (Mitasova & Mitas, 1993). Append <tension> and "
300 "optionally append /<scale> for length-scale [Default is given output spacing].");
301 GMT_Usage (API, 3, "p: Spherical surface spline (Parker, 1994); automatically sets -D4.");
302 GMT_Usage (API, 3, "q: Spherical surface spline in tension (Wessel & Becker, 2008); automatically sets -D4. Append <tension>. "
303 "Optionally, append +e<error> to change maximum error in series truncation [%g] and "
304 "+n<n> to change the (odd) number of precalculated nodes for spline interpolation [%d].", SQ_TRUNC_ERROR, SQ_N_NODES);
305 GMT_Usage (API, 1, "\n-T<maskgrid>");
306 GMT_Usage (API, -2, "Mask grid file whose values are NaN or 0; its header implicitly sets -R, -I (and -r).");
307 GMT_Usage (API, 1, "\n-W[w]");
308 GMT_Usage (API, -2, "Expect one extra input column with data errors sigma_i. "
309 "Append w to indicate this column carries weights instead "
310 "[Default makes weights via w_i = 1/sigma_i^2] for the least squares solution.");
311 GMT_Usage (API, -2, "Note: weights only have an effect if -C is used.");
312 GMT_Usage (API, 1, "\n-Z<mode>");
313 GMT_Usage (API, -2, "Distance <mode> determines how we calculate distances between (x,y) points. "
314 "Mode 0 applies to Cartesian 1-D spline interpolation:");
315 GMT_Usage (API, 3, "0: x in user units, Cartesian distances.");
316 GMT_Usage (API, -2, "Modes 1-3 apply to Cartesian 2-D surface spline interpolation:");
317 GMT_Usage (API, 3, "1: x,y in user units, Cartesian distances.");
318 GMT_Usage (API, 3, "2: x,y in degrees, flat Earth distances in meters.");
319 GMT_Usage (API, 3, "3: x,y in degrees, spherical distances in meters.");
320 GMT_Usage (API, -2, "Mode 4 applies to 2-D spherical surface spline interpolation:");
321 GMT_Usage (API, 3, "4: x,y in degrees, use cosine of spherical distances in degrees.");
322 GMT_Usage (API, -2, "Mode 5 applies to Cartesian 3-D volume interpolation:");
323 GMT_Usage (API, 3, "5: x,y,z in user units, Cartesian distances.");
324 GMT_Usage (API, -2, "Note: For modes 3-4, use PROJ_ELLIPSOID to select geodesic or great circle arcs.");
325 GMT_Option (API, "V,bi");
326 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Default is 2-5 input columns depending on dimensionality (see -Z) and weights (see -W).");
327 GMT_Option (API, "d,e,f,g,h,i,o,q,r,s,w,x,:,.");
328
329 return (GMT_MODULE_USAGE);
330 }
331
greenspline_pre_parser(struct GMT_CTRL * GMT,struct GMT_OPTION * options)332 GMT_LOCAL unsigned int greenspline_pre_parser (struct GMT_CTRL *GMT, struct GMT_OPTION *options) {
333 /* Help GMT_Parse know if -R is geographic based on -Z mode and return dimension */
334 unsigned int dim = 0;
335 struct GMT_OPTION *opt = NULL;
336 for (opt = options; opt; opt = opt->next) { /* Look for -Z only */
337 if (opt->option != 'Z' || opt->arg[0] == '\0') continue;
338 switch (opt->arg[0]) {
339 case '0': dim = 1; break;
340 case '1': dim = 2; break;
341 case '2': case '3': case '4':
342 dim = 2;
343 gmt_set_geographic (GMT, GMT_IN);
344 gmt_set_geographic (GMT, GMT_OUT);
345 break;
346 case '5': dim = 3; break;
347 default: dim = 0; break;
348 }
349 }
350 return (dim);
351 }
352
parse(struct GMT_CTRL * GMT,struct GREENSPLINE_CTRL * Ctrl,struct GMT_OPTION * options)353 static int parse (struct GMT_CTRL *GMT, struct GREENSPLINE_CTRL *Ctrl, struct GMT_OPTION *options) {
354 /* This parses the options provided to greenspline and sets parameters in CTRL.
355 * Any GMT common options will override values set previously by other commands.
356 * It also replaces any file names specified as input or output with the data ID
357 * returned when registering these sources/destinations with the API.
358 */
359
360 int n_items;
361 unsigned int n_errors = 0, dimension, k, pos = 0;
362 char txt[6][GMT_LEN64], p[GMT_BUFSIZ] = {""}, *c = NULL, *i = NULL, *r = NULL;
363 struct GMT_OPTION *opt = NULL;
364 struct GMTAPI_CTRL *API = GMT->parent;
365
366 for (opt = options; opt; opt = opt->next) {
367 switch (opt->option) {
368
369 case '<': /* Skip input files */
370 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
371 break;
372
373 case 'R': /* Normally processed internally but must be handled separately since it can take 1,2,3 dimensions */
374 GMT->common.R.active[RSET] = true;
375 Ctrl->R3.dimension = 1; /* At least */
376 if (GMT->current.setting.run_mode == GMT_MODERN && Ctrl->dimension == 2) { /* Watch for multi-item -R string created by gmt_init_module that would not have been parsed */
377 /* This is needed because we are not using gmt_parse_R_option in greenspline */
378 if ((r = strstr (opt->arg, "+G"))) { /* Got grid registration implicitly via history */
379 switch (r[2]) {
380 case 'G': GMT->common.R.registration = GMT_GRID_NODE_REG; break;
381 default: GMT->common.R.registration = GMT_GRID_PIXEL_REG; break;
382 }
383 r[0] = '\0'; /* Chop off this modifier */
384 GMT->common.R.active[GSET] = true;
385 }
386 if ((i = strstr (opt->arg, "+I"))) { /* Got grid increments implicitly via history */
387 Ctrl->I.active = true;
388 k = gmt_getincn (GMT, &i[2], Ctrl->I.inc, 2);
389 i[0] = '\0'; /* Chop off this modifier */
390 GMT->common.R.active[ISET] = true;
391 }
392 }
393 if (opt->arg[0] == 'g' && opt->arg[1] == '\0') { /* Got -Rg */
394 Ctrl->R3.range[0] = 0.0; Ctrl->R3.range[1] = 360.0; Ctrl->R3.range[2] = -90.0; Ctrl->R3.range[3] = 90.0;
395 Ctrl->R3.dimension = 2;
396 break;
397 }
398 if (opt->arg[0] == 'g' && opt->arg[1] == '/') { /* Got -Rg/zmin/zmax */
399 Ctrl->R3.range[0] = 0.0; Ctrl->R3.range[1] = 360.0; Ctrl->R3.range[2] = -90.0; Ctrl->R3.range[3] = 90.0;
400 n_items = sscanf (&opt->arg[2], "%[^/]/%s", txt[4], txt[5]);
401 if (n_items != 2) {
402 GMT_Report (API, GMT_MSG_ERROR, "Option -Rg/z0/z1: Append the z-range\n");
403 n_errors++;
404 }
405 Ctrl->R3.dimension = 3;
406 break;
407 }
408 if (opt->arg[0] == 'd' && opt->arg[1] == '\0') { /* Got -Rd */
409 Ctrl->R3.range[0] = -180.0; Ctrl->R3.range[1] = 180.0; Ctrl->R3.range[2] = -90.0; Ctrl->R3.range[3] = 90.0;
410 Ctrl->R3.dimension = 2;
411 break;
412 }
413 if (opt->arg[0] == 'd' && opt->arg[1] == '/') { /* Got -Rd/zmin/zmax */
414 Ctrl->R3.range[0] = -180.0; Ctrl->R3.range[1] = 180.0; Ctrl->R3.range[2] = -90.0; Ctrl->R3.range[3] = 90.0;
415 n_items = sscanf (&opt->arg[2], "%[^/]/%s", txt[4], txt[5]);
416 if (n_items != 2) {
417 GMT_Report (API, GMT_MSG_ERROR, "Option -Rd/z0/z1: Append the z-range\n");
418 n_errors++;
419 }
420 Ctrl->R3.dimension = 3;
421 break;
422 }
423 if (!gmt_access (GMT, opt->arg, R_OK)) { /* Gave a readable file, presumably a grid */
424 struct GMT_GRID *G = NULL;
425 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, opt->arg, NULL)) == NULL) { /* Get header only */
426 return (API->error);
427 }
428 Ctrl->R3.range[0] = G->header->wesn[XLO]; Ctrl->R3.range[1] = G->header->wesn[XHI];
429 Ctrl->R3.range[2] = G->header->wesn[YLO]; Ctrl->R3.range[3] = G->header->wesn[YHI];
430 Ctrl->R3.inc[GMT_X] = G->header->inc[GMT_X]; Ctrl->R3.inc[GMT_Y] = G->header->inc[GMT_Y];
431 Ctrl->R3.offset = G->header->registration;
432 Ctrl->R3.dimension = 2;
433 Ctrl->R3.mode = true;
434 if (GMT_Destroy_Data (API, &G) != GMT_NOERROR) {
435 return (API->error);
436 }
437 break;
438 }
439 /* Only get here if the above cases did not trip */
440 n_items = sscanf (opt->arg, "%[^/]/%[^/]/%[^/]/%[^/]/%[^/]/%s", txt[0], txt[1], txt[2], txt[3], txt[4], txt[5]);
441 if (!(n_items == 2 || n_items == 4 || n_items == 6)) {
442 GMT_Report (API, GMT_MSG_ERROR, "Option -R: Give 2, 4, or 6 coordinates\n");
443 n_errors++;
444 }
445 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt[0], gmt_M_type (GMT, GMT_IN, GMT_X), true, &Ctrl->R3.range[0]), txt[0]);
446 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt[1], gmt_M_type (GMT, GMT_IN, GMT_X), true, &Ctrl->R3.range[1]), txt[1]);
447 if (n_items > 2) {
448 Ctrl->R3.dimension = 2;
449 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt[2], gmt_M_type (GMT, GMT_IN, GMT_Y), true, &Ctrl->R3.range[2]), txt[2]);
450 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt[3], gmt_M_type (GMT, GMT_IN, GMT_Y), true, &Ctrl->R3.range[3]), txt[3]);
451 }
452 if (n_items == 6) {
453 Ctrl->R3.dimension = 3;
454 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Z), gmt_scanf_arg (GMT, txt[4], gmt_M_type (GMT, GMT_IN, GMT_Z), false, &Ctrl->R3.range[4]), txt[4]);
455 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Z), gmt_scanf_arg (GMT, txt[5], gmt_M_type (GMT, GMT_IN, GMT_Z), false, &Ctrl->R3.range[5]), txt[5]);
456 }
457 if (Ctrl->R3.dimension > 1) gmt_M_memcpy (GMT->common.R.wesn, Ctrl->R3.range, 4, double);
458 if (r) r[0] = '+'; /* Restore modifier */
459 if (i) i[0] = '+'; /* Restore modifier */
460 break;
461
462 /* Processes program-specific parameters */
463
464 case 'A': /* Gradient data: -A<gradientfile>+f<format> */
465 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
466 Ctrl->A.active = true;
467 if (strchr (opt->arg, ',')) { /* Old syntax: Specified a particular format with -A<mode>,<file> */
468 if (gmt_M_compat_check (API->GMT, 5)) {
469 GMT_Report (API, GMT_MSG_COMPAT, "Option -A<format>,<gradientfile> is deprecated; use -A<gradientfile>+f<format> instead\n");
470 Ctrl->A.mode = (int)(opt->arg[0] - '0');
471 Ctrl->A.file = strdup (&opt->arg[2]);
472 }
473 else {
474 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Expect -A>gradientfile>+f<format>\n");
475 n_errors++;
476 }
477 break;
478 }
479 /* New syntax */
480 if ((c = strstr (opt->arg, "+f")) == NULL) {
481 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Expect -A>gradientfile>+f<format>\n");
482 n_errors++;
483 }
484 else {
485 Ctrl->A.mode = (int)(c[2] - '0');
486 c[0] = '\0'; /* Temporarily chop off the modifier */
487 if (opt->arg[0] == 0) {
488 GMT_Report (API, GMT_MSG_ERROR, "Option -A: No file given\n");
489 n_errors++;
490 }
491 else
492 Ctrl->A.file = strdup (opt->arg);
493 c[0] = '+'; /* Restore the modifier */
494 }
495 break;
496 case 'C': /* Solve by SVD */
497 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
498 Ctrl->C.active = true;
499 if (strchr (opt->arg, '/') && strstr (opt->arg, "+f") == NULL) { /* Old-style file deprecated specification */
500 if (gmt_M_compat_check (API->GMT, 5)) { /* OK */
501 sscanf (&opt->arg[k], "%lf/%s", &Ctrl->C.value, p);
502 Ctrl->C.file = strdup (p);
503 }
504 else {
505 GMT_Report (API, GMT_MSG_ERROR, "Option -C: Expected -C[[n|r|v]<cut>[%]][+c][+f<file>][+i][+n]\n");
506 n_errors++;
507 }
508 break; /* No modifiers for the deprecated syntax */
509 }
510 switch (opt->arg[0]) { /* First check directives (Default is no directive [0] which is the same as -Cr) */
511 case 'n': Ctrl->C.mode = GMT_SVD_EIGEN_NUMBER_CUTOFF; k = 1; break;
512 case 'r': Ctrl->C.mode = GMT_SVD_EIGEN_RATIO_CUTOFF; k = 1; break;
513 case 'v': Ctrl->C.mode = GMT_SVD_EIGEN_VARIANCE_CUTOFF; k = 1; break;
514 default: /* No directive, probably part of a number of modifier, just ignore */
515 k = 0; break;
516 }
517 if ((c = gmt_first_modifier (GMT, &opt->arg[k], "cifmMn"))) { /* Process any modifiers */
518 pos = 0; /* Reset to start of new word */
519 while (gmt_getmodopt (GMT, 'C', c, "cifmMn", &pos, p, &n_errors) && n_errors == 0) {
520 switch (p[0]) {
521 case 'c': Ctrl->C.history |= GMT_SVD_CUMULATIVE; break;
522 case 'i': Ctrl->C.history |= GMT_SVD_INCREMENTAL; break;
523 case 'f': Ctrl->C.file = strdup (&p[1]); break;
524 case 'm': Ctrl->C.history = GMT_SVD_INCREMENTAL; break;
525 case 'M': Ctrl->C.history = GMT_SVD_CUMULATIVE; break;
526 case 'n': Ctrl->C.dryrun = true; break;
527 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
528 }
529 }
530 c[0] = '\0';
531 }
532 if (opt->arg[k]) { /* See if we got any value argument */
533 if (strchr (opt->arg, '%')) { /* Got percentages of largest eigenvalues */
534 Ctrl->C.value = 0.01 * atof (&opt->arg[k]);
535 if (Ctrl->C.mode == GMT_SVD_EIGEN_NUMBER_CUTOFF) Ctrl->C.mode = GMT_SVD_EIGEN_PERCENT_CUTOFF; /* else it is GMT_SVD_EIGEN_VARIANCE_CUTOFF */
536 if (Ctrl->C.value < 0.0 || Ctrl->C.value > 1.0) {
537 GMT_Report (API, GMT_MSG_ERROR, "Option -C: Percentages must be in 0-100%% range!\n");
538 n_errors++;
539 }
540 }
541 else { /* Got regular cutoff value */
542 Ctrl->C.value = atof (&opt->arg[k]);
543 if (Ctrl->C.mode == GMT_SVD_EIGEN_NUMBER_CUTOFF && Ctrl->C.value >= 0.0 && Ctrl->C.value < 1.0) /* Old style fraction given instead */
544 Ctrl->C.mode = GMT_SVD_EIGEN_PERCENT_CUTOFF;
545 }
546 }
547 if (Ctrl->C.value < 0.0) Ctrl->C.dryrun = true, Ctrl->C.value = 0.0; /* Check for deprecated syntax giving negative value */
548 break;
549 case 'D':
550 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
551 if (gmt_M_compat_check (API->GMT, 6) && strlen (opt->arg) == 1) { /* Old -D<mode> option supported for backwards compatibility (now -Z) */
552 Ctrl->Z.active = true;
553 Ctrl->Z.mode = atoi (opt->arg); /* Since I added 0 to be 1-D later so now this is mode -1 */
554 }
555 else { /* Give grid or cube information */
556 Ctrl->D.active = true;
557 Ctrl->D.information = strdup (opt->arg);
558 }
559 break;
560 case 'E': /* Evaluate misfit -E[<file>]*/
561 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
562 Ctrl->E.active = true;
563 if (opt->arg) {
564 Ctrl->E.file = strdup (opt->arg);
565 Ctrl->E.mode = 1;
566 }
567 break;
568 case 'G': /* Output file */
569 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
570 Ctrl->G.active = true;
571 Ctrl->G.file = strdup (opt->arg);
572 break;
573 case 'I': /* Table or grid spacings */
574 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
575 Ctrl->I.active = true;
576 k = gmt_getincn (GMT, opt->arg, Ctrl->I.inc, 3);
577 if (k < 1) {
578 gmt_inc_syntax (GMT, 'I', 1);
579 n_errors++;
580 }
581 if (Ctrl->I.inc[GMT_Y] == 0.0) Ctrl->I.inc[GMT_Y] = Ctrl->I.inc[GMT_X];
582 if (Ctrl->I.inc[GMT_Z] == 0.0) Ctrl->I.inc[GMT_Z] = Ctrl->I.inc[GMT_X];
583 break;
584 case 'L': /* Leave trend alone */
585 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
586 Ctrl->L.active = true;
587 break;
588 case 'M': /* Read or write list of Green's function forces [NOT IMPLEMENTED YET] */
589 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
590 Ctrl->M.active = true;
591 Ctrl->M.file = strdup (opt->arg);
592 break;
593 case 'N': /* Output locations */
594 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
595 Ctrl->N.active = true;
596 if (opt->arg[0]) Ctrl->N.file = strdup (opt->arg);
597 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->N.file))) n_errors++;
598 break;
599 case 'Q': /* Directional derivative */
600 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
601 Ctrl->Q.active = true;
602 if (strchr (opt->arg, '/')) { /* Got 3-D vector components */
603 n_items = sscanf (opt->arg, "%lf/%lf/%lf", &Ctrl->Q.dir[0], &Ctrl->Q.dir[1], &Ctrl->Q.dir[2]);
604 if (n_items != 3) {
605 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Append azimuth (2-D) or x/y/z components (3-D)\n");
606 n_errors++;
607 }
608 gmt_normalize3v (GMT, Ctrl->Q.dir); /* Normalize to unit vector */
609 }
610 else if (opt->arg[0]) /* 2-D azimuth */
611 Ctrl->Q.az = atof(opt->arg);
612 else {
613 GMT_Report (API, GMT_MSG_ERROR, "Option -Q: Append azimuth (2-D) or x/y/z components (3-D)\n");
614 n_errors++;
615 }
616 break;
617 case 'S': /* Spline selection */
618 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
619 Ctrl->S.active = true;
620 Ctrl->S.arg = strdup (opt->arg);
621 switch (opt->arg[0]) {
622 case 'l': /* Cartesian linear spline in 1-D or 2-D (bilinear) */
623 Ctrl->S.mode = LINEAR_1D;
624 break;
625 case 'c': /* Cartesian minimum curvature spline */
626 Ctrl->S.mode = SANDWELL_1987_1D;
627 break;
628 case 't': /* Cartesian minimum curvature spline in tension */
629 Ctrl->S.mode = WESSEL_BERCOVICI_1998_1D;
630 if (strchr (opt->arg, '/'))
631 sscanf (&opt->arg[1], "%lf/%lf", &Ctrl->S.value[0], &Ctrl->S.value[1]);
632 else
633 Ctrl->S.value[0] = atof (&opt->arg[1]);
634 break;
635 case 'r': /* Regularized minimum curvature spline in tension in 2-D or 3-D */
636 Ctrl->S.mode = MITASOVA_MITAS_1993_2D;
637 if (strchr (opt->arg, '/'))
638 sscanf (&opt->arg[1], "%lf/%lf", &Ctrl->S.value[0], &Ctrl->S.value[1]);
639 else
640 Ctrl->S.value[0] = atof (&opt->arg[1]);
641 break;
642 case 'p': /* Spherical minimum curvature spline */
643 Ctrl->S.mode = PARKER_1994;
644 break;
645 case 'Q': /* Spherical minimum curvature spline in tension */
646 if (gmt_M_compat_check (GMT, 4)) {
647 GMT_Report (API, GMT_MSG_COMPAT, "Option -SQ is deprecated; see -Sq syntax instead.\n");
648 Ctrl->S.mode = WESSEL_BECKER_2008;
649 if (strchr (opt->arg, '/')) {
650 n_items = sscanf (&opt->arg[1], "%lf/%lf/%lf/%lf", &Ctrl->S.value[0], &Ctrl->S.value[1], &Ctrl->S.rval[0], &Ctrl->S.rval[1]);
651 if (n_items == 2) {
652 Ctrl->S.rval[0] = -1.0;
653 Ctrl->S.rval[1] = +1.0;
654 }
655 }
656 else
657 Ctrl->S.value[0] = atof (&opt->arg[1]);
658 }
659 else {
660 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Append c|l|t|g|p|q\n");
661 n_errors++;
662 }
663 break;
664 case 'q': /* Spherical minimum curvature spline in tension */
665 Ctrl->S.mode = WESSEL_BECKER_2008;
666 Ctrl->S.value[0] = atof (&opt->arg[1]);
667 if (Ctrl->S.value[0] == 0.0) /* Switch to Parker_1994 since tension is zero */
668 Ctrl->S.mode = PARKER_1994;
669 if ((c = strchr (opt->arg, '+')) != NULL) {
670 while (gmt_strtok (c, "+", &pos, p)) {
671 switch (p[0]) {
672 case 'e': Ctrl->S.value[2] = atof (&p[1]); break; /* Change the truncation error limit */
673 case 'n': Ctrl->S.value[3] = atof (&p[1]); break; /* Change the number of nodes for the spline lookup */
674 case 'l': Ctrl->S.rval[0] = atof (&p[1]); break; /* Min value for spline, undocumented for testing only */
675 case 'u': Ctrl->S.rval[1] = atof (&p[1]); break; /* Max value for spline, undocumented for testing only */
676 default:
677 GMT_Report (API, GMT_MSG_ERROR, "Option -Sq: Unknown modifier %s\n", p);
678 n_errors++;
679 break;
680 }
681 }
682 }
683 break;
684 default:
685 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Append c|l|t|g|p|q[<params>]\n");
686 n_errors++;
687 break;
688 }
689 break;
690 case 'T': /* Input mask grid */
691 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
692 Ctrl->T.active = true;
693 if (opt->arg[0]) Ctrl->T.file = strdup (opt->arg);
694 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->T.file)))
695 n_errors++;
696 else { /* Obtain -R -I -r from file */
697 struct GMT_GRID *G = NULL;
698 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->T.file, NULL)) == NULL) { /* Get header only */
699 return (API->error);
700 }
701 Ctrl->R3.range[0] = G->header->wesn[XLO]; Ctrl->R3.range[1] = G->header->wesn[XHI];
702 Ctrl->R3.range[2] = G->header->wesn[YLO]; Ctrl->R3.range[3] = G->header->wesn[YHI];
703 Ctrl->R3.inc[GMT_X] = G->header->inc[GMT_X]; Ctrl->R3.inc[GMT_Y] = G->header->inc[GMT_Y];
704 Ctrl->R3.offset = G->header->registration;
705 Ctrl->R3.dimension = 2;
706 Ctrl->R3.mode = true;
707 if (GMT_Destroy_Data (API, &G) != GMT_NOERROR) {
708 return (API->error);
709 }
710 GMT->common.R.active[RSET] = true;
711 }
712 break;
713 case 'W': /* Expect data uncertainty or weights in last column */
714 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
715 Ctrl->W.active = true;
716 if (opt->arg[0] == 'w') Ctrl->W.mode = GSP_GOT_WEIGHTS; /* Got weights instead of sigmas */
717 break;
718 case 'Z': /* Distance mode */
719 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
720 Ctrl->Z.active = true;
721 Ctrl->Z.mode = atoi (opt->arg); /* Since I added 0 to be 1-D later so now this is mode -1 */
722 break;
723 #ifdef DEBUG
724 case '0': /* Dump matrices */
725 Ctrl->debug.active = true;
726 break;
727 case '+': /* Turn on TEST mode */
728 TEST = true;
729 break;
730 #endif
731 default: /* Report bad options */
732 n_errors += gmt_default_error (GMT, opt->option);
733 break;
734 }
735 }
736
737 if (Ctrl->M.active) { /* Determine if this is for reading or writing Green's function forces */
738 if (Ctrl->S.active) /* Writing, since -S was given */
739 Ctrl->M.mode = GMT_OUT;
740 else if (gmt_access (GMT, Ctrl->M.file, F_OK)) {
741 GMT_Report (API, GMT_MSG_ERROR, "-M option given but file %s not found\n", Ctrl->M.file);
742 n_errors++;
743 }
744 else /* Read in previous Green's function forces */
745 Ctrl->M.mode = GMT_IN;
746 }
747
748 if (Ctrl->S.mode == WESSEL_BECKER_2008) { /* Check that nodes is an odd integer */
749 double fn = rint (Ctrl->S.value[3]);
750 int64_t n = lrint (fn);
751 if (!doubleAlmostEqual (Ctrl->S.value[3], fn) || ((n%2) == 0)) {
752 GMT_Report (API, GMT_MSG_ERROR, "Option -Sq option +n<N> modifier: <N> must be an odd integer\n");
753 n_errors++;
754 }
755 if (Ctrl->S.value[2] < 0.0 || Ctrl->S.value[2] > 1.0e-4) {
756 GMT_Report (API, GMT_MSG_ERROR, "Option -Sq option +e<err> modifier: <err> must be positive and < 1.0e-4\n");
757 n_errors++;
758 }
759 }
760 if (Ctrl->S.mode == PARKER_1994 || Ctrl->S.mode == WESSEL_BECKER_2008) Ctrl->Z.mode = 4; /* Automatically set */
761 dimension = (Ctrl->Z.mode == 0) ? 1 : ((Ctrl->Z.mode == 5) ? 3 : 2);
762 if (dimension == 2 && Ctrl->R3.mode) { /* Set -R via a gridfile */
763 /* Here, -R<grdfile> was used and we will use the settings supplied by the grid file (unless overridden) */
764 if (!Ctrl->I.active) { /* -I was not set separately; set indirectly */
765 Ctrl->I.inc[GMT_X] = Ctrl->R3.inc[GMT_X];
766 Ctrl->I.inc[GMT_Y] = Ctrl->R3.inc[GMT_Y];
767 Ctrl->I.active = true;
768 }
769 /* Here, -r means toggle the grids registration */
770 if (GMT->common.R.active[GSET]) {
771 GMT->common.R.active[GSET] = !Ctrl->R3.offset;
772 GMT->common.R.registration = !Ctrl->R3.offset;
773 }
774 else {
775 GMT->common.R.active[GSET] = Ctrl->R3.offset;
776 GMT->common.R.registration = Ctrl->R3.offset;
777 }
778 }
779
780 n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && gmt_access (GMT, Ctrl->A.file, R_OK), "Option -A: Cannot read file %s!\n", Ctrl->A.file);
781 n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && Ctrl->A.mode > 5, "Option -A: format must be in 0-5 range\n");
782 n_errors += gmt_M_check_condition (GMT, !(GMT->common.R.active[RSET] || Ctrl->N.active || Ctrl->T.active), "No output locations specified (use either [-R -I], -N, or -T)\n");
783 n_errors += gmt_M_check_condition (GMT, Ctrl->R3.mode && dimension != 2, "The -R<gridfile> or -T<gridfile> option only applies to 2-D gridding\n");
784 n_errors += gmt_M_check_condition (GMT, Ctrl->C.history && dimension != 2, "The -C +c+i modifiers only apply to 2-D gridding\n");
785 n_errors += gmt_M_check_condition (GMT, Ctrl->C.history && strchr (Ctrl->G.file, '%') == NULL && strchr (Ctrl->G.file, '.') == NULL, "Option -G: When -C +i|c is used your grid file must have an extension\n");
786 n_errors += gmt_M_check_condition (GMT, Ctrl->C.history && Ctrl->E.active && Ctrl->E.mode == 0, "Option -E: When -C +i|c is used you must supply a file via -E\n");
787 #ifdef DEBUG
788 n_errors += gmt_M_check_condition (GMT, !TEST && !Ctrl->N.active && Ctrl->R3.dimension != dimension, "The -R and -Z options disagree on the dimension\n");
789 #else
790 n_errors += gmt_M_check_condition (GMT, !Ctrl->N.active && Ctrl->R3.dimension != dimension, "The -R and -Z options disagree on the dimension\n");
791 #endif
792 n_errors += gmt_check_binary_io (GMT, dimension + 1);
793 n_errors += gmt_M_check_condition (GMT, Ctrl->S.value[0] < 0.0 || Ctrl->S.value[0] >= 1.0, "Option -S: Tension must be in range 0 <= t < 1\n");
794 n_errors += gmt_M_check_condition (GMT, !(Ctrl->S.mode == PARKER_1994 || Ctrl->S.mode == WESSEL_BECKER_2008) && Ctrl->Z.mode == 3, "Option -Sc|t|r: Cannot select -Z3\n");
795 n_errors += gmt_M_check_condition (GMT, Ctrl->S.mode == LINEAR_1D && Ctrl->Z.mode > 3, "Option -Sl: Cannot select -Z4 or higher\n");
796 n_errors += gmt_M_check_condition (GMT, Ctrl->I.active && (Ctrl->I.inc[GMT_X] <= 0.0 || (dimension > 1 && Ctrl->I.inc[GMT_Y] <= 0.0) || (dimension == 3 && Ctrl->I.inc[GMT_Z] <= 0.0)), "Option -I: Must specify positive increment(s)\n");
797 n_errors += gmt_M_check_condition (GMT, dimension == 2 && !Ctrl->N.active && !(Ctrl->G.active || Ctrl->G.file), "Option -G: Must specify output grid file name\n");
798 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->C.dryrun && !Ctrl->C.file, "Option -C: Must specify file name for eigenvalues if +n is set");
799 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && !Ctrl->T.file, "Option -T: Must specify mask grid file name\n");
800 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && dimension != 2, "Option -T: Only applies to 2-D gridding\n");
801 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->N.file, "Option -N: Must specify node file name\n");
802 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->N.file && gmt_access (GMT, Ctrl->N.file, R_OK), "Option -N: Cannot read file %s!\n", Ctrl->N.file);
803 n_errors += gmt_M_check_condition (GMT, (Ctrl->I.active + GMT->common.R.active[RSET]) == 1 && dimension == 2, "Must specify -R, -I, [-r], -G for gridding\n");
804 n_errors += gmt_M_check_condition (GMT, dimension == 3 && Ctrl->G.active && API->external && strchr (Ctrl->G.file, '%'), "Option -G: Cannot contain format-specifiers when not used on the command line\n");
805
806 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
807 }
808
809 #ifdef DEBUG
810 /* Dump a table of x, G, dGdx for test purposes [requires option -+ and compilation with -DDEBUG] */
greenspline_dump_green(struct GMT_CTRL * GMT,double (* G)(struct GMT_CTRL *,double,double *,struct GREENSPLINE_LOOKUP *),double (* D)(struct GMT_CTRL *,double,double *,struct GREENSPLINE_LOOKUP *),double par[],double x0,double x1,int N,struct GREENSPLINE_LOOKUP * Lz,struct GREENSPLINE_LOOKUP * Lg)811 GMT_LOCAL void greenspline_dump_green (struct GMT_CTRL *GMT, double (*G) (struct GMT_CTRL *, double, double *, struct GREENSPLINE_LOOKUP *), double (*D) (struct GMT_CTRL *, double, double *, struct GREENSPLINE_LOOKUP *), double par[], double x0, double x1, int N, struct GREENSPLINE_LOOKUP *Lz, struct GREENSPLINE_LOOKUP *Lg) {
812 int i;
813 double x, dx, dy, y, t, ry, rdy;
814 double min_y, max_y, min_dy, max_dy;
815
816 min_y = min_dy = DBL_MAX;
817 max_y = max_dy = -DBL_MAX;
818
819 dx = (x1 - x0) / (N - 1);
820 for (i = 0; i < N; i++) {
821 x = x0 + i * dx;
822 y = G (GMT, x, par, Lz);
823 dy = D (GMT, x, par, Lg);
824 if (y < min_y) min_y = y;
825 if (y > max_y) max_y = y;
826 if (dy < min_dy) min_dy = dy;
827 if (dy > max_dy) max_dy = dy;
828 }
829 ry = max_y - min_y;
830 rdy = max_dy - min_dy;
831 for (i = 0; i < N; i++) {
832 x = x0 + i * dx;
833 t = (x0 < 0.0) ? acosd (x) : x;
834 y = G (GMT, x, par, Lz);
835 dy = D (GMT, x, par, Lg);
836 dy = (rdy > 0.0) ? (dy - min_dy)/rdy : 1.0;
837 printf ("%g\t%g\t%g\t%g\n", x, (y - min_y) / ry, dy, t);
838 }
839 }
840 #endif
841
842 /* Below are all the individual Green's functions. Note that most of them take an argument
843 * that is unused except in the spline lookup version of WB_08. */
844
845 /*---------------------- ONE DIMENSION ---------------------- */
846 /* Basic linear spline (bilinear in 2-D) */
847
greenspline_spline1d_linear(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)848 GMT_LOCAL double greenspline_spline1d_linear (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
849 /* Dumb linear spline */
850 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
851 return (r); /* Just regular spline; par not used */
852 }
853
greenspline_grad_spline1d_linear(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)854 GMT_LOCAL double greenspline_grad_spline1d_linear (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
855 /* d/dr of r is 1 */
856 gmt_M_unused(GMT); gmt_M_unused(r); gmt_M_unused(par); gmt_M_unused(unused);
857 return (1.0);
858 }
859
860 /* greenspline_spline1d_sandwell computes the Green function for a 1-d spline
861 * as per Sandwell [1987], G(r) = r^3. All r must be >= 0.
862 */
863
greenspline_spline1d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)864 GMT_LOCAL double greenspline_spline1d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
865 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
866 if (r == 0.0) return (0.0);
867
868 return (pow (r, 3.0)); /* Just regular spline; par not used */
869 }
870
greenspline_grad_spline1d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)871 GMT_LOCAL double greenspline_grad_spline1d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
872 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
873 return (r); /* Just regular spline; par not used */
874 }
875
greenspline_spline1d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)876 GMT_LOCAL double greenspline_spline1d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
877 /* greenspline_spline1d_Wessel_Bercovici computes the Green function for a 1-d spline
878 * in tension as per Wessel and Bercovici [1988], G(u) = exp(-u) + u - 1,
879 * where u = par[0] * r and par[0] = sqrt (t/(1-t)).
880 * All r must be >= 0. par[0] = c
881 */
882 double cx;
883 gmt_M_unused(GMT); gmt_M_unused(unused);
884
885 if (r == 0.0) return (0.0);
886
887 cx = par[0] * r;
888 return (exp (-cx) + cx - 1.0);
889 }
890
greenspline_grad_spline1d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)891 GMT_LOCAL double greenspline_grad_spline1d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
892 double cx;
893 gmt_M_unused(GMT); gmt_M_unused(unused);
894
895 if (r == 0.0) return (0.0);
896
897 cx = par[0] * r;
898 return (1.0 - exp (-cx));
899 }
900
901 /*---------------------- TWO DIMENSIONS ---------------------- */
902
greenspline_spline2d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)903 GMT_LOCAL double greenspline_spline2d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
904 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
905 if (r == 0.0) return (0.0);
906
907 return (r * r * (log (r) - 1.0)); /* Just regular spline; par not used */
908 }
909
greenspline_grad_spline2d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)910 GMT_LOCAL double greenspline_grad_spline2d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
911 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
912 if (r == 0.0) return (0.0);
913
914 return (r * (2.0 * log (r) - 1.0)); /* Just regular spline; par not used */
915 }
916
917 /* greenspline_spline2d_Wessel_Bercovici computes the Green function for a 2-d spline
918 * in tension as per Wessel and Bercovici [1988], G(u) = K(u) - log(u),
919 * where u = par[0] * r and par[0] = sqrt (t/(1-t)).
920 * K is the modified Bessel function K of order zero.
921 * All r must be >= 0.
922 * par[0] = c
923 * par[1] = 2/c
924 */
925
greenspline_spline2d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)926 GMT_LOCAL double greenspline_spline2d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
927 double y, z, cx, t, g;
928 gmt_M_unused(GMT); gmt_M_unused(unused);
929
930 if (r == 0.0) return (0.0);
931
932 cx = par[0] * r;
933 if (r <= par[1]) {
934 y = 0.25 * (t = cx * cx);
935 z = t / 14.0625;
936 g = (-log(0.5*cx) * (z * (3.5156229 + z * (3.0899424 + z * (1.2067492 + z * (0.2659732 +
937 z * (0.360768e-1 + z * 0.45813e-2))))))) + (y * (0.42278420 + y * (0.23069756 +
938 y * (0.3488590e-1 + y * (0.262698e-2 + y * (0.10750e-3 + y * 0.74e-5))))));
939 }
940 else {
941 y = par[1] / r;
942 g = (exp (-cx) / sqrt (cx)) * (1.25331414 + y * (-0.7832358e-1 + y * (0.2189568e-1 +
943 y * (-0.1062446e-1 + y * (0.587872e-2 + y * (-0.251540e-2 + y * 0.53208e-3))))))
944 + log (cx) - M_LOG_2 + M_GAMMA;
945 }
946 return (g);
947 }
948
greenspline_grad_spline2d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)949 GMT_LOCAL double greenspline_grad_spline2d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
950 double y, z, cx, t, dgdr;
951 gmt_M_unused(GMT); gmt_M_unused(unused);
952
953 if (r == 0.0) return (0.0);
954
955 cx = par[0] * r;
956 if (r <= par[1]) {
957 y = 0.25 * (t = cx * cx);
958 z = t / 14.0625;
959 dgdr = -((log(0.5*cx) * (cx * (0.5 + z * (0.87890594 + z * (0.51498869 + z * (0.15084934 +
960 z * (0.2658733e-1 + z * (0.301532e-2 + z * 0.32411e-3)))))))) + (1.0/cx) * (y * (0.15443144 +
961 y * (-0.67278579 + y * (-0.18156897 + y * (-0.1919402e-1 + y * (-0.110404e-2 + y * (-0.4686e-4))))))));
962 }
963 else {
964 y = par[1] / r;
965 dgdr = 0.5*y - ((exp (-cx) / sqrt (cx)) * (1.25331414 + y * (0.23498619 + y * (-0.3655620e-1 +
966 y * (0.1504268e-1 + y * (-0.780353e-2 + y * (0.325614e-2 + y * (-0.68245e-3))))))));
967 }
968 return (dgdr*par[0]);
969 }
970
971 /* greenspline_spline2d_Mitasova_Mitas computes the regularized Green function for a 2-d spline
972 * in tension as per Mitasova and Mitas [1993], G(u) = log (u) + Ei(u),
973 * where u = par[1] * r^2. All r must be >= 0.
974 * par[0] = phi (par[0] unused by function)
975 * par[1] = phi^2/4
976 */
977
greenspline_spline2d_Mitasova_Mitas(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)978 GMT_LOCAL double greenspline_spline2d_Mitasova_Mitas (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
979 double x, z, g, En, Ed;
980 gmt_M_unused(GMT); gmt_M_unused(unused);
981
982 if (r == 0.0) return (0.0);
983
984 x = z = par[1] * r * r;
985 if (z <= 1.0) {
986 g = 0.99999193 * x;
987 x *= x;
988 g -= 0.24991055 * x;
989 x *= x;
990 g += 0.05519968 * x;
991 x *= x;
992 g -= 0.00976004 * x;
993 x *= x;
994 g += 0.00107857 * x;
995 }
996 else {
997 g = log (x) + M_GAMMA;
998 En = 0.2677737343 + 8.6347608925 * x;
999 Ed = 3.9584869228 + 21.0996530827 * x;
1000 x *= x;
1001 En += 18.0590169730 * x;
1002 Ed += 25.6329561486 * x;
1003 x *= x;
1004 En += 8.5733287401 * x;
1005 Ed += 9.5733223454 * x;
1006 x *= x;
1007 En += x;
1008 Ed += x;
1009 g += (En / Ed) / (z * exp(z));
1010 }
1011 return (g);
1012 }
1013
greenspline_grad_spline2d_Mitasova_Mitas(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1014 GMT_LOCAL double greenspline_grad_spline2d_Mitasova_Mitas (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1015 double u, dgdr;
1016 gmt_M_unused(GMT); gmt_M_unused(unused);
1017
1018 if (r == 0.0) return (0.0);
1019
1020 u = par[1] * r * r;
1021 dgdr = 2.0 * (1.0 - exp (-u))/r;
1022 return (dgdr);
1023 }
1024
1025 /*---------------------- TWO DIMENSIONS (SPHERE) ---------------------- */
1026
1027 /* greenspline_spline2d_Parker computes the Green function for a 2-d surface spline
1028 * as per Parker [1994], G(x) = dilog(),
1029 * where x is cosine of distances. All x must be -1 <= x <= +1.
1030 * Parameters passed are:
1031 * par[0] = 6/M_PI^2 (to normalize results)
1032 */
1033
greenspline_spline2d_Parker(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * unused)1034 GMT_LOCAL double greenspline_spline2d_Parker (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *unused) {
1035 /* Normalized to 0-1 */
1036 gmt_M_unused(unused);
1037 if (x == +1.0) return (1.0);
1038 if (x == -1.0) return (0.0);
1039 return (gmt_dilog (GMT, 0.5 - 0.5 * x) * par[0]);
1040 }
1041
greenspline_grad_spline2d_Parker(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * unused)1042 GMT_LOCAL double greenspline_grad_spline2d_Parker (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *unused) {
1043 /* Normalized to 0-1 */
1044 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
1045 if (x == +1.0 || x == -1.0) return (0.0);
1046 return (log (0.5 - 0.5 * x) * sqrt ((1.0 - x) / (1.0 + x)));
1047 }
1048
greenspline_series_prepare(struct GMT_CTRL * GMT,double p,unsigned int L,struct GREENSPLINE_LOOKUP * Lz,struct GREENSPLINE_LOOKUP * Lg)1049 GMT_LOCAL void greenspline_series_prepare (struct GMT_CTRL *GMT, double p, unsigned int L, struct GREENSPLINE_LOOKUP *Lz, struct GREENSPLINE_LOOKUP *Lg) {
1050 /* Precalculate Legendre series terms involving various powers/ratios of l and p */
1051 unsigned int l;
1052 double pp, t1, t2;
1053 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Precalculate max %u terms for Legendre summation\n", L+1);
1054 Lz->A = gmt_M_memory (GMT, NULL, L+1, double);
1055 Lz->B = gmt_M_memory (GMT, NULL, L+1, double);
1056 Lz->C = gmt_M_memory (GMT, NULL, L+1, double);
1057 if (Lg) Lg->A = gmt_M_memory (GMT, NULL, L+1, double);
1058 pp = p * p;
1059 for (l = 1; l <= L; l++) {
1060 t1 = l + 1.0;
1061 t2 = 2.0 * l + 1.0;
1062 Lz->A[l] = t2 * pp / (l*t1*(l*t1+pp));
1063 Lz->B[l] = t2 / t1;
1064 Lz->C[l] = l / t1;
1065 if (Lg) Lg->A[l] = t2 * pp / (t1*(l*t1+pp));
1066 }
1067 if (Lg) Lg->B = Lz->B; /* Lg needs the same B,C coefficients as Lz, so just share them via the link */
1068 if (Lg) Lg->C = Lz->C;
1069 }
1070
greenspline_get_max_L(struct GMT_CTRL * GMT,double p,double err)1071 GMT_LOCAL unsigned int greenspline_get_max_L (struct GMT_CTRL *GMT, double p, double err) {
1072 /* Return max L needed given p and truncation err for Parker's simplified loop expression */
1073 gmt_M_unused(GMT);
1074 return ((unsigned int)lrint (p / sqrt (err) + 10.0));
1075 }
1076
greenspline_free_lookup(struct GMT_CTRL * GMT,struct GREENSPLINE_LOOKUP ** Lptr,unsigned int mode)1077 GMT_LOCAL void greenspline_free_lookup (struct GMT_CTRL *GMT, struct GREENSPLINE_LOOKUP **Lptr, unsigned int mode) {
1078 /* Free all items allocated under the lookup structures.
1079 * mode = 0 means Lz and mode = 1 means Lg; the latter has no B & C arrays */
1080 struct GREENSPLINE_LOOKUP *L = *Lptr;
1081 if (L == NULL) return; /* Nothing to free */
1082 gmt_M_free (GMT, L->y);
1083 gmt_M_free (GMT, L->c);
1084 gmt_M_free (GMT, L->A);
1085 if (mode == 0) { /* Only Lz has A,B,C while Lg borrows B,C */
1086 gmt_M_free (GMT, L->B);
1087 gmt_M_free (GMT, L->C);
1088 }
1089 gmt_M_free (GMT, L);
1090 *Lptr = NULL;
1091 }
1092
greenspline_get_L(double x,double p,double err)1093 GMT_LOCAL unsigned int greenspline_get_L (double x, double p, double err) {
1094 /* Determines the truncation order L_max given x, p, and err.
1095 * See ParkerNotesJan2013.pdf in gurudocs for explanations and derivation */
1096 unsigned int L_max;
1097 double s, c, pp, Lf, gam, lam;
1098 /* Get all the trig functions needed; x is cos(theta), but we need */
1099 /* s = sin (theta/2); c=cos(theta/2); */
1100 /* cos(theta) = x = 1.0 - 2.0 *s^2, so */
1101 s = sqrt (0.5 * (1 - x));
1102 c = sqrt (0.5 * (1 + x));
1103 pp = p * p;
1104
1105 /* Approximate the highest order L_max we need to sum given specified err limit */
1106 if (s == 0.0)
1107 Lf = p / sqrt (err);
1108 else if (c == 0.0)
1109 Lf = pow (pp / err, 0.333333333);
1110 else {
1111 gam = 0.00633262522 * pow (pp * s * s / err, 2.0) / c;
1112 if (gam <= 1.0) lam = pow (gam / (1.0 + pow (gam, 0.4)), 0.2);
1113 else lam = pow (gam / (1.0 + 1.0 / pow (gam, 0.285714286)), 0.142857143);
1114 Lf = 1.75 * lam / s;
1115 }
1116
1117 Lf = MIN (p / sqrt (err), Lf);
1118 L_max = (unsigned int)lrint (Lf + 10.0);
1119 return (L_max);
1120 }
1121
greenspline_spline2d_Wessel_Becker_Revised(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * Lz)1122 GMT_LOCAL double greenspline_spline2d_Wessel_Becker_Revised (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *Lz) {
1123 /* Evaluate g = M_PI * Pv(-x)/sin (v*x) - log (1-x) via series approximation */
1124 unsigned int L_max, l;
1125 double P0, P1, P2, S;
1126 gmt_M_unused(GMT);
1127
1128 L_max = greenspline_get_L (x, par[0], par[1]); /* Highest order needed in sum given x, p, and err */
1129
1130 /* pp = par[0] * par[0]; */
1131 P0 = 1.0; /* Initialize the P0 and P1 Legendre polynomials */
1132 P1 = x;
1133 S = par[2]; /* The constant l = 0 term was computed during setup */
1134
1135 /* Sum the Legendre series */
1136 for (l = 1; l <= L_max; l++) {
1137 /* All the coeffs in l have been precomputed by greenspline_series_prepare.
1138 * S += (2*l+1)*pp/(l*(l+1)*(l*(l+1)+pp)) * P1;
1139 * P2 = x*(2*l+1)/(l+1) * P1 - l*P0/(l+1); */
1140 S += Lz->A[l] * P1; /* Update sum */
1141 P2 = x * Lz->B[l] * P1 - Lz->C[l] * P0; /* Get next Legendre polynomial value */
1142 P0 = P1;
1143 P1 = P2;
1144 }
1145
1146 return (S);
1147 }
1148
greenspline_grad_spline2d_Wessel_Becker_Revised(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * Lg)1149 GMT_LOCAL double greenspline_grad_spline2d_Wessel_Becker_Revised (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *Lg) {
1150 /* Evaluate g = -M_PI * (v+1)*[x*Pv(-x)+Pv+1(-x)]/(sin (v*x)*sin(theta)) - 1/(1-x) via series approximation */
1151 unsigned int L_max, l;
1152 double sin_theta, P0, P1, P2, S;
1153 gmt_M_unused(GMT);
1154
1155 if (fabs(x) == 1.0) return 1.0;
1156 L_max = greenspline_get_L (x, par[0], par[1]); /* Highest order needed in sum given x, p, and err */
1157 sin_theta = sqrt (1.0 - x * x);
1158 P0 = 1.0; /* Initialize the P0 and P1 Legendre polynomials */
1159 P1 = x;
1160 S = 0.0; /* Initialize sum */
1161
1162 /* Sum the Legendre series */
1163 for (l = 1; l <= L_max; l++) {
1164 /* All the coeffs in l have been precomputed by greenspline_series_prepare.
1165 * S += (2*l+1)*pp/((l+1)*(l*(l+1)+pp)) * (P1 * x - P0);
1166 * P2 = x*(2*l+1)/(l+1) * P1 - l*P0/(l+1); */
1167 S += Lg->A[l] * (P1 * x - P0); /* Update sum */
1168 P2 = x * Lg->B[l] * P1 - Lg->C[l] * P0; /* Get next Legendre polynomial value */
1169 P0 = P1;
1170 P1 = P2;
1171 }
1172 S /= sin_theta;
1173
1174 return (S);
1175 }
1176
1177 /* Given the lookup tables, this is how we use these functions
1178 * Here, par[7] is number of points in spline
1179 * par[8] is spline spacing dx
1180 * par[9] is inverse, 1/dx
1181 * par[10] is min x
1182 */
1183
greenspline_csplint(double * y,double * c,double b,double h2,uint64_t klo)1184 GMT_LOCAL double greenspline_csplint (double *y, double *c, double b, double h2, uint64_t klo) {
1185 /* Special version of support_greenspline_csplint where x is equidistant with spacing squared h2
1186 * and b is the fractional distance relative to x[klo], so x itself is not needed here. */
1187 uint64_t khi;
1188 double a, yp;
1189
1190 khi = klo + 1;
1191 a = 1.0 - b; /* Fractional distance from next node */
1192 yp = a * y[klo] + b * y[khi] + ((a*a*a - a) * c[klo] + (b*b*b - b) * c[khi]) * h2 / 6.0;
1193
1194 return (yp);
1195 }
1196
greenspline_spline2d_lookup(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * L)1197 GMT_LOCAL double greenspline_spline2d_lookup (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *L) {
1198 /* Given x, look up nearest node xx[k] <= x and do cubic spline interpolation */
1199 uint64_t k;
1200 double f, f0, df, y;
1201 gmt_M_unused(GMT);
1202
1203 f = (x - par[10]) * par[9]; /* Floating point index */
1204 f0 = floor (f);
1205 df = f - f0;
1206 k = lrint (f0);
1207 if (df == 0.0) return (L->y[k]); /* Right on a node */
1208 y = greenspline_csplint (L->y, L->c, df, par[4], k); /* Call special cubic spline evaluator */
1209 return (y);
1210 }
1211
greenspline_spline2d_Wessel_Becker_lookup(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * L)1212 GMT_LOCAL double greenspline_spline2d_Wessel_Becker_lookup (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *L) {
1213 return (greenspline_spline2d_lookup (GMT, x, par, L));
1214 }
1215
greenspline_grad_spline2d_Wessel_Becker_lookup(struct GMT_CTRL * GMT,double x,double par[],struct GREENSPLINE_LOOKUP * L)1216 GMT_LOCAL double greenspline_grad_spline2d_Wessel_Becker_lookup (struct GMT_CTRL *GMT, double x, double par[], struct GREENSPLINE_LOOKUP *L) {
1217 return (greenspline_spline2d_lookup (GMT, x, par, L));
1218 }
1219
greenspline_spline2d_Wessel_Becker_splineinit(struct GMT_CTRL * GMT,double par[],double * x,struct GREENSPLINE_LOOKUP * L)1220 GMT_LOCAL void greenspline_spline2d_Wessel_Becker_splineinit (struct GMT_CTRL *GMT, double par[], double *x, struct GREENSPLINE_LOOKUP *L) {
1221 /* Set up cubic spline interpolation given the precomputed x,y values of the function */
1222 gmt_M_unused(GMT);
1223 gmt_M_unused(par);
1224 L->c = gmt_M_memory (GMT, NULL, 3*L->n, double);
1225 gmtlib_cspline (GMT, x, L->y, L->n, L->c);
1226 }
1227
greenspline_spline2d_Wessel_Becker_init(struct GMT_CTRL * GMT,double par[],struct GREENSPLINE_LOOKUP * Lz,struct GREENSPLINE_LOOKUP * Lg)1228 GMT_LOCAL void greenspline_spline2d_Wessel_Becker_init (struct GMT_CTRL *GMT, double par[], struct GREENSPLINE_LOOKUP *Lz, struct GREENSPLINE_LOOKUP *Lg) {
1229 uint64_t i, n_columns;
1230 double *x = NULL;
1231 #ifdef DUMP
1232 FILE *fp = NULL;
1233 uint64_t n_out;
1234 double out[3];
1235 fp = fopen ("greenspline.b", "wb");
1236 n_out = (Lg) ? 3 : 2;
1237 #endif
1238 n_columns = lrint (par[7]);
1239 Lz->n = n_columns;
1240 x = gmt_M_memory (GMT, NULL, n_columns, double);
1241 Lz->y = gmt_M_memory (GMT, NULL, n_columns, double);
1242 if (Lg) {
1243 Lg->y = gmt_M_memory (GMT, NULL, n_columns, double);
1244 Lg->n = n_columns;
1245 }
1246 for (i = 0; i < n_columns; i++) {
1247 x[i] = par[10] + i * par[8];
1248 Lz->y[i] = greenspline_spline2d_Wessel_Becker_Revised (GMT, x[i], par, Lz);
1249 if (Lg) Lg->y[i] = greenspline_grad_spline2d_Wessel_Becker_Revised (GMT, x[i], par, Lg);
1250 #ifdef DUMP
1251 out[0] = x[i]; out[1] = Lz->y[i]; if (Lg) out[2] = Lg->y[i];
1252 fwrite (out, sizeof (double), n_out, fp);
1253 #endif
1254 }
1255 #ifdef DUMP
1256 fclose (fp);
1257 #endif
1258 greenspline_spline2d_Wessel_Becker_splineinit (GMT, par, x, Lz);
1259 if (Lg) {
1260 if (x[0] == -1.0) Lg->y[0] = 2.0*Lg->y[1] - Lg->y[2]; /* Linear interpolation from 2 nearest nodes */
1261 greenspline_spline2d_Wessel_Becker_splineinit (GMT, par, x, Lg);
1262 }
1263 gmt_M_free (GMT, x); /* Done with x array */
1264 }
1265
1266 /*---------------------- THREE DIMENSIONS ---------------------- */
1267
1268 /* greenspline_spline3d_sandwell computes the Green function for a 3-d spline
1269 * as per Sandwell [1987], G(r) = r. All r must be >= 0.
1270 */
1271
greenspline_spline3d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1272 GMT_LOCAL double greenspline_spline3d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1273 gmt_M_unused(GMT); gmt_M_unused(par); gmt_M_unused(unused);
1274 if (r == 0.0) return (0.0);
1275
1276 return (r); /* Just regular spline; par not used */
1277 }
1278
greenspline_grad_spline3d_sandwell(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1279 GMT_LOCAL double greenspline_grad_spline3d_sandwell (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1280 gmt_M_unused(GMT); gmt_M_unused(r); gmt_M_unused(par); gmt_M_unused(unused);
1281 return (1.0); /* Just regular spline; par not used */
1282 }
1283
greenspline_spline3d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1284 GMT_LOCAL double greenspline_spline3d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1285 /* greenspline_spline1d_Wessel_Bercovici computes the Green function for a 3-d spline
1286 * in tension as per Wessel and Bercovici [1988], G(u) = [exp(-u) -1]/u,
1287 * where u = par[0] * r and par[0] = sqrt (t/(1-t)).
1288 * All r must be >= 0. par[0] = c
1289 */
1290 double cx;
1291 gmt_M_unused(GMT); gmt_M_unused(unused);
1292
1293 if (r == 0.0) return (0.0);
1294
1295 cx = par[0] * r;
1296 return (((exp (-cx) - 1.0) / cx) + 1.0);
1297 }
1298
greenspline_grad_spline3d_Wessel_Bercovici(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1299 GMT_LOCAL double greenspline_grad_spline3d_Wessel_Bercovici (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1300 double cx;
1301 gmt_M_unused(GMT); gmt_M_unused(unused);
1302
1303 if (r == 0.0) return (0.0);
1304
1305 cx = par[0] * r;
1306 return ((1.0 - exp (-cx) * (cx + 1.0)) / (cx * r));
1307 }
1308
1309 /* greenspline_spline3d_Mitasova_Mitas computes the regularized Green function for a 3-d spline
1310 * in tension as per Mitasova and Mitas [1993], G(u) = erf (u/2)/u - 1/sqrt(pi),
1311 * where u = par[0] * r. All r must be >= 0. par[0] = phi
1312 */
1313
greenspline_spline3d_Mitasova_Mitas(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1314 GMT_LOCAL double greenspline_spline3d_Mitasova_Mitas (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1315 double x;
1316 gmt_M_unused(GMT); gmt_M_unused(unused);
1317
1318 if (r == 0.0) return (0.0);
1319
1320 x = par[0] * r;
1321 return ((erf (0.5 * x) / x) - M_INV_SQRT_PI);
1322 }
1323
greenspline_grad_spline3d_Mitasova_Mitas(struct GMT_CTRL * GMT,double r,double par[],struct GREENSPLINE_LOOKUP * unused)1324 GMT_LOCAL double greenspline_grad_spline3d_Mitasova_Mitas (struct GMT_CTRL *GMT, double r, double par[], struct GREENSPLINE_LOOKUP *unused) {
1325 double u, dgdr;
1326 gmt_M_unused(GMT); gmt_M_unused(unused);
1327
1328 if (r == 0.0) return (0.0);
1329
1330 u = par[0] * r;
1331 dgdr = ((u/M_SQRT_PI) * exp (-u*u) - erf (0.5 * u)) / (u * r);
1332 return (dgdr);
1333 }
1334
1335 /* GENERAL NUMERICAL FUNCTIONS */
1336
1337 /* Normalization parameters are stored in the coeff array which holds up to GSP_LENGTH terms
1338 * coeff[GSP_MEAN_X]: The mean x coordinate
1339 * coeff[GSP_MEAN_Y]: The mean y coordinate
1340 * coeff[GSP_MEAN_Z]: The mean w coordinate
1341 * coeff[GSP_SLP_X]: The linear x slope
1342 * coeff[GSP_SLP_Y]: The linear y slope
1343 * coeff[GSP_RANGE]: The largest |range| of the detrended data
1344 */
1345
greenspline_undo_normalization(double * X,double w_norm,unsigned int mode,double * coeff,unsigned int dim)1346 GMT_LOCAL double greenspline_undo_normalization (double *X, double w_norm, unsigned int mode, double *coeff, unsigned int dim) {
1347 if (mode & GREENSPLINE_NORM) w_norm *= coeff[GSP_RANGE]; /* Scale back up by residual data range (ir we normalized) */
1348 w_norm += coeff[GSP_MEAN_Z]; /* Add in mean data value plus minimum residual value (ir we normalized by range) */
1349 if (mode & GREENSPLINE_TREND) { /* Restore residual trend */
1350 w_norm += coeff[GSP_SLP_X] * (X[GMT_X] - coeff[GSP_MEAN_X]);
1351 if (dim == 2) w_norm += coeff[GSP_SLP_Y] * (X[GMT_Y] - coeff[GSP_MEAN_Y]);
1352 }
1353 return (w_norm);
1354 }
1355
greenspline_do_normalization_1d(struct GMTAPI_CTRL * API,double ** X,double * obs,uint64_t n,unsigned int mode,double * coeff)1356 GMT_LOCAL void greenspline_do_normalization_1d (struct GMTAPI_CTRL *API, double **X, double *obs, uint64_t n, unsigned int mode, double *coeff) {
1357 /* We always remove/restore the mean observation value. mode is a combination of bitflags that affects what we do:
1358 * Bit GREENSPLINE_TREND will also remove linear trend
1359 * Bit GREENSPLINE_NORM will normalize residuals by range
1360 */
1361
1362 uint64_t i;
1363 double d, min = DBL_MAX, max = -DBL_MAX;
1364
1365 gmt_M_memset (coeff, GSP_LENGTH, double);
1366 for (i = 0; i < n; i++) { /* Find mean w-value */
1367 coeff[GSP_MEAN_Z] += obs[i];
1368 if ((mode & GREENSPLINE_TREND) == 0) continue; /* No linear trend to model */
1369 coeff[GSP_MEAN_X] += X[i][GMT_X];
1370 }
1371 coeff[GSP_MEAN_Z] /= n;
1372
1373 if (mode & GREENSPLINE_TREND) { /* Solve for LS linear trend using deviations from (0, 0, 0) */
1374 double xx, zz, sxx, sxz;
1375 sxx = sxz = 0.0;
1376 coeff[GSP_MEAN_X] /= n;
1377 for (i = 0; i < n; i++) {
1378 xx = X[i][GMT_X] - coeff[GSP_MEAN_X];
1379 zz = obs[i] - coeff[GSP_MEAN_Z];
1380 sxx += (xx * xx);
1381 sxz += (xx * zz);
1382 }
1383 if (sxx != 0.0) coeff[GSP_SLP_X] = sxz/ sxx;
1384 }
1385
1386 /* Remove linear trend (or mean) */
1387
1388 for (i = 0; i < n; i++) { /* Get residuals and find range */
1389 obs[i] -= coeff[GSP_MEAN_Z]; /* Always remove the mean data value */
1390 if (mode & GREENSPLINE_TREND) obs[i] -= (coeff[GSP_SLP_X] * (X[i][GMT_X] - coeff[GSP_MEAN_X]));
1391 if (obs[i] < min) min = obs[i];
1392 if (obs[i] > max) max = obs[i];
1393 }
1394 if (mode & GREENSPLINE_NORM) { /* Normalize by range */
1395 coeff[GSP_RANGE] = MAX (fabs(min), fabs(max)); /* Determine range */
1396 d = (coeff[GSP_RANGE] == 0.0) ? 1.0 : 1.0 / coeff[GSP_RANGE];
1397 for (i = 0; i < n; i++) obs[i] *= d; /* Normalize 0-1 */
1398 }
1399
1400 /* Recover obs(x) = w_norm(x) * coeff[GSP_RANGE] + coeff[GSP_MEAN_Z] + coeff[GSP_SLP_X]*(x-coeff[GSP_MEAN_X]) */
1401 GMT_Report (API, GMT_MSG_INFORMATION, "1-D Normalization coefficients: zoff = %g slope = %g xmean = %g range = %g\n", coeff[GSP_MEAN_Z], coeff[GSP_SLP_X], coeff[GSP_MEAN_X], coeff[GSP_RANGE]);
1402 }
1403
greenspline_do_normalization(struct GMTAPI_CTRL * API,double ** X,double * obs,uint64_t n,unsigned int mode,unsigned int dim,double * coeff)1404 GMT_LOCAL void greenspline_do_normalization (struct GMTAPI_CTRL *API, double **X, double *obs, uint64_t n, unsigned int mode, unsigned int dim, double *coeff) {
1405 /* We always remove/restore the mean observation value. mode is a combination of bitflags that affects what we do:
1406 * Bit GREENSPLINE_TREND will also remove linear trend
1407 * Bit GREENSPLINE_NORM will normalize residuals by range
1408 */
1409
1410 uint64_t i;
1411 double d, min = DBL_MAX, max = -DBL_MAX;
1412 char *type[4] = {"Remove mean\n", "Normalization mode: Remove %d-D linear trend\n", "Remove mean and normalize data\n", "Normalization mode: Remove %d-D linear trend and normalize data\n"};
1413 if (mode % 2)
1414 GMT_Report (API, GMT_MSG_INFORMATION, type[mode], dim);
1415 else
1416 GMT_Report (API, GMT_MSG_INFORMATION, "Normalization mode: %s\n", type[mode]);
1417 if (dim == 1) { /* 1-D trend or mean only */
1418 greenspline_do_normalization_1d (API, X, obs, n, mode, coeff);
1419 return;
1420 }
1421 gmt_M_memset (coeff, GSP_LENGTH, double);
1422 for (i = 0; i < n; i++) { /* Find mean z-value */
1423 coeff[GSP_MEAN_Z] += obs[i];
1424 if ((mode & GREENSPLINE_TREND) == 0) continue; /* Else we also sum up x and y to get their means */
1425 coeff[GSP_MEAN_X] += X[i][GMT_X];
1426 coeff[GSP_MEAN_Y] += X[i][GMT_Y];
1427 }
1428 coeff[GSP_MEAN_Z] /= n; /* Average z value to remove/restore */
1429
1430 if (mode & GREENSPLINE_TREND) { /* Solve for LS plane using deviations from (0, 0, 0) */
1431 double xx, yy, zz, sxx, sxy, sxz, syy, syz;
1432 sxx = sxy = sxz = syy = syz = 0.0;
1433 coeff[GSP_MEAN_X] /= n; /* Mean x */
1434 coeff[GSP_MEAN_Y] /= n; /* Mean y */
1435 for (i = 0; i < n; i++) {
1436
1437 xx = X[i][GMT_X] - coeff[GSP_MEAN_X];
1438 yy = X[i][GMT_Y] - coeff[GSP_MEAN_Y];
1439 zz = obs[i] - coeff[GSP_MEAN_Z];
1440 /* xx,yy,zz are residuals relative to (0,0,0) */
1441 sxx += (xx * xx);
1442 sxz += (xx * zz);
1443 sxy += (xx * yy);
1444 syy += (yy * yy);
1445 syz += (yy * zz);
1446 }
1447
1448 d = sxx*syy - sxy*sxy;
1449 if (d != 0.0) {
1450 coeff[GSP_SLP_X] = (sxz*syy - sxy*syz)/d;
1451 coeff[GSP_SLP_Y] = (sxx*syz - sxy*sxz)/d;
1452 }
1453 }
1454
1455 /* Remove plane (or just mean) */
1456
1457 for (i = 0; i < n; i++) { /* Also find min/max or residuals in the process */
1458 obs[i] -= coeff[GSP_MEAN_Z]; /* Always remove mean data value */
1459 if (mode & GREENSPLINE_TREND) obs[i] -= (coeff[GSP_SLP_X] * (X[i][GMT_X] - coeff[GSP_MEAN_X]) + coeff[GSP_SLP_Y] * (X[i][GMT_Y] - coeff[GSP_MEAN_Y]));
1460 if (obs[i] < min) min = obs[i];
1461 if (obs[i] > max) max = obs[i];
1462 }
1463 if (mode & GREENSPLINE_NORM) { /* Normalize by range */
1464 coeff[GSP_RANGE] = MAX (fabs(min), fabs(max)); /* Determine range */
1465 d = (coeff[GSP_RANGE] == 0.0) ? 1.0 : 1.0 / coeff[GSP_RANGE];
1466 for (i = 0; i < n; i++) obs[i] *= d; /* Normalize 0-1 */
1467 }
1468
1469 /* Recover obs(x,y) = w_norm(x,y) * coeff[GSP_RANGE] + coeff[GSP_MEAN_Z] + coeff[GSP_SLP_X]*(x-coeff[GSP_MEAN_X]) + coeff[GSP_SLP_Y]*(y-coeff[GSP_MEAN_Y]) */
1470 GMT_Report (API, GMT_MSG_INFORMATION, "2-D Normalization coefficients: zoff = %g xslope = %g xmean = %g yslope = %g ymean = %g data range = %g\n",
1471 coeff[GSP_MEAN_Z], coeff[GSP_SLP_X], coeff[GSP_MEAN_X], coeff[GSP_SLP_Y], coeff[GSP_MEAN_Y], coeff[GSP_RANGE]);
1472 }
1473
greenspline_get_radius(struct GMT_CTRL * GMT,double * X0,double * X1,unsigned int dim)1474 GMT_LOCAL double greenspline_get_radius (struct GMT_CTRL *GMT, double *X0, double *X1, unsigned int dim) {
1475 double r = 0.0;
1476 /* Get distance between the two points */
1477 switch (dim) {
1478 case 1: /* 1-D, just get x difference */
1479 r = fabs (X0[GMT_X] - X1[GMT_X]);
1480 break;
1481 case 2: /* 2-D Cartesian or spherical surface in meters */
1482 r = gmt_distance (GMT, X0[GMT_X], X0[GMT_Y], X1[GMT_X], X1[GMT_Y]);
1483 break;
1484 case 3: /* 3-D Cartesian */
1485 r = hypot (X0[GMT_X] - X1[GMT_X], X0[GMT_Y] - X1[GMT_Y]);
1486 r = hypot (r, X0[GMT_Z] - X1[GMT_Z]);
1487 break;
1488 }
1489 return (r);
1490 }
1491
greenspline_get_dircosine(struct GMT_CTRL * GMT,double * D,double * X0,double * X1,unsigned int dim,bool baz)1492 GMT_LOCAL double greenspline_get_dircosine (struct GMT_CTRL *GMT, double *D, double *X0, double *X1, unsigned int dim, bool baz) {
1493 /* D, the directional cosines of the observed gradient:
1494 * X0: Observation point.
1495 * X1: Prediction point.
1496 * Compute N, the direction cosine of X1-X2, then C = D dot N.
1497 */
1498 int ii;
1499 double az, C = 0.0, N[3];
1500
1501 switch (dim) {
1502 case 1: /* 1-D, always 1 */
1503 C = 1.0;
1504 break;
1505 case 2: /* 2-D */
1506 az = gmt_az_backaz (GMT, X0[GMT_X], X0[GMT_Y], X1[GMT_X], X1[GMT_Y], baz);
1507 sincosd (az, &N[GMT_X], &N[GMT_Y]);
1508 for (ii = 0; ii < 2; ii++) C += D[ii] * N[ii]; /* Dot product of 2-D unit vectors */
1509 break;
1510 case 3: /* 3-D */
1511 for (ii = 0; ii < 3; ii++) N[ii] = X1[ii] - X0[ii]; /* Difference vector */
1512 gmt_normalize3v (GMT, N); /* Normalize to unit vector */
1513 C = gmt_dot3v (GMT, D, N); /* Dot product of 3-D unit vectors */
1514 if (baz) C = -C; /* The opposite direction for X0-X1 */
1515 break;
1516 }
1517 return (C);
1518 }
1519
greenspline_dump_system(double * A,double * b,uint64_t nm,char * string)1520 GMT_LOCAL void greenspline_dump_system (double *A, double *b, uint64_t nm, char *string) {
1521 /* Dump an A | b system to stderr for debugging */
1522 uint64_t row, col, ij;
1523 fprintf (stderr, "\n%s\n", string);
1524 for (row = 0; row < nm; row++) {
1525 ij = row * nm;
1526 fprintf (stderr, "%12.6f", A[ij++]);
1527 for (col = 1; col < nm; col++, ij++) fprintf (stderr, "\t%12.6f", A[ij]);
1528 fprintf (stderr, "\t||\t%12.6f\n", b[row]);
1529 }
1530 }
1531
greenspline_set_filename(char * name,unsigned int k,unsigned int width,unsigned int mode,char * file)1532 GMT_LOCAL void greenspline_set_filename (char *name, unsigned int k, unsigned int width, unsigned int mode, char *file) {
1533 /* Turn name, eigenvalue number k, precision width and mode into a filename, e.g.,
1534 * ("solution.grd", 33, 3, GMT_SVD_INCREMENTAL, file) will give solution_inc_033.grd */
1535 unsigned int s = strlen (name) - 1;
1536 static char *type[3] = {"", "inc", "cum"};
1537 while (name[s] != '.') s--; /* Wind backwards to start of extension */
1538 name[s] = '\0'; /* Temporarily chop off extension */
1539 sprintf (file, "%s_%s_%*.*d.%s", name, type[mode], width, width, k, &name[s+1]);
1540 name[s] = '.'; /* Restore original name */
1541 }
1542
1543 #define bailout(code) {gmt_M_free_options (mode); return (code);}
1544 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
1545
GMT_greenspline(void * V_API,int mode,void * args)1546 EXTERN_MSC int GMT_greenspline (void *V_API, int mode, void *args) {
1547 uint64_t col, row, n_read, p, k, i, j, seg, m, n, nm, n_ok = 0, ij, ji, ii, n_duplicates = 0, n_skip = 0;
1548 unsigned int dimension = 0, normalize = 0, n_cols, n_layers = 1, w_col, L_Max = 0;
1549 size_t n_alloc;
1550 int error = GMT_NOERROR, out_ID, way, n_columns, n_use;
1551 bool delete_grid = false, check_longitude, skip, write_3D_records = false;
1552
1553 char *method[N_METHODS] = {"Minimum curvature Cartesian spline [1-D]",
1554 "Minimum curvature Cartesian spline [2-D]",
1555 "Minimum curvature Cartesian spline [3-D]",
1556 "Continuous curvature Cartesian spline in tension [1-D]",
1557 "Continuous curvature Cartesian spline in tension [2-D]",
1558 "Continuous curvature Cartesian spline in tension [3-D]",
1559 "Regularized Cartesian spline in tension [2-D]",
1560 "Regularized Cartesian spline in tension [3-D]",
1561 "Minimum curvature spherical spline",
1562 "Continuous curvature spherical spline in tension",
1563 "Linear Cartesian spline [1-D]",
1564 "Bilinear Cartesian spline [2-D]"};
1565
1566 gmt_grdfloat *data = NULL;
1567
1568 double *v = NULL, *s = NULL, *b = NULL, *ssave = NULL;
1569 double *obs = NULL, **D = NULL, **X = NULL, *alpha = NULL, *in = NULL, *orig_obs = NULL;
1570 double mem, part, C, p_val, r, par[N_PARAMS], norm[GSP_LENGTH], az = 0, grad;
1571 double *A = NULL, *A_orig = NULL, r_min, r_max, err_sum = 0.0, var_sum = 0.0;
1572 #ifdef DEBUG
1573 double x0 = 0.0, x1 = 5.0;
1574 #endif
1575
1576 FILE *fp = NULL;
1577
1578 double (*G) (struct GMT_CTRL *, double, double *, struct GREENSPLINE_LOOKUP *) = NULL; /* Pointer to chosen Green's function */
1579 double (*dGdr) (struct GMT_CTRL *, double, double *, struct GREENSPLINE_LOOKUP *) = NULL; /* Pointer to chosen gradient of Green's function */
1580
1581 struct GMT_GRID *Grid = NULL, *Out = NULL;
1582 struct GMT_GRID_HEADER *header = NULL;
1583 struct GMT_CUBE *Cube = NULL; /* Structure to hold output cube if 3-D interpolation */
1584
1585 struct GREENSPLINE_LOOKUP *Lz = NULL, *Lg = NULL;
1586 struct GMT_DATATABLE *T = NULL;
1587 struct GMT_DATASET *Nin = NULL;
1588 struct GMT_GRID_INFO info;
1589 struct GMT_RECORD *In = NULL, *Rec = NULL;
1590 struct GREENSPLINE_CTRL *Ctrl = NULL;
1591 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
1592 struct GMT_OPTION *options = NULL;
1593 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1594
1595 /*----------------------- Standard module initialization and parsing ----------------------*/
1596
1597 if (API == NULL) return (GMT_NOT_A_SESSION);
1598 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
1599 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
1600
1601 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
1602
1603 /* Parse the command-line arguments */
1604
1605 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 */
1606 dimension = greenspline_pre_parser (GMT, options); /* Check -Z and possibly change default to geographic data mode before -R is parsed */
1607 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
1608 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
1609 Ctrl->dimension = dimension;
1610 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
1611
1612 /*---------------------------- This is the greenspline main code ----------------------------*/
1613
1614 gmt_enable_threads (GMT); /* Set number of active threads, if supported */
1615 dimension = (Ctrl->Z.mode == 0) ? 1 : ((Ctrl->Z.mode == 5) ? 3 : 2);
1616 gmt_M_memset (par, N_PARAMS, double);
1617 gmt_M_memset (norm, GSP_LENGTH, double);
1618 gmt_M_memset (&info, 1, struct GMT_GRID_INFO);
1619
1620 write_3D_records = (dimension == 3 && !Ctrl->G.active); /* Just so it is only true if 3-D and no output filename was given */
1621
1622 /* As many S.mode reflects 1-D after parse, here we increment S.mode if dimension is 2 or 3 */
1623 if (Ctrl->S.mode == SANDWELL_1987_1D || Ctrl->S.mode == WESSEL_BERCOVICI_1998_1D) Ctrl->S.mode += (dimension - 1);
1624 if (Ctrl->S.mode == LINEAR_1D) Ctrl->S.mode += (dimension - 1);
1625 if (Ctrl->S.mode == MITASOVA_MITAS_1993_2D) Ctrl->S.mode += (dimension - 2);
1626
1627 way = gmt_M_is_spherical (GMT) ? GMT_GREATCIRCLE : GMT_GEODESIC;
1628 Ctrl->Z.mode--; /* Since I added 0 to be 1-D later so now it is -1 */
1629 switch (Ctrl->Z.mode) { /* Set pointers to 2-D distance functions */
1630 case -1: /* Cartesian 1-D x data */
1631 normalize = GREENSPLINE_TREND + GREENSPLINE_NORM;
1632 break;
1633 case 0: /* Cartesian 2-D x,y data */
1634 error = gmt_init_distaz (GMT, 'X', 0, GMT_MAP_DIST);
1635 normalize = GREENSPLINE_TREND + GREENSPLINE_NORM;
1636 break;
1637 case 1: /* 2-D lon, lat data, but scale to Cartesian flat earth km */
1638 gmt_set_geographic (GMT, GMT_IN);
1639 gmt_set_geographic (GMT, GMT_OUT);
1640 error = gmt_init_distaz (GMT, 'k', GMT_FLATEARTH, GMT_MAP_DIST);
1641 normalize = GREENSPLINE_TREND + GREENSPLINE_NORM;
1642 break;
1643 case 2: /* 2-D lon, lat data, use spherical distances in km (geodesic if PROJ_ELLIPSOID is nor sphere) */
1644 gmt_set_geographic (GMT, GMT_IN);
1645 gmt_set_geographic (GMT, GMT_OUT);
1646 error = gmt_init_distaz (GMT, 'k', way, GMT_MAP_DIST);
1647 normalize = GREENSPLINE_NORM;
1648 break;
1649 case 3: /* 2-D lon, lat data, and Green's function needs cosine of spherical or geodesic distance */
1650 gmt_set_geographic (GMT, GMT_IN);
1651 gmt_set_geographic (GMT, GMT_OUT);
1652 error = gmt_init_distaz (GMT, 'S', way, GMT_MAP_DIST);
1653 normalize = GREENSPLINE_NORM;
1654 break;
1655 case 4: /* 3-D Cartesian x,y,z data handled separately */
1656 normalize = GREENSPLINE_NORM;
1657 break;
1658 default: /* Cannot happen unless we make a bug */
1659 GMT_Report (API, GMT_MSG_ERROR, "BUG since D (=%d) cannot be outside 0-5 range\n", Ctrl->Z.mode+1);
1660 break;
1661 }
1662 if (error == GMT_NOT_A_VALID_TYPE) Return (error);
1663
1664 if (Ctrl->Z.mode <= 1 && Ctrl->L.active)
1665 normalize = GREENSPLINE_NORM; /* Do not de-plane, just remove mean and normalize */
1666 else if (Ctrl->Z.mode > 1 && Ctrl->L.active)
1667 GMT_Report (API, GMT_MSG_ERROR, "-L ignored for -Z modes 3-5\n");
1668
1669 if (Ctrl->Q.active && dimension == 2) sincosd (Ctrl->Q.az, &Ctrl->Q.dir[GMT_X], &Ctrl->Q.dir[GMT_Y]);
1670
1671 /* Now we are ready to take on some input values */
1672
1673 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
1674 Return (API->error);
1675 }
1676 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
1677 Return (API->error);
1678 }
1679
1680 n_cols = (Ctrl->W.active) ? dimension + 1 : dimension; /* So X[k][dimension] holds the weight if -W is active */
1681 if ((error = GMT_Set_Columns (API, GMT_IN, n_cols+1, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR)
1682 Return (error);
1683 n_alloc = GMT_INITIAL_MEM_ROW_ALLOC;
1684 X = gmt_M_memory (GMT, NULL, n_alloc, double *);
1685 obs = gmt_M_memory (GMT, NULL, n_alloc, double);
1686 check_longitude = (dimension == 2 && (Ctrl->Z.mode == 1 || Ctrl->Z.mode == 2));
1687 w_col = dimension + 1; /* Where weights would be in input, if given */
1688
1689 GMT_Report (API, GMT_MSG_INFORMATION, "Read input data and check for data constraint duplicates\n");
1690 n = m = n_read = 0;
1691 r_min = DBL_MAX; r_max = -DBL_MAX;
1692 do { /* Keep returning records until we reach EOF */
1693 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
1694 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
1695 for (p = 0; p < n; p++) gmt_M_free (GMT, X[p]);
1696 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1697 Return (GMT_RUNTIME_ERROR);
1698 }
1699 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
1700 break;
1701 continue; /* Go back and read the next record */
1702 }
1703 if (In->data == NULL) {
1704 gmt_quit_bad_record (API, In);
1705 for (p = 0; p < n; p++) gmt_M_free (GMT, X[p]);
1706 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1707 Return (API->error);
1708 }
1709
1710 /* Data record to process */
1711 in = In->data; /* Only need to process numerical part here */
1712
1713 if (check_longitude) {
1714 /* Ensure geographic longitudes fit the range since the normalization function expects it */
1715 if (in[GMT_X] < Ctrl->R3.range[XLO] && (in[GMT_X] + 360.0) < Ctrl->R3.range[XHI]) in[GMT_X] += 360.0;
1716 else if (in[GMT_X] > Ctrl->R3.range[XHI] && (in[GMT_X] - 360.0) > Ctrl->R3.range[XLO]) in[GMT_X] -= 360.0;
1717 }
1718
1719 if (X[n] == NULL) X[n] = gmt_M_memory (GMT, NULL, n_cols, double); /* Allocate space for this constraint */
1720 for (k = 0; k < dimension; k++) X[n][k] = in[k]; /* Get coordinates + optional weights (if -W) */
1721 /* Check for duplicates */
1722 skip = false;
1723 for (i = 0; !skip && i < n; i++) {
1724 r = greenspline_get_radius (GMT, X[i], X[n], dimension);
1725 if (gmt_M_is_zero (r)) { /* Duplicates will give zero point separation */
1726 if (doubleAlmostEqualZero (in[dimension], obs[i])) {
1727 GMT_Report (API, GMT_MSG_WARNING,
1728 "Data constraint %" PRIu64 " is identical to %" PRIu64 " and will be skipped\n", n_read, i);
1729 skip = true;
1730 n_skip++;
1731 }
1732 else {
1733 GMT_Report (API, GMT_MSG_ERROR,
1734 "Data constraint %" PRIu64 " and %" PRIu64 " occupy the same location but differ"
1735 " in observation (%.12g vs %.12g)\n", n_read, i, in[dimension], obs[i]);
1736 n_duplicates++;
1737 }
1738 }
1739 else {
1740 if (r < r_min) r_min = r;
1741 if (r > r_max) r_max = r;
1742 }
1743 }
1744 n_read++;
1745 if (skip) continue; /* Current point was a duplicate of a previous point */
1746
1747 if (Ctrl->W.active) { /* Planning a weighted solution */
1748 if (Ctrl->W.mode == GSP_GOT_SIG) { /* Got sigma, must convert to weight */
1749 err_sum += in[w_col] * in[w_col]; /* Sum up variance first */
1750 X[n][dimension] = 1.0 / in[w_col]; /* We will square later */
1751 }
1752 else /* Got weight, use as is, no squaring later */
1753 X[n][dimension] = in[w_col];
1754 }
1755 var_sum += in[dimension] * in[dimension]; /* Sum up data variance */
1756 obs[n++] = in[dimension];
1757
1758 if (n == n_alloc) { /* Get more memory */
1759 n_alloc <<= 1;
1760 X = gmt_M_memory (GMT, X, n_alloc, double *);
1761 obs = gmt_M_memory (GMT, obs, n_alloc, double);
1762 }
1763 } while (true);
1764
1765 if (n_skip && n < n_alloc && X[n]) /* If we end with a skip then we have allocated one too many */
1766 gmt_M_free (GMT, X[n]);
1767
1768 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
1769 for (p = 0; p < n; p++) gmt_M_free (GMT, X[p]);
1770 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1771 Return (API->error);
1772 }
1773
1774 X = gmt_M_memory (GMT, X, n, double *);
1775 obs = gmt_M_memory (GMT, obs, n, double);
1776 nm = n;
1777 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " unique data constraints\n", n);
1778 if (n_skip) GMT_Report (API, GMT_MSG_WARNING, "Skipped %" PRIu64 " data constraints as duplicates\n", n_skip);
1779
1780 if (Ctrl->W.active && Ctrl->W.mode == GSP_GOT_SIG) { /* Got data uncertainties */
1781 err_sum = sqrt (err_sum / nm); /* Mean data uncertainty */
1782 GMT_Report (API, GMT_MSG_INFORMATION, "Mean data uncertainty is %g\n", err_sum);
1783 }
1784
1785 if (Ctrl->A.active) { /* Read gradient constraints from file */
1786 unsigned int n_A_cols = 0;
1787 struct GMT_DATASET *Din = NULL;
1788 struct GMT_DATASEGMENT *Slp = NULL;
1789 switch (dimension) {
1790 case 1: /* 1-D */
1791 switch (Ctrl->A.mode) {
1792 case 0: n_A_cols = 2; break;/* x, slope */
1793 default:
1794 GMT_Report (API, GMT_MSG_ERROR, "Bad gradient mode selected for 1-D data (%d) - aborting!\n", Ctrl->A.mode);
1795 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1796 Return (GMT_DATA_READ_ERROR);
1797 break;
1798 }
1799 break;
1800 case 2: /* 2-D */
1801 switch (Ctrl->A.mode) {
1802 case 1: n_A_cols = 4; break; /* (x, y, az, gradient) */
1803 case 2: n_A_cols = 4; break; /* (x, y, gradient, azimuth) */
1804 case 3: n_A_cols = 4; break; /* (x, y, direction, gradient) */
1805 case 4: n_A_cols = 4; break; /* (x, y, gx, gy) */
1806 case 5: n_A_cols = 5; break; /* (x, y, nx, ny, gradient) */
1807 default:
1808 GMT_Report (API, GMT_MSG_ERROR, "Bad gradient mode selected for 2-D data (%d) - aborting!\n", Ctrl->A.mode);
1809 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1810 Return (GMT_DATA_READ_ERROR);
1811 break;
1812 }
1813 break;
1814 case 3: /* 3-D */
1815 switch (Ctrl->A.mode) {
1816 case 4: n_A_cols = 6; break; /* (x, y, z, gx, gy, gz) */
1817 case 5: n_A_cols = 7; break; /* (x, y, z, nx, ny, nz, gradient) */
1818 default:
1819 GMT_Report (API, GMT_MSG_ERROR, "Bad gradient mode selected for 3-D data (%d) - aborting!\n", Ctrl->A.mode);
1820 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1821 Return (GMT_DATA_READ_ERROR);
1822 break;
1823 }
1824 break;
1825 }
1826
1827 if (GMT->common.b.active[GMT_IN]) GMT->common.b.ncol[GMT_IN]++; /* Must assume it is just one extra column */
1828 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -C,-F,-L files */
1829 if ((Din = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->A.file, NULL)) == NULL) {
1830 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
1831 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1832 Return (API->error);
1833 }
1834 if (Din->n_columns < n_A_cols) {
1835 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least %u are needed\n", (int)Din->n_columns, n_A_cols);
1836 Return (GMT_DIM_TOO_SMALL);
1837 }
1838 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
1839 m = Din->n_records; /* Total number of gradient constraints */
1840 nm += m; /* New total of linear equations to solve */
1841 X = gmt_M_memory (GMT, X, nm, double *);
1842 for (k = n; k < nm; k++) X[k] = gmt_M_memory (GMT, NULL, n_cols, double);
1843 obs = gmt_M_memory (GMT, obs, nm, double);
1844 D = gmt_M_memory (GMT, NULL, m, double *);
1845 for (k = 0; k < m; k++) D[k] = gmt_M_memory (GMT, NULL, n_cols, double);
1846 n_skip = n_read = 0;
1847 for (seg = k = 0, p = n; seg < Din->n_segments; seg++) {
1848 Slp = Din->table[0]->segment[seg];
1849 for (row = 0; row < Slp->n_rows; row++, k++, p++) {
1850 for (ii = 0; ii < n_cols; ii++) X[p][ii] = Slp->data[ii][row];
1851 switch (dimension) {
1852 case 1: /* 1-D: x, slope */
1853 D[k][0] = 1.0; /* Dummy since there is no direction for 1-D spline (the gradient is in the x-y plane) */
1854 obs[p] = Slp->data[dimension][row];
1855 break;
1856 case 2: /* 2-D */
1857 switch (Ctrl->A.mode) {
1858 case 1: /* (x, y, az, gradient) */
1859 az = D2R * Slp->data[2][row];
1860 obs[p] = Slp->data[3][row];
1861 break;
1862 case 2: /* (x, y, gradient, azimuth) */
1863 az = D2R * Slp->data[3][row];
1864 obs[p] = Slp->data[2][row];
1865 break;
1866 case 3: /* (x, y, direction, gradient) */
1867 az = M_PI_2 - D2R * Slp->data[2][row];
1868 obs[p] = Slp->data[3][row];
1869 break;
1870 case 4: /* (x, y, gx, gy) */
1871 az = atan2 (Slp->data[2][row], Slp->data[3][row]); /* Get azimuth of gradient */
1872 obs[p] = hypot (Slp->data[3][row], Slp->data[3][row]); /* Get magnitude of gradient */
1873 break;
1874 case 5: /* (x, y, nx, ny, gradient) */
1875 az = atan2 (Slp->data[2][row], Slp->data[3][row]); /* Get azimuth of gradient */
1876 obs[p] = Slp->data[4][row]; /* Magnitude of gradient */
1877 break;
1878 }
1879 sincos (az, &D[k][GMT_X], &D[k][GMT_Y]);
1880 break;
1881 case 3: /* 3-D */
1882 switch (Ctrl->A.mode) {
1883 case 4: /* (x, y, z, gx, gy, gz) */
1884 for (ii = 0; ii < 3; ii++) D[k][ii] = Slp->data[3+ii][row]; /* Get the gradient vector */
1885 obs[p] = gmt_mag3v (GMT, D[k]); /* This is the gradient magnitude */
1886 gmt_normalize3v (GMT, D[k]); /* These are the direction cosines of the gradient */
1887 break;
1888 case 5: /* (x, y, z, nx, ny, nz, gradient) */
1889 for (ii = 0; ii < 3; ii++) D[k][ii] = Slp->data[3+ii][row]; /* Get the unit vector */
1890 obs[p] = Slp->data[6][row]; /* This is the gradient magnitude */
1891 break;
1892 }
1893 break;
1894 default:
1895 GMT_Report (API, GMT_MSG_ERROR, "Bad dimension selected (%d) - aborting!\n", dimension);
1896 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
1897 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1898 Return (GMT_DATA_READ_ERROR);
1899 break;
1900 }
1901 /* Check for duplicates */
1902 skip = false;
1903 for (i = n; !skip && i < p; i++) {
1904 r = greenspline_get_radius (GMT, X[i], X[p], dimension);
1905 if (gmt_M_is_zero (r)) { /* Duplicates will give zero point separation */
1906 if (doubleAlmostEqualZero (in[dimension], obs[i])) {
1907 GMT_Report (API, GMT_MSG_WARNING, "Slope constraint %" PRIu64 " is identical to %" PRIu64
1908 " and will be skipped\n", n_read, i-n);
1909 skip = true;
1910 n_skip++;
1911 }
1912 else {
1913 GMT_Report (API, GMT_MSG_ERROR, "Slope constraint %" PRIu64 " and %" PRIu64
1914 " occupy the same location but differ in observation (%.12g vs %.12g)\n", n_read, i-n, obs[p], obs[i]);
1915 n_duplicates++;
1916 }
1917 }
1918 else {
1919 if (r < r_min) r_min = r;
1920 if (r > r_max) r_max = r;
1921 }
1922 }
1923 n_read++;
1924 if (skip) p--; /* Current point was a duplicate of a previous point; reduce p since it will be incremented in the loop */
1925 }
1926 }
1927 if (GMT_Destroy_Data (API, &Din) != GMT_NOERROR) {
1928 Return (API->error);
1929 }
1930 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " unique slope constraints\n", m);
1931 if (n_skip) GMT_Report (API, GMT_MSG_WARNING, "Skipped %" PRIu64 " slope constraints as duplicates\n", n_skip);
1932 }
1933
1934 /* Check for duplicates which would result in a singular matrix system; also update min/max radius */
1935
1936 GMT_Report (API, GMT_MSG_INFORMATION, "Distance between the closest constraints: %.12g]\n", r_min);
1937 GMT_Report (API, GMT_MSG_INFORMATION, "Distance between most distant constraints: %.12g]\n", r_max);
1938
1939 if (n_duplicates) { /* These differ in observation value so need to be averaged, medianed, or whatever first */
1940 if (!Ctrl->C.active || gmt_M_is_zero (Ctrl->C.value)) {
1941 GMT_Report (API, GMT_MSG_ERROR,
1942 "Found %" PRIu64 " data constraint duplicates with different observation values\n", n_duplicates);
1943 GMT_Report (API, GMT_MSG_ERROR,
1944 "You must reconcile duplicates before running greenspline since they will result in a singular matrix\n");
1945 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
1946 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1947 if (m) {
1948 for (p = 0; p < m; p++) gmt_M_free (GMT, D[p]);
1949 gmt_M_free (GMT, D);
1950 }
1951 Return (GMT_DATA_READ_ERROR);
1952 }
1953 else {
1954 GMT_Report (API, GMT_MSG_WARNING,
1955 "Found %" PRIu64 " data constraint duplicates with different observation values\n", n_duplicates);
1956 GMT_Report (API, GMT_MSG_WARNING, "Expect some eigenvalues to be identically zero\n");
1957 }
1958 }
1959
1960 if (m > 0 && (normalize & GREENSPLINE_TREND)) {
1961 normalize = GREENSPLINE_NORM; /* Only allow taking out data mean for mixed z/slope data */
1962 GMT_Report (API, GMT_MSG_WARNING, "Can only remove/restore mean z in mixed {z, grad(z)} data sets\n");
1963 }
1964
1965 if (m == 0)
1966 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " data points, yielding a %" PRIu64 " by %" PRIu64 " set of linear equations\n",
1967 n, nm, nm);
1968 else
1969 GMT_Report (API, GMT_MSG_INFORMATION, "Found %" PRIu64 " data points and %" PRIu64 " gradients, yielding a %" PRIu64 " by %"
1970 PRIu64 " set of linear equations\n", n, m, nm, nm);
1971
1972 if (Ctrl->T.file) { /* Existing grid that will have zeros and NaNs, only */
1973 if ((Grid = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->T.file, NULL)) == NULL) { /* Get header only */
1974 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1975 Return (API->error);
1976 }
1977 if (!(Grid->header->wesn[XLO] == Ctrl->R3.range[0] && Grid->header->wesn[XHI] == Ctrl->R3.range[1] &&
1978 Grid->header->wesn[YLO] == Ctrl->R3.range[2] && Grid->header->wesn[YHI] == Ctrl->R3.range[3])) {
1979 GMT_Report (API, GMT_MSG_ERROR, "The mask grid does not match your specified region\n");
1980 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1981 Return (GMT_RUNTIME_ERROR);
1982 }
1983 if (! (Grid->header->inc[GMT_X] == Ctrl->I.inc[GMT_X] && Grid->header->inc[GMT_Y] == Ctrl->I.inc[GMT_Y])) {
1984 GMT_Report (API, GMT_MSG_ERROR, "The mask grid resolution does not match your specified grid spacing\n");
1985 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1986 Return (GMT_RUNTIME_ERROR);
1987 }
1988 if (! (Grid->header->registration == GMT->common.R.registration)) {
1989 GMT_Report (API, GMT_MSG_ERROR, "The mask grid registration does not match your specified grid registration\n");
1990 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1991 Return (GMT_RUNTIME_ERROR);
1992 }
1993 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->T.file, Grid) == NULL) { /* Get data */
1994 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
1995 Return (API->error);
1996 }
1997 (void)gmt_set_outgrid (GMT, Ctrl->T.file, false, 0, Grid, &Out); /* true if input is a read-only array; otherwise Out is just a pointer to Grid */
1998 n_ok = Grid->header->nm;
1999 gmt_M_grd_loop (GMT, Grid, row, col, ij) if (gmt_M_is_fnan (Grid->data[ij])) n_ok--;
2000 }
2001 else if (Ctrl->N.active) { /* Read output locations from file */
2002 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -C,-F,-L files */
2003 if ((Nin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) {
2004 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
2005 Return (API->error);
2006 }
2007 if (Nin->n_columns < dimension) {
2008 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least %u are needed\n",
2009 Ctrl->N.file, (int)Nin->n_columns, dimension);
2010 gmt_M_free (GMT, X); gmt_M_free (GMT, obs);
2011 Return (GMT_DIM_TOO_SMALL);
2012 }
2013 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
2014 T = Nin->table[0];
2015 }
2016 else { /* Fill in an equidistant output table, grid, or cube */
2017 if (dimension == 1) { /* Dummy grid to hold the 1-D info */
2018 if ((Grid = gmt_create_grid (GMT)) == NULL) Return (API->error);
2019 delete_grid = true;
2020 Grid->header->wesn[XLO] = Ctrl->R3.range[XLO]; Grid->header->wesn[XHI] = Ctrl->R3.range[XHI];
2021 Grid->header->registration = GMT->common.R.registration;
2022 Grid->header->inc[GMT_X] = Ctrl->I.inc[GMT_X];
2023 Grid->header->n_rows = 1; /* So that output logic will work for 1-D which only has columns */
2024 n_ok = Grid->header->n_columns = gmt_M_grd_get_nx (GMT, Grid->header);
2025 header = Grid->header;
2026 data = gmt_M_memory (GMT, NULL, Grid->header->n_columns, gmt_grdfloat);
2027 }
2028 else if (dimension == 2) { /* Need a full-fledged Grid creation since we are writing it to who knows where */
2029 if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->R3.range, Ctrl->I.inc, \
2030 GMT->common.R.registration, GMT_NOTSET, NULL)) == NULL) Return (API->error);
2031 n_ok = Grid->header->nm;
2032 header = Grid->header;
2033 data = Grid->data; /* Pointer to the float 2-D grid */
2034 if (Ctrl->D.active && gmt_decode_grd_h_info (GMT, Ctrl->D.information, Grid->header)) {
2035 Return (GMT_PARSE_ERROR);
2036 }
2037 }
2038 else { /* 3-D cube needed */
2039 if ((Cube = GMT_Create_Data (API, GMT_IS_CUBE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, NULL, Ctrl->R3.range, Ctrl->I.inc, \
2040 GMT->common.R.registration, GMT_NOTSET, NULL)) == NULL) Return (API->error);
2041 n_layers = Cube->header->n_bands;
2042 n_ok = Cube->header->nm * n_layers;
2043 header = Cube->header;
2044 data = Cube->data; /* Pointer to the float 3-D cube */
2045 if (Ctrl->D.active && gmt_decode_cube_h_info (GMT, Ctrl->D.information, Cube)) {
2046 Return (GMT_PARSE_ERROR);
2047 }
2048 }
2049 Out = Grid; /* Just pointer since we created Grid above (except for cube) */
2050 }
2051
2052 switch (Ctrl->S.mode) { /* Assign pointers to Green's functions and the gradient and set up required parameters */
2053 case LINEAR_1D:
2054 case LINEAR_2D:
2055 G = &greenspline_spline1d_linear;
2056 dGdr = &greenspline_grad_spline1d_linear;
2057 break;
2058 case SANDWELL_1987_1D:
2059 G = &greenspline_spline1d_sandwell;
2060 dGdr = &greenspline_grad_spline1d_sandwell;
2061 break;
2062 case SANDWELL_1987_2D:
2063 G = &greenspline_spline2d_sandwell;
2064 dGdr = &greenspline_grad_spline2d_sandwell;
2065 break;
2066 case SANDWELL_1987_3D:
2067 G = &greenspline_spline3d_sandwell;
2068 dGdr = &greenspline_grad_spline3d_sandwell;
2069 break;
2070 case WESSEL_BERCOVICI_1998_1D:
2071 if (Ctrl->S.value[1] == 0.0 && Grid->header->inc[GMT_X] > 0.0) Ctrl->S.value[1] = Grid->header->inc[GMT_X];
2072 if (Ctrl->S.value[1] == 0.0) Ctrl->S.value[1] = 1.0;
2073 par[0] = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])) / Ctrl->S.value[1];
2074 par[1] = 2.0 / par[0];
2075 G = &greenspline_spline1d_Wessel_Bercovici;
2076 dGdr = &greenspline_grad_spline1d_Wessel_Bercovici;
2077 break;
2078 case WESSEL_BERCOVICI_1998_2D:
2079 if (Ctrl->S.value[1] == 0.0 && Grid->header->inc[GMT_X] > 0.0)
2080 Ctrl->S.value[1] = 0.5 * (Grid->header->inc[GMT_X] + Grid->header->inc[GMT_Y]);
2081 if (Ctrl->S.value[1] == 0.0) Ctrl->S.value[1] = 1.0;
2082 par[0] = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])) / Ctrl->S.value[1];
2083 par[1] = 2.0 / par[0];
2084 G = &greenspline_spline2d_Wessel_Bercovici;
2085 dGdr = &greenspline_grad_spline2d_Wessel_Bercovici;
2086 break;
2087 case WESSEL_BERCOVICI_1998_3D:
2088 if (Ctrl->S.value[1] == 0.0 && Grid->header->inc[GMT_X] > 0.0)
2089 Ctrl->S.value[1] = (Grid->header->inc[GMT_X] + Grid->header->inc[GMT_Y] + Cube->z_inc) / 3.0;
2090 if (Ctrl->S.value[1] == 0.0) Ctrl->S.value[1] = 1.0;
2091 par[0] = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])) / Ctrl->S.value[1];
2092 par[1] = 2.0 / par[0];
2093 G = &greenspline_spline3d_Wessel_Bercovici;
2094 dGdr = &greenspline_grad_spline3d_Wessel_Bercovici;
2095 break;
2096 case MITASOVA_MITAS_1993_2D:
2097 /* par[0] = Ctrl->S.value[0]; */
2098 if (Ctrl->S.value[1] == 0.0) Ctrl->S.value[1] = 1.0;
2099 p_val = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])) / Ctrl->S.value[1];
2100 GMT_Report (API, GMT_MSG_DEBUG, "p_val = %g\n", p_val);
2101 par[0] = p_val;
2102 par[1] = 0.25 * par[0] * par[0];
2103 G = &greenspline_spline2d_Mitasova_Mitas;
2104 dGdr = &greenspline_grad_spline2d_Mitasova_Mitas;
2105 break;
2106 case MITASOVA_MITAS_1993_3D:
2107 /* par[0] = Ctrl->S.value[0]; */
2108 if (Ctrl->S.value[1] == 0.0) Ctrl->S.value[1] = 1.0;
2109 p_val = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])) / Ctrl->S.value[1];
2110 GMT_Report (API, GMT_MSG_DEBUG, "p_val = %g\n", p_val);
2111 par[0] = p_val;
2112 par[1] = 0.25 * par[0] * par[0];
2113 G = &greenspline_spline3d_Mitasova_Mitas;
2114 dGdr = &greenspline_grad_spline3d_Mitasova_Mitas;
2115 break;
2116 case PARKER_1994:
2117 par[0] = 6.0 / (M_PI*M_PI);
2118 G = &greenspline_spline2d_Parker;
2119 dGdr = &greenspline_grad_spline2d_Parker;
2120 #ifdef DEBUG
2121 if (TEST) x0 = -1.0, x1 = 1.0;
2122 #endif
2123 break;
2124 case WESSEL_BECKER_2008:
2125 par[0] = sqrt (Ctrl->S.value[0] / (1.0 - Ctrl->S.value[0])); /* The p value */
2126 par[1] = Ctrl->S.value[2]; /* The truncation error */
2127 par[2] = -log (2.0) + (par[0]*par[0] - 1.0) / (par[0]*par[0]); /* Precalculate the constant for the l = 0 term here */
2128 Lz = gmt_M_memory (GMT, NULL, 1, struct GREENSPLINE_LOOKUP);
2129 #ifdef DEBUG
2130 if (TEST) Lg = gmt_M_memory (GMT, NULL, 1, struct GREENSPLINE_LOOKUP);
2131 else
2132 #endif
2133 if (Ctrl->A.active) Lg = gmt_M_memory (GMT, NULL, 1, struct GREENSPLINE_LOOKUP);
2134 L_Max = greenspline_get_max_L (GMT, par[0], par[1]);
2135 GMT_Report (API, GMT_MSG_DEBUG, "New scheme p = %g, err = %g, L_Max = %u\n", par[0], par[1], L_Max);
2136 greenspline_series_prepare (GMT, par[0], L_Max, Lz, Lg);
2137 /* Set up the cubic spline lookup/interpolation */
2138 par[7] = Ctrl->S.value[3];
2139 n_columns = irint (par[7]);
2140 par[8] = (Ctrl->S.rval[1] - Ctrl->S.rval[0]) / (par[7] - 1.0);
2141 par[9] = 1.0 / par[8];
2142 par[10] = Ctrl->S.rval[0];
2143 par[4] = par[8] * par[8]; /* Spline spacing squared, needed by greenspline_csplint */
2144
2145 GMT_Report (API, GMT_MSG_INFORMATION, "Precalculate -Sq lookup table with %d items from %g to %g\n", n_columns, Ctrl->S.rval[0], Ctrl->S.rval[1]);
2146 greenspline_spline2d_Wessel_Becker_init (GMT, par, Lz, Lg);
2147 G = &greenspline_spline2d_Wessel_Becker_lookup;
2148 dGdr = &greenspline_grad_spline2d_Wessel_Becker_lookup;
2149 #ifdef DEBUG
2150 if (TEST) x0 = -1.0, x1 = 1.0;
2151 #endif
2152 break;
2153 }
2154
2155 #ifdef DEBUG
2156 if (TEST) {
2157 GMT_Report (API, GMT_MSG_WARNING, "greenspline running in TEST mode for %s\n", method[Ctrl->S.mode]);
2158 printf ("# %s\n#x\tG\tdG/dx\tt\n", method[Ctrl->S.mode]);
2159 greenspline_dump_green (GMT, G, dGdr, par, x0, x1, 10001, Lz, Lg);
2160 gmt_free_grid (GMT, &Grid, dimension == 2);
2161 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
2162 greenspline_free_lookup (GMT, &Lz, 0);
2163 greenspline_free_lookup (GMT, &Lg, 1);
2164 Return (0);
2165 }
2166 #endif
2167
2168 if (dimension == 1) gmt_increase_abstime_format_precision (GMT, GMT_X, Ctrl->I.inc[GMT_X]); /* In case we need more sub-second precision output */
2169
2170 /* Remove mean (or LS plane) from data (we will add it back later) */
2171
2172 if (Ctrl->E.active) { /* Need to duplicate the data since SVD destroys it. */
2173 orig_obs = gmt_M_memory (GMT, NULL, nm, double);
2174 gmt_M_memcpy (orig_obs, obs, nm, double);
2175 }
2176
2177 greenspline_do_normalization (API, X, obs, n, normalize, dimension, norm);
2178
2179 /* Set up linear system Ax = obs. To clarify, the matrix A will be
2180 * of size nm by nm, where nm = n + m. Again, n is the number of
2181 * value constraints and m is the number of gradient constraints.
2182 * for most problems m will be 0.
2183 * The loops below takes advantage of the fact that A will be symmetrical
2184 * (except for terms involving gradients where A_ij = -A_ji). So we
2185 * start the loop over columns as col = row and deal with A)ij and A_ji
2186 * at the same time since we can evaluate the same costly G() function
2187 * [or dGdr () function)] once.
2188 */
2189
2190 mem = (double)nm * (double)nm * (double)sizeof (double); /* In bytes */
2191 GMT_Report (API, GMT_MSG_INFORMATION, "Square matrix A (size %d x %d) requires %s\n", (int)nm, (int)nm, gmt_memory_use ((size_t)mem, 1));
2192 A = gmt_M_memory (GMT, NULL, nm * nm, double);
2193
2194 GMT_Report (API, GMT_MSG_INFORMATION, "Build square linear system Ax = b using %s\n", method[Ctrl->S.mode]);
2195
2196 for (row = 0; row < nm; row++) { /* For each value or slope constraint */
2197 for (col = row; col < nm; col++) {
2198 ij = row * nm + col;
2199 ji = col * nm + row;
2200 r = greenspline_get_radius (GMT, X[col], X[row], dimension);
2201 if (row < n) { /* Value constraint (so entire row uses G) */
2202 A[ij] = G (GMT, r, par, Lz);
2203 if (ij == ji) /* Do the diagonal terms only once */
2204 continue;
2205 if (col < n)
2206 A[ji] = A[ij];
2207 else {
2208 /* Get D, the directional cosine between the two points */
2209 /* Then get C = gmt_dot3v (GMT, D, dataD); */
2210 /* A[ji] = dGdr (r, par, Lg) * C; */
2211 C = greenspline_get_dircosine (GMT, D[col-n], X[col], X[row], dimension, false);
2212 grad = dGdr (GMT, r, par, Lg);
2213 A[ji] = grad * C;
2214 }
2215 }
2216 else if (col > n) { /* Remaining gradient constraints (entire row uses dGdr) */
2217 if (ij == ji) continue; /* Diagonal gradient term from a point to itself is zero */
2218 C = greenspline_get_dircosine (GMT, D[row-n], X[col], X[row], dimension, true);
2219 grad = dGdr (GMT, r, par, Lg);
2220 A[ij] = grad * C;
2221 C = greenspline_get_dircosine (GMT, D[col-n], X[col], X[row], dimension, false);
2222 A[ji] = grad * C;
2223 }
2224 }
2225 }
2226
2227 if (Ctrl->debug.active) greenspline_dump_system (A, obs, nm, "A Matrix row || obs"); /* Dump the A | b system under debug */
2228 if (Ctrl->E.active && Ctrl->C.history == GMT_SVD_NO_HISTORY) { /* Needed A to evaluate misfit later as predict = A_orig * x */
2229 A_orig = gmt_M_memory (GMT, NULL, nm * nm, double);
2230 gmt_M_memcpy (A_orig, A, nm * nm, double);
2231 }
2232
2233 if (Ctrl->W.active) {
2234 /* Here we have requested an approximate fit instead of an exact interpolation.
2235 * For exact interpolation the weights do not matter, but here they do. Since
2236 * we wish to solve via SVD we must convert our unweighted A*x = b linear system
2237 * to the weighted W*A*x = W*b whose normal equations are [A'*S*A]*x = A'*S*b,
2238 * where S = W*W, the squared weights. This is N*x = r, where N = A'*S*A and
2239 * r = A'*S*b. Thus, we do these multiplication and store N and r in the
2240 * original A and obs vectors so that the continuation of the code can work as is.
2241 * Weighted solution idea credit: Leo Uieda. Jan 14, 2018.
2242 */
2243
2244 double *At = NULL, *AtS = NULL, *S = NULL; /* Need temporary work space */
2245
2246 GMT_Report (API, GMT_MSG_INFORMATION, "Forming weighted normal equations A'SAx = A'Sb -> Nx = r\n");
2247 At = gmt_M_memory (GMT, NULL, nm * nm, double);
2248 AtS = gmt_M_memory (GMT, NULL, nm * nm, double);
2249 S = gmt_M_memory (GMT, NULL, nm, double);
2250 /* 1. Transpose A and set diagonal matrix with squared weights (here a vector) S */
2251 GMT_Report (API, GMT_MSG_INFORMATION, "Create S = W'*W diagonal matrix, A', and compute A' * S\n");
2252 for (row = 0; row < nm; row++) {
2253 /* Set the diagonal using (=1/sigma^2) if given sigmas or the weights as given */
2254 S[row] = (Ctrl->W.mode == GSP_GOT_SIG) ? X[row][dimension] * X[row][dimension] : X[row][dimension];
2255 for (col = 0; col < nm; col++) {
2256 ij = row * nm + col;
2257 ji = col * nm + row;
2258 At[ji] = A[ij];
2259 }
2260 }
2261 /* 2. Compute AtS = At * S. This means scaling all terms in At columns by the corresponding S entry */
2262 for (row = ij = 0; row < nm; row++) {
2263 for (col = 0; col < nm; col++, ij++)
2264 AtS[ij] = At[ij] * S[col];
2265 }
2266 /* 3. Compute r = AtS * obs (but we recycle S to hold r) */
2267 GMT_Report (API, GMT_MSG_DEBUG, "Compute r = A'*S*b\n");
2268 gmt_matrix_matrix_mult (GMT, AtS, obs, nm, nm, 1U, S);
2269 /* 4. Compute N = AtS * A (but we recycle At to hold N) */
2270 GMT_Report (API, GMT_MSG_DEBUG, "Compute N = A'*S*A\n");
2271 gmt_matrix_matrix_mult (GMT, AtS, A, nm, nm, nm, At);
2272 /* Now free A, AtS and obs and let "A" be N and "obs" be r; these are the weighted normal equations */
2273 gmt_M_free (GMT, A); gmt_M_free (GMT, AtS); gmt_M_free (GMT, obs);
2274 A = At; obs = S;
2275 if (Ctrl->debug.active) greenspline_dump_system (A, obs, nm, "Normal equation N row || r");
2276 }
2277
2278 if (Ctrl->C.active) { /* Solve using SVD */
2279 int error;
2280
2281 GMT_Report (API, GMT_MSG_INFORMATION, "Solve linear equations by Singular Value Decomposition\n");
2282 #ifndef HAVE_LAPACK
2283 GMT_Report (API, GMT_MSG_WARNING, "Note: SVD solution without LAPACK will be very, very slow.\n");
2284 GMT_Report (API, GMT_MSG_WARNING, "We strongly recommend you install LAPACK and recompile GMT.\n");
2285 #endif
2286 v = gmt_M_memory (GMT, NULL, nm * nm, double);
2287 s = gmt_M_memory (GMT, NULL, nm, double);
2288 if ((error = gmt_svdcmp (GMT, A, (unsigned int)nm, (unsigned int)nm, s, v)) != 0) Return (error);
2289 if (Ctrl->C.history) { /* Keep copy of original singular values */
2290 ssave = gmt_M_memory (GMT, NULL, nm, double);
2291 gmt_M_memcpy (ssave, s, nm, double);
2292 }
2293 if (Ctrl->C.file) { /* Save the singular values for study */
2294 struct GMT_SINGULAR_VALUE *eigen = gmt_sort_svd_values (GMT, s, nm);
2295 uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, nm, 2};
2296 unsigned int col_type[3];
2297 struct GMT_DATASET *E = NULL;
2298 if ((E = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, e_dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2299 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving singular values\n");
2300 Return (API->error);
2301 }
2302 gmt_M_memcpy (col_type, GMT->current.io.col_type[GMT_OUT], 2, unsigned int); /* Save previous x/y output col types */
2303 GMT->current.io.col_type[GMT_OUT][GMT_X] = GMT->current.io.col_type[GMT_OUT][GMT_Y] = GMT_IS_FLOAT;
2304 for (i = 0; i < nm; i++) {
2305 E->table[0]->segment[0]->data[GMT_X][i] = i; /* Let 0 be x-value of the first eigenvalue */
2306 E->table[0]->segment[0]->data[GMT_Y][i] = eigen[i].value;
2307 }
2308 E->table[0]->segment[0]->n_rows = nm;
2309 if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, "index\teigenvalue", E)) {
2310 Return (API->error);
2311 }
2312 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, Ctrl->C.file, E) != GMT_NOERROR) {
2313 Return (API->error);
2314 }
2315 gmt_M_memcpy (GMT->current.io.col_type[GMT_OUT], col_type, 2, unsigned int); /* Restore output col types */
2316 GMT_Report (API, GMT_MSG_INFORMATION, "Eigenvalues saved to %s\n", Ctrl->C.file);
2317 gmt_M_free (GMT, eigen);
2318
2319 if (Ctrl->C.dryrun) { /* We are done */
2320 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
2321 gmt_M_free (GMT, X);
2322 gmt_M_free (GMT, s);
2323 gmt_M_free (GMT, v);
2324 gmt_M_free (GMT, A);
2325 gmt_M_free (GMT, obs);
2326 if (dimension == 2) gmt_free_grid (GMT, &Grid, true);
2327 Return (GMT_NOERROR);
2328 }
2329 }
2330 b = gmt_M_memory (GMT, NULL, nm, double);
2331 gmt_M_memcpy (b, obs, nm, double);
2332 n_use = gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, Ctrl->C.value, Ctrl->C.mode);
2333 if (n_use == -1) Return (GMT_RUNTIME_ERROR);
2334 GMT_Report (API, GMT_MSG_INFORMATION, "[%d of %" PRIu64 " eigen-values used]\n", n_use, nm);
2335
2336 if (Ctrl->C.history == GMT_SVD_NO_HISTORY) {
2337 gmt_M_free (GMT, s);
2338 gmt_M_free (GMT, v);
2339 gmt_M_free (GMT, b);
2340 }
2341 }
2342 else { /* Gauss-Jordan elimination */
2343 int error;
2344 if (gmt_M_is_zero (r_min)) {
2345 GMT_Report (API, GMT_MSG_ERROR, "Your matrix is singular because you have duplicate data constraints\n");
2346 GMT_Report (API, GMT_MSG_ERROR, "Preprocess your data with one of the blockm* modules to eliminate them\n");
2347
2348 }
2349 GMT_Report (API, GMT_MSG_INFORMATION, "Solve linear equations by Gauss-Jordan elimination\n");
2350 if (!Ctrl->W.active) GMT_Report (API, GMT_MSG_INFORMATION, "In the absence of weights, the solution will be exact\n");
2351 if ((error = gmt_gaussjordan (GMT, A, (unsigned int)nm, obs)) != 0) {
2352 GMT_Report (API, GMT_MSG_ERROR, "You probably have nearly duplicate data constraints\n");
2353 GMT_Report (API, GMT_MSG_ERROR, "Preprocess your data with one of the blockm* modules\n");
2354 Return (error);
2355 }
2356 }
2357 alpha = obs; /* Just a different name since the obs vector now holds the alpha factors */
2358
2359 if (Ctrl->M.active && Ctrl->M.mode == GMT_OUT) {
2360 /* EXPERIMENTAL and not completed - need normalization information, trend etc */
2361 bool was = GMT->current.setting.io_header[GMT_OUT]; /* Current setting */
2362 uint64_t m_dim[GMT_DIM_SIZE] = {1, 1, 0, 1}; /* Do not allocate any rows */
2363 char header[GMT_LEN64] = {""};
2364 struct GMT_DATASET *M = NULL;
2365
2366 if ((M = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, m_dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2367 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving misfit estimates\n");
2368 Return (API->error);
2369 }
2370 M->table[0]->segment[0]->n_rows = nm;
2371 M->table[0]->segment[0]->data[GMT_X] = alpha;
2372
2373 sprintf (header, "N: %" PRIu64 " S: %s G: %s", nm, (Ctrl->C.active) ? "SVD" : "G-J", Ctrl->S.arg);
2374 gmt_insert_tableheader (GMT, M->table[0], header);
2375 GMT->current.setting.io_header[GMT_OUT] = true;
2376
2377 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, Ctrl->M.file, M) != GMT_NOERROR) {
2378 Return (API->error);
2379 }
2380 M->table[0]->segment[0]->data[GMT_X] = NULL; /* Since we did not allocate that array */
2381 GMT->current.setting.io_header[GMT_OUT] = was; /* Restore default */
2382 }
2383
2384 if (Ctrl->C.history == GMT_SVD_NO_HISTORY) gmt_M_free (GMT, A);
2385
2386 if (Ctrl->E.active && Ctrl->C.history == GMT_SVD_NO_HISTORY) {
2387 double value, mean = 0.0, std = 0.0, rms = 0.0, dev, chi2, chi2_sum = 0, pvar_sum = 0.0, *predicted = NULL;
2388 uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, nm, dimension+3+Ctrl->W.active};
2389 unsigned int m = 0;
2390 struct GMT_DATASET *E = NULL;
2391 struct GMT_DATASEGMENT *S = NULL;
2392 if (Ctrl->E.mode) { /* Want to write out prediction errors */
2393 if ((E = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, e_dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2394 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving misfit estimates\n");
2395 Return (API->error);
2396 }
2397 S = E->table[0]->segment[0];
2398 S->n_rows = nm;
2399 }
2400 predicted = gmt_M_memory (GMT, NULL, nm, double); /* To hold predictions */
2401 gmt_matrix_matrix_mult (GMT, A_orig, alpha, nm, nm, 1U, predicted); /* predicted = A * alpha are normalized predictions at data points */
2402 for (j = 0; j < nm; j++) { /* For each data constraint */
2403 predicted[j] = greenspline_undo_normalization (X[j], predicted[j], normalize, norm, dimension); /* undo normalization first */
2404 pvar_sum += predicted[j] * predicted[j]; /* Sum of predicted variance */
2405 dev = value = orig_obs[j] - predicted[j]; /* Deviation between observed and predicted */
2406 rms += dev * dev; /* Accumulate rms sum */
2407 if (Ctrl->W.active) { /* If data had uncertainties we also compute the chi2 sum */
2408 chi2 = pow (dev * X[j][dimension], 2.0);
2409 chi2_sum += chi2;
2410 }
2411 /* Use Welford (1962) algorithm to compute mean and variance */
2412 m++;
2413 dev = value - mean;
2414 mean += dev / m;
2415 std += dev * (value - mean);
2416 if (Ctrl->E.mode) { /* Store assessment for each observation in misfit table */
2417 for (p = 0; p < dimension; p++) /* Duplicate point coordinates */
2418 S->data[p][j] = X[j][p];
2419 S->data[p++][j] = orig_obs[j];
2420 S->data[p++][j] = predicted[j];
2421 S->data[p][j] = value;
2422 if (Ctrl->W.active)
2423 S->data[++p][j] = chi2;
2424 }
2425 }
2426 rms = sqrt (rms / nm);
2427 std = (m > 1) ? sqrt (std / (m-1.0)) : GMT->session.d_NaN;
2428 if (Ctrl->W.active)
2429 GMT_Report (API, GMT_MSG_INFORMATION, "Misfit evaluation: N = %u\tMean = %g\tStd.dev = %g\tRMS = %g\tChi^2 = %g\n", nm, mean, std, rms, chi2_sum);
2430 else
2431 GMT_Report (API, GMT_MSG_INFORMATION, "Misfit evaluation: N = %u\tMean = %g\tStd.dev = %g\tRMS = %g\n", nm, mean, std, rms);
2432 GMT_Report (API, GMT_MSG_INFORMATION, "Variance evaluation: Data = %g\tModel = %g\tExplained = %5.1lf %%\n", var_sum, pvar_sum, 100.0 * pvar_sum / var_sum);
2433 gmt_M_free (GMT, orig_obs); gmt_M_free (GMT, predicted); gmt_M_free (GMT, A_orig);
2434 if (Ctrl->E.mode) { /* Want to write out prediction errors */
2435 char header[GMT_LEN64] = {""};
2436 if (gmt_M_x_is_lon (GMT, GMT_IN))
2437 sprintf (header, "#lon\t");
2438 else
2439 sprintf (header, "#x\t");
2440 if (dimension > 1) strcat (header, gmt_M_y_is_lat (GMT, GMT_IN) ? "lat\t" : "y\t");
2441 if (dimension > 2) strcat (header, "z\t");
2442 strcat (header, "obs\tpredict\tdev");
2443 if (Ctrl->W.active) strcat (header, "\tchi2");
2444 if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, E)) {
2445 Return (API->error);
2446 }
2447 gmt_set_tableheader (API->GMT, GMT_OUT, true); /* So header is written */
2448 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, Ctrl->E.file, E) != GMT_NOERROR) {
2449 Return (API->error);
2450 }
2451 }
2452 }
2453
2454 if (Ctrl->N.file || dimension == 1 || write_3D_records) Rec = gmt_new_record (GMT, NULL, NULL);
2455
2456 if (Ctrl->N.file) { /* Specified nodes only */
2457 unsigned int wmode = GMT_ADD_DEFAULT;
2458 double out[4];
2459
2460 /* Must register Ctrl->G.file first since we are going to writing rec-by-rec */
2461 if (Ctrl->G.active) {
2462 if ((out_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_OUT, NULL, Ctrl->G.file)) == GMT_NOTSET)
2463 Return (API->error);
2464 wmode = GMT_ADD_EXISTING;
2465 }
2466 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, wmode, 0, options) != GMT_NOERROR) { /* Establishes output */
2467 Return (API->error);
2468 }
2469 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
2470 Return (API->error);
2471 }
2472 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
2473 Return (API->error);
2474 }
2475 if ((error = GMT_Set_Columns (API, GMT_OUT, dimension + 1, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
2476 Return (error);
2477 }
2478 gmt_M_memset (out, 4, double);
2479 Rec->data = out;
2480 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline at %" PRIu64 " given locations\n", T->n_records);
2481 for (seg = 0; seg < T->n_segments; seg++) {
2482 for (row = 0; row < T->segment[seg]->n_rows; row++) {
2483 for (ii = 0; ii < dimension; ii++) out[ii] = T->segment[seg]->data[ii][row];
2484 if (T->segment[seg]->text) Rec->text = T->segment[seg]->text[row];
2485 out[dimension] = 0.0;
2486 for (p = 0; p < nm; p++) {
2487 r = greenspline_get_radius (GMT, out, X[p], dimension);
2488 if (Ctrl->Q.active) {
2489 C = greenspline_get_dircosine (GMT, Ctrl->Q.dir, out, X[p], dimension, false);
2490 part = dGdr (GMT, r, par, Lz) * C;
2491 }
2492 else
2493 part = G (GMT, r, par, Lz);
2494 out[dimension] += alpha[p] * part;
2495 }
2496 out[dimension] = greenspline_undo_normalization (out, out[dimension], normalize, norm, dimension);
2497 GMT_Put_Record (API, GMT_WRITE_DATA, Rec);
2498 }
2499 }
2500 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
2501 Return (API->error);
2502 }
2503 if (GMT_Destroy_Data (API, &Nin) != GMT_NOERROR) {
2504 Return (API->error);
2505 }
2506 gmt_fclose (GMT, fp);
2507 }
2508 else { /* Output on equidistant lattice */
2509 uint64_t nz_off, nxy;
2510 unsigned int layer, wmode = GMT_ADD_DEFAULT;
2511 double *xp = NULL, *yp = NULL, wp, V[4] = {0.0, 0.0, 0.0, 0.0};
2512 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline at %" PRIu64 " equidistant output locations\n", n_ok);
2513 /* Precalculate coordinates */
2514 xp = gmt_grd_coord (GMT, header, GMT_X);
2515 if (dimension > 1) yp = gmt_grd_coord (GMT, header, GMT_Y);
2516 nxy = header->size; /* Will only be used for 3-D anyway when there are layers */
2517 if (dimension == 1 || write_3D_records) { /* Write ASCII table to named file or stdout for 1-D or stdout for 3-D */
2518 GMT->common.b.ncol[GMT_OUT] = dimension + 1;
2519 if (Ctrl->G.active) {
2520 if ((out_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_OUT, NULL, Ctrl->G.file)) == GMT_NOTSET) {
2521 gmt_M_free (GMT, xp);
2522 Return (error);
2523 }
2524 wmode = GMT_ADD_EXISTING;
2525 }
2526 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, wmode, 0, options) != GMT_NOERROR) { /* Establishes output */
2527 gmt_M_free (GMT, xp);
2528 Return (API->error);
2529 }
2530 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
2531 gmt_M_free (GMT, xp);
2532 Return (API->error);
2533 }
2534 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
2535 gmt_M_free (GMT, xp);
2536 Return (API->error);
2537 }
2538
2539 } /* Else we are writing a grid or cube */
2540 if (Ctrl->C.history) { /* 2-D only: Write out grid after adding contribution for each eigenvalue separately */
2541 /* Note: Because the SVD decomposition is not sorting the vectors from largest to smallest eigenvalue the
2542 * gmt_solve_svd sets to zero those we don't want but we must still loop over its full length to ensure we
2543 * include the eigenvalues we want. */
2544 unsigned int width = urint (floor (log10 ((double)n_use))) + 1; /* Width of maximum integer needed */
2545 int64_t e, col, row, p; /* On Windows, the 'for' index variables must be signed, so redefine these 3 inside this block only */
2546 gmt_grdfloat *current = NULL, *previous = NULL;
2547 double l2_sum_n = 0.0, l2_sum_e = 0.0, predicted;
2548 static char *mkind[3] = {"", "Incremental", "Cumulative"};
2549 char file[PATH_MAX] = {""};
2550 struct GMT_SINGULAR_VALUE *eigen = NULL;
2551 struct GMT_DATASET *E = NULL;
2552 struct GMT_DATASEGMENT *S = NULL;
2553
2554 current = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat);
2555 if (Ctrl->C.history & GMT_SVD_INCREMENTAL)
2556 previous = gmt_M_memory_aligned (GMT, NULL, Out->header->size, gmt_grdfloat);
2557 gmt_grd_init (GMT, Out->header, options, true);
2558 if (Ctrl->E.active) { /* Want to write out misfit as function of eigenvalue */
2559 uint64_t e_dim[GMT_DIM_SIZE] = {1, 1, n_use, 4+Ctrl->W.active};
2560 eigen = gmt_sort_svd_values (GMT, ssave, nm); /* Get sorted eigenvalues */
2561 if ((E = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_NONE, 0, e_dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2562 GMT_Report (API, GMT_MSG_ERROR, "Unable to create a data set for saving misfit estimates per eigenvector\n");
2563 Return (API->error);
2564 }
2565 for (k = 0; k < nm; k++) /* Get sum of squared eigenvalues */
2566 l2_sum_n += eigen[k].value * eigen[k].value;
2567 S = E->table[0]->segment[0];
2568 S->n_rows = n_use;
2569 }
2570
2571 for (e = 0; e < n_use; e++) { /* Only loop over the first n_use eigenvalues (if restricted) */
2572 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)e);
2573 gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */
2574 (void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)e, GMT_SVD_EIGEN_NUMBER_CUTOFF);
2575 /* obs (hence alpha) now has the solution for the coefficients based on the first e eigenvalues */
2576 if (Ctrl->Q.active) { /* Derivatives of solution */
2577 #ifdef _OPENMP
2578 #pragma omp parallel for private(row,V,col,ij,p,wp,r,C,part) shared(Grid,yp,xp,nm,GMT,Ctrl,X,G,par,Lz,alpha,Out,normalize,norm)
2579 #endif
2580 for (row = 0; row < Grid->header->n_rows; row++) {
2581 V[GMT_Y] = yp[row];
2582 for (col = 0; col < Grid->header->n_columns; col++) {
2583 ij = gmt_M_ijp (Grid->header, row, col);
2584 if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */
2585 V[GMT_X] = xp[col];
2586 /* Here, V holds the current output coordinates */
2587 for (p = 0, wp = 0.0; p < nm; p++) {
2588 r = greenspline_get_radius (GMT, V, X[p], 2U);
2589 C = greenspline_get_dircosine (GMT, Ctrl->Q.dir, V, X[p], 2U, false);
2590 part = dGdr (GMT, r, par, Lz) * C;
2591 wp += alpha[p] * part;
2592 }
2593 V[GMT_Z] = greenspline_undo_normalization (V, wp, normalize, norm, 2U);
2594 Out->data[ij] = (gmt_grdfloat)V[GMT_Z];
2595 }
2596 }
2597 }
2598 else { /* Surface solution */
2599 if (Ctrl->E.active) { /* Compute the history of model misfit */
2600 double dev, rms = 0.0, chi2_sum = 0.0;
2601 for (j = 0; j < nm; j++) { /* For each data constraint */
2602 for (p = 0, wp = 0.0; p < nm; p++) { /* Add contribution for each data constraint */
2603 r = greenspline_get_radius (GMT, X[j], X[p], 2U);
2604 part = G (GMT, r, par, Lz);
2605 wp += alpha[p] * part; /* Just add this scaled Green's function */
2606 }
2607 predicted = greenspline_undo_normalization (X[j], wp, normalize, norm, dimension); /* Undo normalization first */
2608 dev = orig_obs[j] - predicted; /* Deviation between observed and predicted */
2609 dev *= dev; /* Squared misfit */
2610 rms += dev; /* Accumulate rms sum */
2611 if (Ctrl->W.active) { /* If data had uncertainties we also compute the chi2 sum */
2612 double chi2 = dev * pow (X[j][2], 2.0);
2613 chi2_sum += chi2;
2614 }
2615 }
2616 rms = sqrt (rms / nm);
2617 l2_sum_e += eigen[e].value * eigen[e].value;
2618 if (Ctrl->W.active)
2619 GMT_Report (API, GMT_MSG_INFORMATION, "Cumulative data misfit for eigenvalue # %d: rms = %lg chi2 = %lg\n", (int)e, rms, chi2_sum);
2620 else
2621 GMT_Report (API, GMT_MSG_INFORMATION, "Cumulative data misfit for eigenvalue # %d: rms = %lg\n", (int)e, rms);
2622 S->data[0][e] = e; /* Eigenvalue number (starting at 0) */
2623 S->data[1][e] = eigen[e].value; /* Eigenvalue, from largest to smallest */
2624 S->data[2][e] = 100.0 * l2_sum_e / l2_sum_n; /* Percent of model variance */
2625 S->data[3][e] = rms; /* RMS misfit for this solution */
2626 if (Ctrl->W.active) S->data[4][e] = chi2_sum; /* Chi^2 sum for this solution */
2627 }
2628 #ifdef _OPENMP
2629 #pragma omp parallel for private(row,V,col,ij,p,wp,r,part) shared(Grid,yp,xp,nm,GMT,X,G,par,Lz,alpha,Out,normalize,norm)
2630 #endif
2631 for (row = 0; row < Grid->header->n_rows; row++) {
2632 V[GMT_Y] = yp[row];
2633 for (col = 0; col < Grid->header->n_columns; col++) {
2634 ij = gmt_M_ijp (Grid->header, row, col);
2635 if (gmt_M_is_fnan (Grid->data[ij])) continue; /* Only do solution where mask is not NaN */
2636 V[GMT_X] = xp[col];
2637 /* Here, V holds the current output coordinates */
2638 for (p = 0, wp = 0.0; p < nm; p++) {
2639 r = greenspline_get_radius (GMT, V, X[p], 2U);
2640 part = G (GMT, r, par, Lz);
2641 wp += alpha[p] * part;
2642 }
2643 Out->data[ij] = (gmt_grdfloat)greenspline_undo_normalization (V, wp, normalize, norm, 2U);
2644 }
2645 } /* End of row-loop [OpenMP] */
2646 }
2647 gmt_M_memcpy (current, Out->data, Out->header->size, gmt_grdfloat); /* Save current solution */
2648
2649 if (Ctrl->C.history & GMT_SVD_CUMULATIVE) { /* Write out the cumulative solution first */
2650 if (strchr (Ctrl->G.file, '%')) /* Gave a template, use it to write one of the two types of grids */
2651 sprintf (file, Ctrl->G.file, (int)e+1);
2652 else /* Create the appropriate cumulative gridfile name from static file name */
2653 greenspline_set_filename (Ctrl->G.file, e, width, GMT_SVD_CUMULATIVE, file);
2654 snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s). %s contribution for eigenvalue # %d", method[Ctrl->S.mode], Ctrl->S.arg, mkind[GMT_SVD_CUMULATIVE], (int)e+1);
2655 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out))
2656 Return (API->error); /* Update solution for e eigenvalues only */
2657 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Out) != GMT_NOERROR)
2658 Return (API->error);
2659 }
2660 if (Ctrl->C.history & GMT_SVD_INCREMENTAL) { /* Want to write out incremental solution due to this eigenvalue */
2661 gmt_M_grd_loop (GMT, Out, row, col, ij) Out->data[ij] = current[ij] - previous[ij]; /* Incremental improvement since last time */
2662 gmt_M_memcpy (previous, current, Out->header->size, gmt_grdfloat); /* Save current solution which will be previous for next eigenvalue */
2663 if (strchr (Ctrl->G.file, '%')) /* Gave a template, use it to write one of the two types of grids */
2664 sprintf (file, Ctrl->G.file, (int)e+1);
2665 else /* Create the appropriate cumulative gridfile name from static file name */
2666 greenspline_set_filename (Ctrl->G.file, e, width, GMT_SVD_INCREMENTAL, file);
2667 snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s). %s contribution for eigenvalue # %d", method[Ctrl->S.mode], Ctrl->S.arg, mkind[GMT_SVD_INCREMENTAL], (int)e+1);
2668 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out))
2669 Return (API->error); /* Update solution for e eigenvalues only */
2670 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, file, Out) != GMT_NOERROR)
2671 Return (API->error);
2672 }
2673 }
2674 if (Ctrl->E.active) { /* Compute the history of model misfit as rms */
2675 char header[GMT_LEN64] = {""};
2676 sprintf (header, "# eigenno\teigenval\tvar_percent\trms");
2677 if (Ctrl->W.active) strcat (header, "\tchi2");
2678 if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, E)) {
2679 Return (API->error);
2680 }
2681 gmt_set_tableheader (API->GMT, GMT_OUT, true); /* So header is written */
2682 for (k = 0; k < 5; k++) GMT->current.io.col_type[GMT_OUT][k] = GMT_IS_FLOAT; /* Set plain float column types */
2683 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_WRITE_SET, NULL, Ctrl->E.file, E) != GMT_NOERROR) {
2684 Return (API->error);
2685 }
2686 gmt_M_free (GMT, eigen);
2687 gmt_M_free (GMT, orig_obs);
2688 }
2689
2690 /* Free temporary arrays */
2691 gmt_M_free_aligned (GMT, current);
2692 if (Ctrl->C.history & GMT_SVD_INCREMENTAL)
2693 gmt_M_free_aligned (GMT, previous);
2694 gmt_M_free (GMT, A);
2695 gmt_M_free (GMT, s);
2696 gmt_M_free (GMT, v);
2697 gmt_M_free (GMT, b);
2698 gmt_M_free (GMT, ssave);
2699 }
2700 else { /* Just compute the final interpolation */
2701 if (dimension == 1 || write_3D_records) Rec->data = V; /* For rec-by-rec output */
2702 for (layer = 0, nz_off = 0; layer < n_layers; layer++, nz_off += nxy) { /* Might be dummy loop of 1 layer unless 3-D */
2703 int64_t col, row, p; /* On Windows, the 'for' index variables must be signed, so redefine these 3 inside this block only */
2704 double z_layer = (dimension == 3) ? Cube->z[layer] : 0.0;
2705 if (Ctrl->Q.active) { /* Derivatives of solution */
2706 #ifdef _OPENMP
2707 #pragma omp parallel for private(row,V,ij,col,p,wp,r,C,part) shared(header,dimension,yp,z_layer,nz_off,data,xp,nm,GMT,X,Ctrl,dGdr,par,Lz,alpha,normalize,norm)
2708 #endif
2709 for (row = 0; row < header->n_rows; row++) { /* This would be a dummy loop for 1 row if 1-D data */
2710 if (dimension > 1) V[GMT_Y] = yp[row];
2711 if (dimension == 3) V[GMT_Z] = z_layer;
2712 ij = (dimension > 1) ? gmt_M_ijp (header, row, 0) + nz_off : 0;
2713 for (col = 0; col < header->n_columns; col++, ij++) { /* This loop is always active for 1,2,3D */
2714 if (dimension == 2 && gmt_M_is_fnan (data[ij])) continue; /* Only do solution where mask is not NaN */
2715 V[GMT_X] = xp[col];
2716 /* Here, V holds the current output coordinates */
2717 for (p = 0, wp = 0.0; p < (int64_t)nm; p++) { /* Loop over Green's function components */
2718 r = greenspline_get_radius (GMT, V, X[p], dimension);
2719 C = greenspline_get_dircosine (GMT, Ctrl->Q.dir, V, X[p], dimension, false);
2720 part = dGdr (GMT, r, par, Lz) * C;
2721 wp += alpha[p] * part;
2722 }
2723 data[ij] = (gmt_grdfloat)greenspline_undo_normalization (V, wp, normalize, norm, dimension);
2724 }
2725 } /* End of row-loop [OpenMP] */
2726 }
2727 else { /* Regular surface */
2728 #ifdef _OPENMP
2729 #pragma omp parallel for private(row,V,ij,col,p,wp,r,part) shared(header,dimension,yp,z_layer,nz_off,data,xp,nm,GMT,X,G,par,Lz,alpha,normalize,norm)
2730 #endif
2731 for (row = 0; row < header->n_rows; row++) { /* This would be a dummy loop for 1 row if 1-D data */
2732 if (dimension > 1) V[GMT_Y] = yp[row];
2733 if (dimension == 3) V[GMT_Z] = z_layer;
2734 ij = (dimension > 1) ? gmt_M_ijp (header, row, 0) + nz_off : 0;
2735 for (col = 0; col < header->n_columns; col++, ij++) { /* This loop is always active for 1,2,3D */
2736 if (dimension == 2 && gmt_M_is_fnan (data[ij])) continue; /* Only do solution where mask is not NaN */
2737 V[GMT_X] = xp[col];
2738 /* Here, V holds the current output coordinates */
2739 for (p = 0, wp = 0.0; p < (int64_t)nm; p++) { /* Loop over Green's function components */
2740 r = greenspline_get_radius (GMT, V, X[p], dimension);
2741 part = G (GMT, r, par, Lz);
2742 wp += alpha[p] * part;
2743 }
2744 data[ij] = (gmt_grdfloat)greenspline_undo_normalization (V, wp, normalize, norm, dimension);
2745 }
2746 } /* End of row-loop [OpenMP] */
2747 }
2748
2749 if (write_3D_records) { /* Must dump this slice of the 3-D cube as ASCII slices as a backwards compatibility option */
2750 V[GMT_Z] = z_layer;
2751 for (row = 0; row < header->n_rows; row++) {
2752 V[GMT_Y] = yp[row];
2753 for (col = 0; col < header->n_columns; col++) {
2754 V[GMT_X] = xp[col];
2755 ij = gmt_M_ijp (header, row, col) + nz_off;
2756 V[dimension] = data[ij];
2757 GMT_Put_Record (API, GMT_WRITE_DATA, Rec);
2758 }
2759 }
2760 }
2761 } /* End of layer loop */
2762
2763 /* Time to write output */
2764 if (dimension == 1) { /* Must dump 1-D records */
2765 for (col = 0; col < header->n_columns; col++) {
2766 V[GMT_X] = xp[col];
2767 V[GMT_Y] = data[col];
2768 GMT_Put_Record (API, GMT_WRITE_DATA, Rec);
2769 }
2770 gmt_M_free (GMT, data);
2771 }
2772 else if (dimension == 2) { /* Write the 2-D grid */
2773 gmt_grd_init (GMT, Out->header, options, true);
2774 if (Ctrl->D.active && gmt_decode_grd_h_info (GMT, Ctrl->D.information, Out->header)) {
2775 Return (GMT_PARSE_ERROR);
2776 }
2777 snprintf (Out->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s)", method[Ctrl->S.mode], Ctrl->S.arg);
2778 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Out)) Return (API->error);
2779 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Out) != GMT_NOERROR) {
2780 Return (API->error);
2781 }
2782 }
2783 else if (dimension == 3 && !write_3D_records) { /* Write the 3-D cube */
2784 gmt_grd_init (GMT, Cube->header, options, true);
2785 snprintf (Cube->header->remark, GMT_GRID_REMARK_LEN160, "%s (-S%s)", method[Ctrl->S.mode], Ctrl->S.arg);
2786 if (GMT_Set_Comment (API, GMT_IS_CUBE, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Cube)) Return (API->error);
2787 if (GMT_Write_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Cube))
2788 Return (EXIT_FAILURE);
2789 }
2790 }
2791 if (delete_grid) /* No longer required for 1-D and 3-D */
2792 gmt_free_grid (GMT, &Grid, false);
2793 if (dimension == 3) GMT_Destroy_Data (API, &Cube); /* Done with the output cube */
2794
2795 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
2796 Return (API->error);
2797 }
2798 gmt_M_free (GMT, xp);
2799 if (dimension > 1) gmt_M_free (GMT, yp);
2800 }
2801
2802 /* Clean up */
2803
2804 for (p = 0; p < nm; p++) gmt_M_free (GMT, X[p]);
2805 gmt_M_free (GMT, X);
2806 gmt_M_free (GMT, obs);
2807 if (m) {
2808 for (p = 0; p < m; p++) gmt_M_free (GMT, D[p]);
2809 gmt_M_free (GMT, D);
2810 }
2811 if (Rec) gmt_M_free (GMT, Rec);
2812 greenspline_free_lookup (GMT, &Lz, 0);
2813 greenspline_free_lookup (GMT, &Lg, 1);
2814
2815 Return (GMT_NOERROR);
2816 }
2817