1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 2000-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 * originater reads file of seamount locations and tries to match each
18 * seamount with a probable hotspot by drawing flowines back in time and
19 * keeping track of which hotspot is closest to each flowline. It then
20 * reports the closest hotspot, the stage of the flowline involved, the
21 * implied pseudo-age of the seamount, and the minimum distance between
22 * the flowline and hotspot (in km).
23 *
24 * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
25 * Date: 29-DEC-1999
26 * Version: 1.0
27 *
28 *-------------------------------------------------------------------------
29 * An ASCII stage pole (Euler) file must have following format:
30 *
31 * 1. Any number of comment lines starting with # in first column
32 * 2. Any number of blank lines (just carriage return, no spaces)
33 * 2. Any number of stage pole records which each have the format:
34 * lon(deg) lat(deg) tstart(Ma) tstop(Ma) ccw-angle(deg)
35 * 3. stage records must go from oldest to youngest rotation
36 * 4. Note tstart is larger (older) that tstop for each record
37 * 5. No gaps allowed: tstart must equal the previous records tstop
38 *
39 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
40 *
41 * # Time in Ma, angles in degrees
42 * # lon lat tstart tend ccw-angle
43 * 165 85 150 100 24.0
44 * 284 36 100 74 15.0
45 * 265 22 74 65 7.5
46 * 253 17 65 42 14.0
47 * 285 68 42 0 34.0
48 *
49 * AN ASCII total reconstruction file must have the following format:
50 *
51 * 1. Any number of comment lines starting with # in first column
52 * 2. Any number of blank lines (just carriage return, no spaces)
53 * 2. Any number of finite pole records which each have the format:
54 * lon(deg) lat(deg) tstop(Ma) ccw-angle(deg)
55 * 3. total reconstruction rotations must go from youngest to oldest
56 *
57 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
58 *
59 * # Time in Ma, angles in degrees
60 * #longitude latitude time(My) angle(deg)
61 * 285.00000 68.00000 42.0000 34.0000
62 * 275.66205 53.05082 65.0000 43.5361
63 * 276.02501 48.34232 74.0000 50.0405
64 * 279.86436 46.30610 100.0000 64.7066
65 * 265.37800 55.69932 150.0000 82.9957
66 *
67 * ASCII seamount location file(s) must have the following format:
68 *
69 * 1. Any number of comment lines starting with # in first column
70 * 2. Any number of blank lines (just carriage return, no spaces)
71 * 3. For special header records, see -H
72 * 4. Any number of data records which each have the format:
73 * lon lat height radius crustal_age (or lat lon ..., see -: option).
74 * crustal_age in Ma, height and radius are not used by originater but
75 * are used by hotspotter.
76 *
77 * Binary files cannot have header records, and data fields must all be
78 * either single or double precision (see -bi option). Output file will
79 * be ASCII since it contains a text string (hotspot ID).
80 *
81 * The file with a list of hotspots must have the following format:
82 *
83 * 1. Any number of comment lines starting with # in first column
84 * 2. Any number of blank lines (just carriage return, no spaces)
85 * 2. Any number of hotspot records which each have the format:
86 * lon(deg) lat(deg) id name
87 * the id is a 3-character tag (e.g., HWI), the name is the
88 * full name of the hotspot (e.g., Hawaii).
89 *
90 * Example:
91 *
92 * # Partial list (Pacific) of HotSpots from Table 1 of Yamaji, 1992
93 * #Lon Lat Abbreviation Hotspot_name
94 * 167 3 CRL Caroline
95 * 230 46 COB Cobb
96 * 205 20 HWI Hawaii
97 * 221.9 -50.9 LSV Louisville
98 * 220 -29 MDN MacDonald
99 * 221 -11 MRQ Marquesas
100 * 231 -27 PTC Pitcairn
101 * 254 -27 SLG Sala y Gomez
102 * 192 -15 SAM Samoa
103 * 212 -18 SOC Society
104 */
105
106 #include "gmt_dev.h"
107 #include "spotter.h"
108
109 #define THIS_MODULE_CLASSIC_NAME "originater"
110 #define THIS_MODULE_MODERN_NAME "originater"
111 #define THIS_MODULE_LIB "spotter"
112 #define THIS_MODULE_PURPOSE "Associate seamounts with nearest hotspot point sources"
113 #define THIS_MODULE_KEYS "<D{,FD(=,>D}"
114 #define THIS_MODULE_NEEDS ""
115 #define THIS_MODULE_OPTIONS "-:>Vbdehiqs" GMT_OPT("HMm")
116
117 EXTERN_MSC int gmtlib_great_circle_intersection (struct GMT_CTRL *GMT, double A[], double B[], double C[], double X[], double *CX_dist);
118 EXTERN_MSC double gmtlib_great_circle_dist_degree (struct GMT_CTRL *GMT, double lon1, double lat1, double lon2, double lat2);
119
120 #define KM_PR_RAD (R2D * GMT->current.proj.DIST_KM_PR_DEG)
121
122 struct HOTSPOT_ORIGINATOR {
123 struct HOTSPOT *h; /* Pointer to regular HOTSPOT structure */
124 /* Extra variables needed for this program */
125 double np_dist; /* Distance to nearest point on the current flowline */
126 double np_sign; /* "Sign" of this distance (see code) */
127 double np_time; /* Predicted time at nearest point */
128 double np_lon; /* Longitude of nearest point on the current flowline */
129 double np_lat; /* Latitude of nearest point on the current flowline */
130 uint64_t nearest; /* Point id of current flowline node points closest to hotspot */
131 unsigned int stage; /* Stage to which seamount belongs */
132 struct GMT_DATASEGMENT *D; /* Used for drifting hotspots only */
133 };
134
135 struct ORIGINATOR_CTRL { /* All control options for this program (except common args) */
136 /* active is true if the option has been activated */
137 struct ORIGINATER_D { /* -D<factor */
138 bool active;
139 double value;
140 } D;
141 struct ORIGINATER_E { /* -Erotfile[+i] */
142 bool active;
143 bool mode;
144 char *file;
145 } E;
146 struct ORIGINATER_F { /* -Fhotspotfile[+d] */
147 bool active;
148 bool mode;
149 char *file;
150 } F;
151 struct ORIGINATER_L { /* -L */
152 bool active;
153 unsigned int mode;
154 bool degree; /* Report degrees */
155 } L;
156 struct ORIGINATER_N { /* -N */
157 bool active;
158 double t_upper;
159 } N;
160 struct ORIGINATER_Q { /* -Q<tfix> */
161 bool active;
162 double t_fix, r_fix;
163 } Q;
164 struct ORIGINATER_S { /* -S */
165 bool active;
166 unsigned int n;
167 } S;
168 struct ORIGINATER_T { /* -T */
169 bool active;
170 } T;
171 struct ORIGINATER_W { /* -W<max_dist> */
172 bool active;
173 double dist;
174 } W;
175 struct ORIGINATER_Z { /* -Z */
176 bool active;
177 } Z;
178 };
179
New_Ctrl(struct GMT_CTRL * GMT)180 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
181 struct ORIGINATOR_CTRL *C;
182
183 C = gmt_M_memory (GMT, NULL, 1, struct ORIGINATOR_CTRL);
184
185 /* Initialize values whose defaults are not 0/false/NULL */
186
187 C->D.value = 5.0;
188 C->N.t_upper = 180.0;
189 C->S.n = 1;
190 C->W.dist = 1.0e100;
191 return (C);
192 }
193
Free_Ctrl(struct GMT_CTRL * GMT,struct ORIGINATOR_CTRL * C)194 static void Free_Ctrl (struct GMT_CTRL *GMT, struct ORIGINATOR_CTRL *C) { /* Deallocate control structure */
195 if (!C) return;
196 gmt_M_str_free (C->E.file);
197 gmt_M_str_free (C->F.file);
198 gmt_M_free (GMT, C);
199 }
200
originater_comp_hs(const void * p1,const void * p2)201 GMT_LOCAL int originater_comp_hs (const void *p1, const void *p2) {
202 const struct HOTSPOT_ORIGINATOR *a = p1, *b = p2;
203
204 if (a->np_dist < b->np_dist) return (-1);
205 if (a->np_dist > b->np_dist) return (+1);
206 return (0);
207 }
208
usage(struct GMTAPI_CTRL * API,int level)209 static int usage (struct GMTAPI_CTRL *API, int level) {
210 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
211 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
212 GMT_Usage (API, 0, "usage: %s [<table>] %s -F<hotspottable>[+d] [-D<d_km>] [-H] [-L[l|t|w]] "
213 "[-N<upper_age>] [-Qr/t] [-S<n_hs>] [-T] [%s] [-W<maxdist>] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
214 name, SPOTTER_E_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_d_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
215
216 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
217
218 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
219 GMT_Usage (API, 1, "\n<table> (in ASCII, binary, or netCDF) has 5 or more columns. If no file(s) is given, "
220 "standard input is read. Expects (x,y,z,r,t) records, with t in Ma.");
221 spotter_rot_usage (API);
222 GMT_Usage (API, 1, "\n-F<hotspottable>[+d]");
223 GMT_Usage (API, -2, "Specify file name for hotspot locations. "
224 "Append +d if we should look for hotspot drift tables. "
225 "If found then we interpolate to get hotspot location as a function of time [fixed].");
226 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
227 GMT_Usage (API, 1, "\n-D<d_km>");
228 GMT_Usage (API, -2, "Set sampling interval in km along tracks [5].");
229 GMT_Usage (API, 1, "\n-L[l|t|w]");
230 GMT_Usage (API, -2, "Output information for closest approach for nearest hotspot only (ignores -S). Select directive:");
231 GMT_Usage (API, 3, "l: Give (lon, lat, time, dist, z).");
232 GMT_Usage (API, 3, "t: Give (time, dist, z) [Default].");
233 GMT_Usage (API, 3, "w: Give (omega, dist, z).");
234 GMT_Usage (API, -2, "Note: dist is in km; use upper case L,T,W to get dist in spherical degrees.");
235 GMT_Usage (API, 1, "\n-N<upper_age>");
236 GMT_Usage (API, -2, "Set age (in m.y.) for seafloor where age == NaN [180].");
237 GMT_Usage (API, 1, "\n-Qr/t");
238 GMT_Usage (API, -2, "Input files has (x,y,z) only. Append constant r/t to append to input record.");
239 GMT_Usage (API, 1, "\n-S<n_hs>");
240 GMT_Usage (API, -2, "Report the <n_hs> closest hotSpots [1].");
241 GMT_Usage (API, 1, "\n-T Truncate seamount ages exceeding the upper age set with -N [no truncation].");
242 GMT_Option (API, "V");
243 GMT_Usage (API, 1, "\n-W<maxdist>");
244 GMT_Usage (API, -2, "Report seamounts whose closest encounter to a hotspot is less than <maxdist> km "
245 "[Default reports for all seamounts].");
246 GMT_Usage (API, 1, "\n-Z Write hotspot ID number rather than hotspot TAG.");
247 GMT_Option (API, "bi5,d,e,h,i,q,s,:,.");
248
249 return (GMT_MODULE_USAGE);
250 }
251
parse(struct GMT_CTRL * GMT,struct ORIGINATOR_CTRL * Ctrl,struct GMT_OPTION * options)252 static int parse (struct GMT_CTRL *GMT, struct ORIGINATOR_CTRL *Ctrl, struct GMT_OPTION *options) {
253 /* This parses the options provided to originater and sets parameters in CTRL.
254 * Any GMT common options will override values set previously by other commands.
255 * It also replaces any file names specified as input or output with the data ID
256 * returned when registering these sources/destinations with the API.
257 */
258
259 unsigned int n_errors = 0, n_input;
260 int k;
261 char *c = NULL;
262 struct GMT_OPTION *opt = NULL;
263 struct GMTAPI_CTRL *API = GMT->parent;
264
265 for (opt = options; opt; opt = opt->next) {
266 switch (opt->option) {
267
268 case '<': /* Skip input files */
269 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
270 break;
271
272 /* Supplemental parameters */
273 case 'C': /* Now done automatically in spotter_init */
274 if (gmt_M_compat_check (GMT, 4))
275 GMT_Report (API, GMT_MSG_COMPAT, "-C is no longer needed as total reconstruction vs stage rotation is detected automatically.\n");
276 else
277 n_errors += gmt_default_error (GMT, opt->option);
278 break;
279 case 'D':
280 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
281 Ctrl->D.active = true;
282 Ctrl->D.value = atof (opt->arg);
283 break;
284 case 'E':
285 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
286 Ctrl->E.active = true; k = 0;
287 if (opt->arg[0] == '+') { Ctrl->E.mode = true; k = 1;}
288 else if ((c = strstr (opt->arg, "+i"))) {Ctrl->E.mode = true; c[0] = '\0';}
289 if (opt->arg[k]) Ctrl->E.file = strdup (&opt->arg[k]);
290 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->E.file))) n_errors++;
291 if (c) c[0] = '+';
292 break;
293 case 'F':
294 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
295 Ctrl->F.active = true; k = 0;
296 if (opt->arg[0] == '+') { Ctrl->F.mode = true; k = 1;}
297 if (opt->arg[k]) Ctrl->F.file = strdup (&opt->arg[k]);
298 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
299 break;
300 case 'L':
301 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
302 Ctrl->L.active = true;
303 switch (opt->arg[0]) {
304 case 'L':
305 Ctrl->L.degree = true;
306 /* Intentionally fall through - to 'l' */
307 case 'l':
308 Ctrl->L.mode = 3;
309 break;
310 case 'w':
311 case 'W':
312 Ctrl->L.mode = 2;
313 break;
314 case 'T':
315 Ctrl->L.degree = true;
316 /* Intentionally fall through - to 't' */
317 case 't':
318 Ctrl->L.mode = 1;
319 break;
320 default:
321 Ctrl->L.mode = 1;
322 break;
323 }
324 break;
325 case 'N':
326 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
327 Ctrl->N.active = true;
328 Ctrl->N.t_upper = atof (opt->arg);
329 break;
330 case 'Q':
331 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
332 Ctrl->Q.active = true;
333 sscanf (opt->arg, "%lg/%lg", &Ctrl->Q.r_fix, &Ctrl->Q.t_fix);
334 break;
335 case 'S':
336 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
337 Ctrl->S.active = true;
338 k = atoi (opt->arg);
339 n_errors += gmt_M_check_condition (GMT, k < 1, "Option -S: Must specify a positive number of hotspots\n");
340 Ctrl->S.n = k;
341 break;
342 case 'T':
343 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
344 Ctrl->T.active = true;
345 break;
346 case 'W':
347 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
348 Ctrl->W.active = true;
349 Ctrl->W.dist = atof (opt->arg);
350 break;
351 case 'Z':
352 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
353 Ctrl->Z.active = true;
354 break;
355 default: /* Report bad options */
356 n_errors += gmt_default_error (GMT, opt->option);
357 break;
358 }
359 }
360
361 n_input = (Ctrl->Q.active) ? 3 : 5;
362 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = n_input;
363
364 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < n_input, "Binary input data (-bi) must have at least %d columns\n", n_input);
365 n_errors += gmt_M_check_condition (GMT, !Ctrl->F.file, "Option -F: Must specify hotspot file\n");
366 n_errors += gmt_M_check_condition (GMT, !Ctrl->E.file, "Option -F: Must specify Euler pole file\n");
367 n_errors += gmt_M_check_condition (GMT, Ctrl->D.value <= 0.0, "Option -D: Must specify a positive interval\n");
368 n_errors += gmt_M_check_condition (GMT, Ctrl->W.dist <= 0.0, "Option -W: Must specify a positive distance in km\n");
369
370 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
371 }
372
373 #define bailout(code) {gmt_M_free_options (mode); return (code);}
374 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
375
GMT_originater(void * V_API,int mode,void * args)376 EXTERN_MSC int GMT_originater (void *V_API, int mode, void *args) {
377 unsigned int n_max_spots, n_input;
378 unsigned int spot, smt, n_stages, n_hotspots, n_read, n_skipped = 0;
379 uint64_t k, kk, np, n_expected_fields, n_out;
380
381 int error = 0, ns;
382 bool better;
383
384 double x_smt, y_smt, z_smt, r_smt, t_smt, t, *c = NULL, *in = NULL, dist, dlon, lon, lat, out[5];
385 double hx_dist, hx_dist_km, dist_NA, dist_NX, del_dist, dt = 0.0, A[3], H[3], N[3], X[3];
386
387 char record[GMT_BUFSIZ] = {""}, buffer[GMT_BUFSIZ] = {""}, fmt1[GMT_BUFSIZ] = {""}, fmt2[GMT_BUFSIZ] = {""};
388
389 struct EULER *p = NULL;
390 struct HOTSPOT *orig_hotspot = NULL;
391 struct HOTSPOT_ORIGINATOR *hotspot = NULL, *hot = NULL;
392 struct GMT_DATASET *F = NULL;
393 struct GMT_RECORD *In = NULL, *Out = NULL;
394 struct GMT_OPTION *ptr = NULL;
395 struct ORIGINATOR_CTRL *Ctrl = NULL;
396 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
397 struct GMT_OPTION *options = NULL;
398 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
399
400 /*----------------------- Standard module initialization and parsing ----------------------*/
401
402 if (API == NULL) return (GMT_NOT_A_SESSION);
403 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
404 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
405
406 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
407
408 /* Parse the command-line arguments */
409
410 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 */
411 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
412 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
413 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
414 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
415
416 /*---------------------------- This is the originater main code ----------------------------*/
417
418 if ((ns = spotter_hotspot_init (GMT, Ctrl->F.file, true, &orig_hotspot)) < 0) { /* Get geocentric hotspot locations */
419 Return (GMT_RUNTIME_ERROR); /* An error message was already issued by spotter_hotspot_init() */
420 }
421
422 n_hotspots = (unsigned int)ns;
423 if (Ctrl->S.n > n_hotspots) {
424 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Give value between 1 and %d\n", n_hotspots);
425 Return (GMT_RUNTIME_ERROR);
426 }
427 n_max_spots = MIN (Ctrl->S.n, n_hotspots);
428
429 hotspot = gmt_M_memory (GMT, NULL, n_hotspots, struct HOTSPOT_ORIGINATOR);
430 for (spot = 0; spot < n_hotspots; spot++) {
431 hotspot[spot].h = &orig_hotspot[spot]; /* Point to the original hotspot structures */
432 hotspot[spot].np_dist = 1.0e100;
433 if (Ctrl->F.mode) { /* See if there is a drift file for this hotspot */
434 char path[PATH_MAX] = {""}, file[GMT_LEN64] = {""};
435 uint64_t row;
436 sprintf (file, "%s_drift.txt", hotspot[spot].h->abbrev);
437 strncpy (path, file, PATH_MAX);
438 if (gmt_access (GMT, path, R_OK)) { /* Not found in current dir or GMT_DATADIR; check if -F gave an explicit directory */
439 if (strchr (Ctrl->F.file, '/')) { /* Filename has leading path so we will use that path */
440 strncpy (path, Ctrl->F.file, PATH_MAX);
441 k = strlen (path);
442 while (k && path[k] != '/') k--; /* Look for last slash */
443 k++; path[k] = 0; /* Truncate anything after last slash */
444 strcat (path, file); /* Prepend path to drift file name */
445 if (gmt_access (GMT, path, R_OK)) continue; /* file do not exist there either */
446 }
447 else /* No directory given, so nothing to do */
448 continue;
449 }
450 /* Here we successfully found a drift file, somewhere */
451 if ((F = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, file, NULL)) == NULL) {
452 Return (API->error);
453 }
454 if (F->n_columns < 3) {
455 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 3 are needed\n", file, (int)F->n_columns);
456 Return (GMT_DIM_TOO_SMALL);
457 }
458 hotspot[spot].D = F->table[0]->segment[0]; /* Only one table with one segment for histories */
459 for (row = 0; row < hotspot[spot].D->n_rows; row++) hotspot[spot].D->data[GMT_Y][row] = gmt_lat_swap (GMT, hotspot[spot].D->data[GMT_Y][row], GMT_LATSWAP_G2O); /* Convert to geocentric */
460 }
461 }
462
463 if ((error = spotter_init (GMT, Ctrl->E.file, &p, 1, false, Ctrl->E.mode, &Ctrl->N.t_upper)) < 0)
464 Return (-error);
465 n_stages = (unsigned int)error;
466
467 hot = gmt_M_memory (GMT, NULL, n_hotspots, struct HOTSPOT_ORIGINATOR);
468
469 sprintf (fmt1, "%s%s%s%s%s%s%s%s%%s", GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out,
470 GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator,
471 GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator);
472 if (Ctrl->Z.active)
473 sprintf (fmt2, "%s%%d%s%%d%s%s%s%s", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out);
474 else
475 sprintf (fmt2, "%s%%s%s%%d%s%s%s%s", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out, GMT->current.setting.io_col_separator, GMT->current.setting.format_float_out);
476 n_input = (Ctrl->Q.active) ? 3 : 5;
477 n_expected_fields = (GMT->common.b.ncol[GMT_IN]) ? GMT->common.b.ncol[GMT_IN] : n_input;
478 if (Ctrl->L.active) {
479 n_out = (Ctrl->L.mode == 3) ? 5 : 3;
480 }
481 else
482 n_out = n_expected_fields;
483 if (n_out == 3)
484 gmt_set_cartesian (GMT, GMT_OUT); /* Since output is no longer lon/lat */
485 if ((error = GMT_Set_Columns (API, GMT_IN, (unsigned int)n_out, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
486 Return (error);
487 }
488 if ((error = GMT_Set_Columns (API, GMT_OUT, (unsigned int)n_expected_fields, GMT_COL_FIX)) != GMT_NOERROR) {
489 Return (error);
490 }
491 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
492 Return (API->error);
493 }
494 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
495 Return (API->error);
496 }
497 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
498 Return (API->error);
499 }
500 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
501 Return (API->error);
502 }
503 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
504 Return (API->error);
505 }
506
507 n_read = smt = 0;
508 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
509
510 do { /* Keep returning records until we reach EOF */
511 n_read++;
512 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
513 if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */
514 Return (GMT_RUNTIME_ERROR);
515 if (gmt_M_rec_is_any_header (GMT)) /* Skip all headers */
516 continue;
517 if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
518 break;
519 }
520
521 if (In->data == NULL) {
522 gmt_quit_bad_record (API, In);
523 Return (API->error);
524 }
525
526 /* Data record to process */
527 in = In->data; /* Only need to process numerical part here */
528 Out->text = In->text;
529
530 if (n_input == 3) { /* Set constant r,t values based on -Q */
531 in[3] = Ctrl->Q.r_fix;
532 in[4] = Ctrl->Q.t_fix;
533 }
534 if (gmt_M_is_dnan (in[4])) /* Age is NaN, assign upper value */
535 t_smt = Ctrl->N.t_upper;
536 else { /* Assign given value, truncate if necessary */
537 t_smt = in[4];
538 if (t_smt > Ctrl->N.t_upper) {
539 if (Ctrl->T.active)
540 t_smt = Ctrl->N.t_upper;
541 else {
542 GMT_Report (API, GMT_MSG_WARNING, "Seamounts near line %d has age (%g) > oldest stage (%g) (skipped)\n", n_read, t_smt, Ctrl->N.t_upper);
543 continue;
544 }
545 }
546 }
547 if (t_smt < 0.0) { /* Negative ages are flags for points to be skipped */
548 n_skipped++;
549 continue;
550 }
551 y_smt = D2R * gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O); /* Convert to geocentric, and radians */
552 x_smt = in[GMT_X] * D2R;
553 z_smt = in[GMT_Z];
554 r_smt = in[3];
555
556 if (!(smt % 10)) GMT_Report (API, GMT_MSG_INFORMATION, "Working on seamount # %5d\r", smt);
557
558 if (spotter_forthtrack (GMT, &x_smt, &y_smt, &t_smt, 1, p, n_stages, Ctrl->D.value, 0.0, 1, NULL, &c) <= 0) {
559 GMT_Report (API, GMT_MSG_ERROR, "Nothing returned from spotter_forthtrack - aborting\n");
560 Return (GMT_RUNTIME_ERROR);
561 }
562
563 np = lrint (c[0]);
564
565 gmt_M_memcpy (hot, hotspot, n_hotspots, struct HOTSPOT_ORIGINATOR);
566
567 for (kk = 0, k = 1; kk < np; kk++, k += 3) { /* For this seamounts track */
568 for (spot = 0; spot < n_hotspots; spot++) { /* For all hotspots */
569 if (hot[spot].D) { /* Must interpolate drifting hotspot location at current time c[k+2] */
570 t = c[k+2]; /* Current time */
571 gmt_intpol (GMT, hot[spot].D->data[GMT_Z], hot[spot].D->data[GMT_X], NULL, hot[spot].D->n_rows, 1, &t, &lon, 0.0, GMT->current.setting.interpolant);
572 gmt_intpol (GMT, hot[spot].D->data[GMT_Z], hot[spot].D->data[GMT_Y], NULL, hot[spot].D->n_rows, 1, &t, &lat, 0.0, GMT->current.setting.interpolant);
573 }
574 else { /* Use the fixed hotspot location */
575 lon = hot[spot].h->lon;
576 lat = hot[spot].h->lat;
577 }
578 /* Compute distance from track location to (moving or fixed) hotspot */
579 dist = gmtlib_great_circle_dist_degree (GMT, lon, lat, R2D * c[k], R2D * c[k+1]);
580 if (!Ctrl->L.degree) dist *= GMT->current.proj.DIST_KM_PR_DEG;
581 if (dist < hot[spot].np_dist) {
582 hot[spot].np_dist = dist;
583 hot[spot].nearest = kk; /* Index of nearest point on the flowline */
584 }
585 }
586 }
587 for (spot = 0; spot < n_hotspots; spot++) {
588
589 if (hot[spot].D) { /* Must interpolate drifting hotspot location at current time c[k+2] */
590 t = c[3*hot[spot].nearest+3]; /* Time of closest approach */
591 gmt_intpol (GMT, hot[spot].D->data[GMT_Z], hot[spot].D->data[GMT_X], NULL, hot[spot].D->n_rows, 1, &t, &lon, 0.0, GMT->current.setting.interpolant);
592 gmt_intpol (GMT, hot[spot].D->data[GMT_Z], hot[spot].D->data[GMT_Y], NULL, hot[spot].D->n_rows, 1, &t, &lat, 0.0, GMT->current.setting.interpolant);
593 }
594 else { /* Use the fixed hotspot location */
595 lon = hot[spot].h->lon;
596 lat = hot[spot].h->lat;
597 }
598 gmt_geo_to_cart (GMT, lat, lon, H, true); /* 3-D Cartesian vector of this hotspot */
599
600 /* Fine-tune the nearest point by considering intermediate points along greatcircle between knot points */
601
602 k = 3 * hot[spot].nearest + 1; /* Corresponding index for x into the (x,y,t) array c */
603 gmt_geo_to_cart (GMT, c[k+1], c[k], N, false); /* 3-D vector of nearest node to this hotspot */
604 better = false;
605 if (hot[spot].nearest > 0) { /* There is a point along the flowline before the nearest node */
606 gmt_geo_to_cart (GMT, c[k-2], c[k-3], A, false); /* 3-D vector of end of this segment */
607 if (gmtlib_great_circle_intersection (GMT, A, N, H, X, &hx_dist) == 0) { /* X is between A and N */
608 hx_dist_km = d_acos (hx_dist) * KM_PR_RAD;
609 if (hx_dist_km < hot[spot].np_dist) { /* This intermediate point is even closer */
610 gmt_cart_to_geo (GMT, &hot[spot].np_lat, &hot[spot].np_lon, X, true);
611 hot[spot].np_dist = hx_dist_km;
612 dist_NA = d_acos (fabs (gmt_dot3v (GMT, A, N))) * KM_PR_RAD;
613 dist_NX = d_acos (fabs (gmt_dot3v (GMT, X, N))) * KM_PR_RAD;
614 del_dist = dist_NA - dist_NX;
615 dt = (del_dist > 0.0) ? (c[k+2] - c[k-1]) * dist_NX / del_dist : 0.0;
616 better = true;
617 }
618 }
619 }
620 if (hot[spot].nearest < (np-1) ) { /* There is a point along the flowline after the nearest node */
621 gmt_geo_to_cart (GMT, c[k+4], c[k+3], A, false); /* 3-D vector of end of this segment */
622 if (gmtlib_great_circle_intersection (GMT, A, N, H, X, &hx_dist) == 0) { /* X is between A and N */
623 hx_dist_km = d_acos (hx_dist) * KM_PR_RAD;
624 if (hx_dist_km < hot[spot].np_dist) { /* This intermediate point is even closer */
625 gmt_cart_to_geo (GMT, &hot[spot].np_lat, &hot[spot].np_lon, X, true);
626 hot[spot].np_dist = hx_dist_km;
627 dist_NA = d_acos (fabs (gmt_dot3v (GMT, A, N))) * KM_PR_RAD;
628 dist_NX = d_acos (fabs (gmt_dot3v (GMT, X, N))) * KM_PR_RAD;
629 del_dist = dist_NA - dist_NX;
630 dt = (del_dist > 0.0) ? (c[k+5] - c[k+2]) * dist_NX / del_dist : 0.0;
631 better = true;
632 }
633 }
634 }
635 if (better) { /* Point closer to hotspot was found between nodes */
636 hot[spot].np_time = c[k+2] + dt; /* Add time adjustment */
637 }
638 else { /* Just use node coordinates */
639 hot[spot].np_lon = c[k] * R2D; /* Longitude of the flowline's closest approach to hotspot */
640 hot[spot].np_lat = c[k+1] * R2D; /* Latitude of the flowline's closest approach to hotspot */
641 hot[spot].np_time = c[k+2]; /* Predicted time at the flowline's closest approach to hotspot */
642 }
643
644 /* Assign sign to distance: If the vector from the hotspot pointing up along the trail is positive
645 * x-axis and y-axis is normal to that, flowlines whose closest approach point's longitude is
646 * further east are said to have negative distance. */
647
648 gmt_M_set_delta_lon (hot[spot].np_lon, lon, dlon);
649 hot[spot].np_sign = copysign (1.0, dlon);
650
651 /* Assign stage id for this point on the flowline */
652
653 k = 0;
654 while (k < n_stages && hot[spot].np_time <= p[k].t_stop) k++;
655 hot[spot].stage = n_stages - (unsigned int)k;
656 if (hot[spot].stage == 0) hot[spot].stage++;
657 }
658
659 if (n_hotspots > 1) qsort (hot, n_hotspots, sizeof (struct HOTSPOT_ORIGINATOR), originater_comp_hs);
660
661 if (hot[0].np_dist < Ctrl->W.dist) {
662 if (Ctrl->L.mode == 1) { /* Want time, dist, z output */
663 out[0] = hot[0].np_time;
664 out[1] = hot[0].np_dist * hot[0].np_sign;
665 out[2] = z_smt;
666 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
667 }
668 else if (Ctrl->L.mode == 2) { /* Want omega, dist, z output */
669 out[0] = spotter_t2w (GMT, p, n_stages, hot[0].np_time);
670 out[1] = hot[0].np_dist * hot[0].np_sign;
671 out[2] = z_smt;
672 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
673 }
674 else if (Ctrl->L.mode == 3) { /* Want x, y, time, dist, z output */
675 out[GMT_X] = hot[0].np_lon;
676 out[GMT_Y] = gmt_lat_swap (GMT, hot[0].np_lat, GMT_LATSWAP_O2G); /* Convert back to geodetic */
677 out[2] = hot[0].np_time;
678 out[3] = hot[0].np_dist * hot[0].np_sign;
679 out[4] = z_smt;
680 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
681 }
682 else { /* Conventional originater output */
683 out[GMT_X] = in[GMT_X]; out[GMT_Y] = in[GMT_Y]; out[GMT_Z] = z_smt;
684 out[3] = r_smt; out[4] = (t_smt == 180.0) ? GMT->session.d_NaN : t_smt;
685 record[0] = '\0';
686 for (spot = 0; spot < n_max_spots; spot++) {
687 if (Ctrl->Z.active)
688 sprintf (buffer, fmt2, hot[spot].h->id, hot[spot].stage, hot[spot].np_time, hot[spot].np_dist);
689 else
690 sprintf (buffer, fmt2, hot[spot].h->abbrev, hot[spot].stage, hot[spot].np_time, hot[spot].np_dist);
691 strcat (record, buffer);
692 }
693 Out->text = record;
694 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
695 }
696 }
697
698 gmt_M_free (GMT, c);
699 smt++;
700 } while (true);
701
702 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
703 Return (API->error);
704 }
705 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
706 Return (API->error);
707 }
708
709 GMT_Report (API, GMT_MSG_INFORMATION, "Working on seamount # %5d\n", smt);
710
711 gmt_M_free (GMT, Out);
712 gmt_M_free (GMT, hotspot);
713 gmt_M_free (GMT, orig_hotspot);
714 gmt_M_free (GMT, hot);
715 gmt_M_free (GMT, p);
716
717 Return (GMT_NOERROR);
718 }
719
720
GMT_originator(void * V_API,int mode,void * args)721 EXTERN_MSC int GMT_originator (void *V_API, int mode, void *args) {
722 /* This was the GMT5 name */
723 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
724 if (gmt_M_compat_check (API->GMT, 4)) {
725 GMT_Report (API, GMT_MSG_COMPAT, "Module originator is deprecated; use originater.\n");
726 return (GMT_Call_Module (API, "originater", mode, args));
727 }
728 GMT_Report (API, GMT_MSG_ERROR, "Shared GMT module not found: originator\n");
729 return (GMT_NOT_A_VALID_MODULE);
730 }
731
732