1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1999-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU Lesser General Public License as published by
7 * the Free Software Foundation; version 3 or any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU Lesser General Public License for more details.
13 *
14 * Contact info: www.generic-mapping-tools.org
15 *--------------------------------------------------------------------*/
16 /*
17 * HOTSPOTTER will (1) read ASCII file(s) with records for each seamount
18 * (2) read an ASCII file with stage or total reconstruction rotations,
19 * (3) convolve the flowline of each seamount with its gravimetric shape, and
20 * (4) build a cumulative volcano amplitude (CVA) grid. The grid is
21 * written out in GMT format and can be processed and plotted with GMT.
22 * HOTSPOTTER is part of the SPOTTER supplemental GMT package and should
23 * be installed accordingly, see README.spotter.
24 *
25 * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
26 * Date: 22-JUN-1999
27 * Version: 1.0
28 *
29 *-------------------------------------------------------------------------
30 * An ASCII stage pole (Euler) file must have following format:
31 *
32 * 1. Any number of comment lines starting with # in first column
33 * 2. Any number of blank lines (just carriage return, no spaces)
34 * 2. Any number of stage pole records which each have the format:
35 * lon(deg) lat(deg) tstart(Ma) tstop(Ma) ccw-angle(deg)
36 * 3. stage records must go from oldest to youngest rotation
37 * 4. Note tstart is larger (older) that tstop for each record
38 * 5. No gaps allowed: tstart must equal the previous records tstop
39 *
40 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
41 *
42 * # Time in Ma, angles in degrees
43 * # lon lat tstart tend ccw-angle
44 * 165 85 150 100 24.0
45 * 284 36 100 74 15.0
46 * 265 22 74 65 7.5
47 * 253 17 65 42 14.0
48 * 285 68 42 0 34.0
49 *
50 * AN ASCII total reconstruction file must have the following format:
51 *
52 * 1. Any number of comment lines starting with # in first column
53 * 2. Any number of blank lines (just carriage return, no spaces)
54 * 2. Any number of finite pole records which each have the format:
55 * lon(deg) lat(deg) tstop(Ma) ccw-angle(deg)
56 * 3. total reconstructions rotations must go from youngest to oldest
57 *
58 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
59 *
60 * # Time in Ma, angles in degrees
61 * #longitude latitude time(My) angle(deg)
62 * 285.00000 68.00000 42.0000 34.0000
63 * 275.66205 53.05082 65.0000 43.5361
64 * 276.02501 48.34232 74.0000 50.0405
65 * 279.86436 46.30610 100.0000 64.7066
66 * 265.37800 55.69932 150.0000 82.9957
67 *
68 *
69 * Seamount data file(s) must have the following format:
70 * 1. Any number of comment lines starting with # in first column
71 * 2. Any number of blank lines (just carriage return, no spaces)
72 * 3. Any number of seamount records which each have the format:
73 * lon(deg) lat(deg) amplitude radius(km) age(Ma)
74 * 4. The amplitude is in user units (m, mGal, km^3, whatever)
75 * 5. Age represents the upper possible age for seamount, which
76 * is usually the age of the seafloor beneath it. If the
77 * crustal age is not known, set it to NaN in the file.
78 * NaN values will be replaced with -N (see usage)
79 *
80 * Example: Wessel & Lyons [1997] Pacific seamounts (just a few records):
81 * # Pacific seamounts > 100 Eotvos amplitude in Vertical Gravity Gradient
82 * # From Wessel & Lyons [1997]
83 * #LON FAA VGG RADIUS CRUST_AGE
84 * 134.38333 0.9436415 120.5 22.97 37.606796
85 * 136.05 7.6325042 102.7 18.67 NaN
86 * 131.28333 1.1423035 129.0 17.16 NaN
87 * .....
88 *
89 *------------------------------------------------------------------------------
90 * REFERENCES:
91 *
92 * -> The hotspotting technique:
93 *
94 * Aslanian, D., L. Geli, and J.-L. Olivet, 1998, Hotspotting called into
95 * question, Nature, 396, 127.
96 * Wessel, P., and L. W. Kroenke, 1997, A geometric technique for
97 * relocating hotspots and refining absolute plate motions, Nature,
98 * 387, 365-369.
99 * Wessel, P., and L. W. Kroenke, 1998a, Factors influencing the locations
100 * of hot spots determined by the hot-spotting technique, Geophys. Res.
101 * Lett., 25, 55-58.
102 * Wessel, P., and L. W. Kroenke, 1998b, The geometric relationship between
103 * hot spots and seamounts: implications for Pacific hot spots, Earth
104 * Planet. Sci. Lett., 158, 1-18.
105 * Wessel, P., and L. W. Kroenke, 1998c, Hotspotting called into question
106 * - Reply, Nature, 396, 127-128.
107 *
108 * -> Seamount data set:
109 *
110 * Wessel, P., and S. Lyons, 1997, Distribution of large Pacific seamounts
111 * from Geosat/ERS-1: Implications for the history of intraplate volcanism,
112 * J. Geophys. Res., 102,22,459-22,475.
113 * Wessel, P., 1997, Sizes and ages of seamounts using remote sensing:
114 * Implications for intraplate volcanism, Science, 277, 802-805.
115 *
116 * -> Plate motion models (stage poles):
117 *
118 * Duncan, R.A., and D. Clague, 1985, Pacific plate motion recorded by linear
119 * volcanic chains, in: A.E.M. Nairn, F. G. Stehli, S. Uyeda (eds.), The
120 * Ocean Basins and Margins, Vol. 7A, Plenum, New York, pp. 89-121.
121 * Wessel and Kroenke, 1997, (see above)
122 *
123 */
124
125 #include "gmt_dev.h"
126 #include "spotter.h"
127
128 #define THIS_MODULE_CLASSIC_NAME "hotspotter"
129 #define THIS_MODULE_MODERN_NAME "hotspotter"
130 #define THIS_MODULE_LIB "spotter"
131 #define THIS_MODULE_PURPOSE "Create CVA grid from seamount locations"
132 #define THIS_MODULE_KEYS "<D{,GG}"
133 #define THIS_MODULE_NEEDS "R"
134 #define THIS_MODULE_OPTIONS "-:>RVbdeghiqr" GMT_OPT("FHMm")
135
136 struct HOTSPOTTER_CTRL { /* All control options for this program (except common args) */
137 /* active is true if the option has been activated */
138 struct HOTSPOTTER_D { /* -D<factor> */
139 bool active;
140 double value; /* Resampling factor */
141 } D;
142 struct HOTSPOTTER_E { /* -E[+]rotfile */
143 bool active;
144 bool mode;
145 char *file;
146 } E;
147 struct HOTSPOTTER_G { /* -Goutfile */
148 bool active;
149 char *file;
150 } G;
151 struct HOTSPOTTER_I { /* -I (for checking only) */
152 bool active;
153 } I;
154 struct HOTSPOTTER_N { /* -N */
155 bool active;
156 double t_upper;
157 } N;
158 struct HOTSPOTTER_S { /* -S */
159 bool active;
160 char *file;
161 } S;
162 struct HOTSPOTTER_T { /* -T<tzero> */
163 bool active;
164 double t_zero; /* Set zero age*/
165 } T;
166 };
167
New_Ctrl(struct GMT_CTRL * GMT)168 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
169 struct HOTSPOTTER_CTRL *C;
170
171 C = gmt_M_memory (GMT, NULL, 1, struct HOTSPOTTER_CTRL);
172
173 /* Initialize values whose defaults are not 0/false/NULL */
174
175 C->D.value = 0.5; /* Sets how many points along flowline per node */
176 C->N.t_upper = 180.0; /* Upper age assigned to seamounts on undated seafloor */
177 return (C);
178 }
179
Free_Ctrl(struct GMT_CTRL * GMT,struct HOTSPOTTER_CTRL * C)180 static void Free_Ctrl (struct GMT_CTRL *GMT, struct HOTSPOTTER_CTRL *C) { /* Deallocate control structure */
181 if (!C) return;
182 gmt_M_str_free (C->E.file);
183 gmt_M_str_free (C->G.file);
184 gmt_M_str_free (C->S.file);
185 gmt_M_free (GMT, C);
186 }
187
usage(struct GMTAPI_CTRL * API,int level)188 static int usage (struct GMTAPI_CTRL *API, int level) {
189 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
190 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
191 GMT_Usage (API, 0, "usage: %s [<table>] %s -G%s %s %s [-D<factor>] [-N<upper_age>] "
192 "[-S] [-T] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
193 name, SPOTTER_E_OPT, GMT_OUTGRID, GMT_Id_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT,
194 GMT_qi_OPT, GMT_r_OPT, GMT_colon_OPT, GMT_PAR_OPT);
195
196 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
197
198 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
199 GMT_Usage (API, 1, "\n<table> (in ASCII, binary, or netCDF) has 5 or more columns. If no file(s) is given, "
200 "standard input is read. Expects (x,y,z,r,t) records, with t in Ma.");
201 spotter_rot_usage (API);
202 gmt_outgrid_syntax (API, 'G', "Specify file name for output CVA grid");
203 GMT_Usage (API, 1, "\n%s", GMT_Id_OPT);
204 GMT_Usage (API, -2, "Specify grid interval(s); Append m [or s] to <dx> and/or <dy> for minutes [or seconds].");
205 GMT_Option (API, "Rg");
206 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
207 GMT_Usage (API, 1, "\n-D<factor>");
208 GMT_Usage (API, -2, "Scale affecting distance between points along flowline [0.5].");
209 GMT_Usage (API, 1, "\n-N<upper_age>");
210 GMT_Usage (API, -2, "Set upper age in m.y. for seamounts whose plate age is NaN [180].");
211 GMT_Usage (API, 1, "\n-S Normalize CVA grid to percentages of the CVA maximum.");
212 GMT_Usage (API, 1, "\n-T Truncate all ages to max age in stage pole model [Default extrapolates].");
213 GMT_Option (API, "V,bi5,di,e,h,i,qi,r,:,.");
214
215 return (GMT_MODULE_USAGE);
216 }
217
parse(struct GMT_CTRL * GMT,struct HOTSPOTTER_CTRL * Ctrl,struct GMT_OPTION * options)218 static int parse (struct GMT_CTRL *GMT, struct HOTSPOTTER_CTRL *Ctrl, struct GMT_OPTION *options) {
219 /* This parses the options provided to hotspotter and sets parameters in CTRL.
220 * Any GMT common options will override values set previously by other commands.
221 * It also replaces any file names specified as input or output with the data ID
222 * returned when registering these sources/destinations with the API.
223 */
224
225 unsigned int n_errors = 0, k;
226 struct GMT_OPTION *opt = NULL;
227 struct GMTAPI_CTRL *API = GMT->parent;
228
229 for (opt = options; opt; opt = opt->next) {
230 switch (opt->option) {
231
232 case '<': /* Skip input files */
233 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;
234 break;
235
236 /* Supplemental parameters */
237 case 'C': /* Now done automatically in spotter_init */
238 if (gmt_M_compat_check (GMT, 4))
239 GMT_Report (API, GMT_MSG_COMPAT, "-C is no longer needed as total reconstruction vs stage rotation is detected automatically.\n");
240 else
241 n_errors += gmt_default_error (GMT, opt->option);
242 break;
243
244 case 'D':
245 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
246 Ctrl->D.active = true;
247 Ctrl->D.value = atof (opt->arg);
248 break;
249 case 'E':
250 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
251 Ctrl->E.active = true; k = 0;
252 if (opt->arg[0] == '+') { Ctrl->E.mode = true; k = 1;}
253 if (opt->arg[k]) Ctrl->E.file = strdup (&opt->arg[k]);
254 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->E.file))) n_errors++;
255 break;
256 case 'G':
257 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
258 Ctrl->G.active = true;
259 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
260 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_REMOTE, &(Ctrl->G.file))) n_errors++;
261 break;
262 case 'I':
263 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
264 Ctrl->I.active = true;
265 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
266 break;
267 case 'N':
268 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
269 Ctrl->N.active = true;
270 Ctrl->N.t_upper = atof (opt->arg);
271 break;
272 case 'S':
273 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
274 Ctrl->S.active = true;
275 break;
276 case 'T':
277 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
278 Ctrl->T.active = true;
279 break;
280 default: /* Report bad options */
281 n_errors += gmt_default_error (GMT, opt->option);
282 break;
283 }
284 }
285
286 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = 5;
287 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
288 n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0, "Option -I: Must specify positive increment(s)\n");
289 n_errors += gmt_M_check_condition (GMT, !(Ctrl->G.active || Ctrl->G.file), "Option -G: Must specify output file\n");
290 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < 5, "Binary input data (-bi) must have at least 5 columns\n");
291
292 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
293 }
294
295 #define bailout(code) {gmt_M_free_options (mode); return (code);}
296 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
297
GMT_hotspotter(void * V_API,int mode,void * args)298 EXTERN_MSC int GMT_hotspotter (void *V_API, int mode, void *args) {
299
300 uint64_t n_smts; /* Number of seamounts read */
301 uint64_t n_track; /* Number of points along a single flowline */
302 uint64_t n_read = 0; /* Number of records read */
303 uint64_t node, node_0; /* Grid indices */
304 uint64_t n_expected_fields;
305 unsigned int n_stages; /* Number of stage rotations (poles) */
306 unsigned int row, col, kx, ky, m;
307 int node_x_width; /* Number of x-nodes covered by the seamount in question (y-dependent) */
308 int node_y_width; /* Number of y-nodes covered by the seamount */
309 int d_col, d_row, col_0, row_0, n_columns, n_rows;
310 int error = 0; /* nonzero when arguments are wrong */
311
312
313 double sampling_int_in_km; /* Sampling interval along flowline (in km) */
314 double x_smt; /* Seamount longitude (input degrees, stored as radians) */
315 double y_smt; /* Seamount latitude (input degrees, stored as radians) */
316 double z_smt; /* Seamount amplitude (in user units) */
317 double r_smt; /* Seamount radius (input km, stored as radians) */
318 double t_smt; /* Seamount upper age (up to age of seafloor) */
319 double norm; /* Normalization factor based on r_smt */
320 double *xpos = NULL, *ypos = NULL; /* Coordinates of the output grid (in radians) */
321 double *c = NULL; /* Array with one flowline */
322 double dx, dy; /* x,y distance from projected seamount center to nearest node */
323 double *latfactor = NULL, *ilatfactor = NULL; /* Arrays of latitudinal-dependent x-scales (cos(lat)-stuff) */
324 double x_part, y_part; /* Components of radius from projected seamount center to current node */
325 double y_part2; /* y_part squared */
326 double i_yinc_r; /* 1/y_inc (in radians) */
327 double r2; /* Radius squared from projected seamount center to current node */
328 double *in = NULL; /* GMT read array */
329 double wesn[4]; /* Map region in geocentric coordinates in radians */
330 double yg; /* Geodetic latitude */
331
332 char *processed_node = NULL; /* Pointer to array with true/false values for each grid node */
333
334 struct GMT_GRID *G = NULL; /* Grid structure for output CVA grid */
335 struct GMT_GRID *G_rad = NULL; /* Same but has radians in header (no grid) */
336 struct GMT_RECORD *In = NULL;
337
338 struct EULER *p = NULL; /* Array of structures with Euler stage rotations */
339 struct GMT_OPTION *ptr = NULL;
340 struct HOTSPOTTER_CTRL *Ctrl = NULL;
341 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
342 struct GMT_OPTION *options = NULL;
343 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
344
345 /*----------------------- Standard module initialization and parsing ----------------------*/
346
347 if (API == NULL) return (GMT_NOT_A_SESSION);
348 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
349 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
350
351 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
352
353 /* Parse the command-line arguments */
354
355 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 */
356 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
357 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
358 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
359 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
360
361 /*---------------------------- This is the hotspotter main code ----------------------------*/
362
363 /* Load in the Euler stage poles */
364
365 if ((error = spotter_init (GMT, Ctrl->E.file, &p, 1, false, Ctrl->E.mode, &Ctrl->N.t_upper)) < 0)
366 Return (-error);
367 n_stages = (unsigned int)error;
368
369 /* Initialize the CVA grid and structure */
370
371 if ((G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
372 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
373
374 /* Assign grid-region variables in radians to avoid conversions inside convolution loop */
375
376 if ((G_rad = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, G)) == NULL) Return (API->error);
377 G_rad->header->inc[GMT_X] = G->header->inc[GMT_X] * D2R;
378 G_rad->header->inc[GMT_Y] = G->header->inc[GMT_Y] * D2R;
379 G_rad->header->wesn[XLO] = G->header->wesn[XLO] * D2R;
380 G_rad->header->wesn[XHI] = G->header->wesn[XHI] * D2R;
381 G_rad->header->wesn[YLO] = G->header->wesn[YLO] * D2R;
382 G_rad->header->wesn[YHI] = G->header->wesn[YHI] * D2R;
383
384 /* Get geocentric wesn region in radians */
385 gmt_M_memcpy (wesn, G_rad->header->wesn, 4, double);
386 wesn[YLO] = D2R * gmt_lat_swap (GMT, R2D * wesn[YLO], GMT_LATSWAP_G2O);
387 wesn[YHI] = D2R * gmt_lat_swap (GMT, R2D * wesn[YHI], GMT_LATSWAP_G2O);
388
389 /* Initialize the CVA grid and structure */
390
391 i_yinc_r = 1.0 / G_rad->header->inc[GMT_Y];
392
393 /* Set flowline sampling interval to 1/2 of the shortest distance between x-nodes */
394
395 sampling_int_in_km = Ctrl->D.value * G_rad->header->inc[GMT_X] * EQ_RAD * ((fabs (G_rad->header->wesn[YHI]) > fabs (G_rad->header->wesn[YLO])) ? cos (G_rad->header->wesn[YHI]) : cos (G_rad->header->wesn[YLO]));
396 GMT_Report (API, GMT_MSG_INFORMATION, "Flowline sampling interval = %.3f km\n", sampling_int_in_km);
397
398 if (Ctrl->T.active) GMT_Report (API, GMT_MSG_WARNING, "Seamount ages truncated to %g\n", Ctrl->N.t_upper);
399
400 n_expected_fields = (GMT->common.b.ncol[GMT_IN]) ? GMT->common.b.ncol[GMT_IN] : 5;
401 if ((error = GMT_Set_Columns (API, GMT_IN, (unsigned int)n_expected_fields, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
402 Return (error);
403 }
404 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
405 Return (API->error);
406 }
407 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
408 Return (API->error);
409 }
410
411 /* Precalculate coordinates xpos[], ypos[] and scale factors(lat) on the grid */
412
413 xpos = gmt_grd_coord (GMT, G_rad->header, GMT_X);
414 ypos = gmt_grd_coord (GMT, G_rad->header, GMT_Y);
415
416 latfactor = gmt_M_memory (GMT, NULL, G->header->n_rows, double);
417 ilatfactor = gmt_M_memory (GMT, NULL, G->header->n_rows, double);
418
419 for (row = 0; row < G->header->n_rows; row++) {
420 latfactor[row] = G_rad->header->inc[GMT_X] * cos (ypos[row]);
421 ilatfactor[row] = 1.0 / latfactor[row];
422 }
423
424 /* Allocate T/F array */
425
426 processed_node = gmt_M_memory (GMT, NULL, G->header->size, char);
427
428 /* Start to read input data */
429
430 n_smts = 0;
431 n_columns = G->header->n_columns; n_rows = G->header->n_rows; /* Signed integers */
432
433 do { /* Keep returning records until we reach EOF */
434 n_read++;
435 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
436 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
437 gmt_M_free (GMT, processed_node);
438 Return (GMT_RUNTIME_ERROR);
439 }
440 if (gmt_M_rec_is_any_header (GMT)) /* Skip all headers */
441 continue;
442 if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
443 break;
444 }
445
446 if (In->data == NULL) {
447 gmt_quit_bad_record (API, In);
448 Return (API->error);
449 }
450
451 /* Data record to process */
452 in = In->data; /* Only need to process numerical part here */
453
454 /* STEP 1: Read information about a single seamount from input record */
455
456 if (gmt_M_is_dnan (in[4])) /* Age is NaN, assign value */
457 t_smt = Ctrl->N.t_upper;
458 else { /* Assign given value, truncate if necessary */
459 t_smt = in[4];
460 if (t_smt > Ctrl->N.t_upper) {
461 if (Ctrl->T.active)
462 t_smt = Ctrl->N.t_upper;
463 else {
464 GMT_Report (API, GMT_MSG_WARNING, "Seamounts near line %" PRIu64 " has age (%g) > oldest stage (%g) (skipped)\n",
465 n_read, t_smt, Ctrl->N.t_upper);
466 continue;
467 }
468 }
469 }
470
471 y_smt = D2R * gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O); /* Convert to geocentric, and radians */
472 x_smt = in[GMT_X] * D2R; /* Seamount positions in RADIANS */
473
474 /* STEP 2: Calculate this seamount's flowline */
475
476 if (spotter_forthtrack (GMT, &x_smt, &y_smt, &t_smt, 1, p, n_stages, sampling_int_in_km, 0.0, 0, NULL, &c) <= 0) {
477 GMT_Report (API, GMT_MSG_ERROR, "Nothing returned from spotter_forthtrack - aborting\n");
478 Return (GMT_RUNTIME_ERROR);
479 }
480
481 /* Do some normalizations here to save processing inside convolution later */
482
483 z_smt = in[GMT_Z];
484 r_smt = in[3];
485 r_smt /= EQ_RAD; /* Converts radius in km to radians */
486 norm = -4.5 / (r_smt * r_smt); /* Gaussian normalization */
487 node_y_width = irint (ceil (i_yinc_r * r_smt)); /* y-node coverage */
488
489 /* STEP 3: Convolve this flowline with seamount shape and add to CVA grid */
490
491 /* Our convolution is approximate: We sample the flowline frequently and use
492 * one of the points on the flowline that are closest to the node. Ideally we
493 * want the nearest distance from each node to the flowline. Later versions may
494 * improve on this situation */
495
496 n_track = lrint (c[0]); /* Number of point pairs making up this flowline */
497
498 gmt_M_memset (processed_node, G->header->size, char); /* Fresh start for this flowline convolution */
499
500 for (m = 0, kx = 1; m < n_track; m++, kx += 2) { /* For each point along flowline */
501
502 ky = kx + 1; /* Index for the y-coordinate */
503
504 /* First throw out points outside specified grid region */
505
506 if (c[ky] < wesn[YLO] || c[ky] > wesn[YHI]) continue; /* Geocentric latitude outside region */
507
508 if (c[kx] > wesn[XHI]) c[kx] -= TWO_PI;
509 while (c[kx] < wesn[XLO]) c[kx] += TWO_PI;
510 if (c[kx] > wesn[XHI]) continue; /* Longitude outside region */
511
512 /* OK, this point is within our region, get node index */
513
514 col = (unsigned int)gmt_M_grd_x_to_col (GMT, c[kx], G_rad->header);
515 yg = gmt_lat_swap (GMT, R2D * c[ky], GMT_LATSWAP_O2G); /* Convert back to geodetic */
516 row = (unsigned int)gmt_M_grd_y_to_row (GMT, yg, G->header);
517 node = gmt_M_ijp (G->header, row, col);
518
519 if (!processed_node[node]) { /* Have not added to the CVA at this node yet */
520
521 /* Shape is z_smt * exp (r^2 * norm) */
522
523 node_x_width = irint (ceil (r_smt * ilatfactor[row]));
524 dx = c[kx] - xpos[col];
525 dy = c[ky] - ypos[row];
526
527 /* Loop over a square that circumscribes this seamount's basal outline */
528
529 for (d_row = -node_y_width, row_0 = row - node_y_width; d_row <= node_y_width; d_row++, row_0++) {
530
531 if (row_0 < 0 || row_0 >= n_rows) continue; /* Outside grid */
532
533 y_part = d_row * G_rad->header->inc[GMT_Y] - dy;
534 y_part2 = y_part * y_part;
535 node_0 = gmt_M_ijp (G->header, row_0, 0);
536
537 for (d_col = -node_x_width, col_0 = col - node_x_width; d_col <= node_x_width; d_col++, col_0++) {
538
539 if (col_0 < 0 || col_0 >= n_columns) continue; /* Outside grid */
540
541 x_part = d_col * latfactor[row] - dx;
542 r2 = (x_part * x_part + y_part2) * norm;
543 G->data[node_0+col_0] += (gmt_grdfloat)(z_smt * exp (r2));
544 }
545 }
546 processed_node[node] = 1; /* Now we have visited this node */
547 }
548 }
549
550 gmt_M_free (GMT, c); /* Free the flowline vector */
551
552 n_smts++; /* Go to next seamount */
553
554 if (!(n_smts%100)) GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5ld seamounts\r", n_smts);
555 } while (true);
556
557 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
558 gmt_M_free (GMT, processed_node); gmt_M_free (GMT, latfactor); gmt_M_free (GMT, p);
559 gmt_M_free (GMT, ilatfactor); gmt_M_free (GMT, xpos); gmt_M_free (GMT, ypos);
560 Return (API->error);
561 }
562
563 /* OK, Done processing, time to write out */
564
565 GMT_Report (API, GMT_MSG_INFORMATION, "Processed %5ld seamounts\n", n_smts);
566
567 if (Ctrl->S.active) { /* Convert CVA values to percent of CVA maximum */
568 uint64_t node;
569 double scale;
570
571 GMT_Report (API, GMT_MSG_INFORMATION, "Normalize CVS grid to percentages of max CVA\n");
572 G->header->z_min = +DBL_MAX;
573 G->header->z_max = -DBL_MAX;
574 gmt_M_grd_loop (GMT, G, row, col, node) { /* Loop over all output nodes */
575 if (gmt_M_is_fnan (G->data[node])) continue;
576 if (G->data[node] < G->header->z_min) G->header->z_min = G->data[node];
577 if (G->data[node] > G->header->z_max) G->header->z_max = G->data[node];
578 }
579 GMT_Report (API, GMT_MSG_INFORMATION, "CVA min/max: %g %g -> ", G->header->z_min, G->header->z_max);
580 scale = 100.0 / G->header->z_max;
581 for (node = 0; node < G->header->size; node++) G->data[node] *= (gmt_grdfloat)scale;
582 G->header->z_min *= scale; G->header->z_max *= scale;
583 GMT_Report (API, GMT_MSG_INFORMATION, "%g %g\n", G->header->z_min, G->header->z_max);
584 }
585 GMT_Report (API, GMT_MSG_INFORMATION, "Write CVA grid %s\n", Ctrl->G.file);
586
587 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G)) Return (API->error);
588 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, G) != GMT_NOERROR) {
589 gmt_M_free (GMT, processed_node); gmt_M_free (GMT, latfactor); gmt_M_free (GMT, p);
590 gmt_M_free (GMT, ilatfactor); gmt_M_free (GMT, xpos); gmt_M_free (GMT, ypos);
591 Return (API->error);
592 }
593
594 /* Clean up memory */
595
596 gmt_M_free (GMT, processed_node); gmt_M_free (GMT, latfactor); gmt_M_free (GMT, p);
597 gmt_M_free (GMT, ilatfactor); gmt_M_free (GMT, xpos); gmt_M_free (GMT, ypos);
598 if (GMT_Destroy_Data (API, &G_rad) != GMT_NOERROR) {
599 GMT_Report (API, GMT_MSG_ERROR, "Failed to free G_rad\n");
600 }
601
602 Return (GMT_NOERROR);
603 }
604