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 * gmtpmodeler will read an lon ,at[, age] table and a plate motion model and
18 * evaluates chosen attributes such as speed, direction, lon/lat origin
19 * for each input point.
20 *
21 * Author: Paul Wessel
22 * Date: 30-NOV-2015
23 * Ver: 5.3
24 */
25
26 #include "gmt_dev.h"
27 #include "spotter.h"
28
29 #define THIS_MODULE_CLASSIC_NAME "gmtpmodeler"
30 #define THIS_MODULE_MODERN_NAME "gmtpmodeler"
31 #define THIS_MODULE_LIB "spotter"
32 #define THIS_MODULE_PURPOSE "Evaluate a plate motion model at given locations"
33 #define THIS_MODULE_KEYS "<D{,FD(,>D}"
34 #define THIS_MODULE_NEEDS ""
35 #define THIS_MODULE_OPTIONS "-:>Vbdefghioqs"
36
37 #define N_PM_ITEMS 9
38 #define PM_AZIM 0
39 #define PM_DIST 1
40 #define PM_STAGE 2
41 #define PM_VEL 3
42 #define PM_OMEGA 4
43 #define PM_DLON 5
44 #define PM_DLAT 6
45 #define PM_LON 7
46 #define PM_LAT 8
47
48 struct GMTPMODELER_CTRL { /* All control options for this program (except common args) */
49 /* active is true if the option has been activated */
50 struct GMTPMODELER_In {
51 bool active;
52 char *file;
53 } In;
54 struct GMTPMODELER_E { /* -E[+]rotfile, -E[+]<ID1>-<ID2>, or -E<lon/lat/angle> */
55 bool active;
56 struct SPOTTER_ROT rot;
57 } E;
58 struct GMTPMODELER_F { /* -Fpolfile */
59 bool active;
60 char *file;
61 } F;
62 struct GMTPMODELER_N { /* -N */
63 bool active;
64 double t_upper;
65 } N;
66 struct GMTPMODELER_S { /* -Sa|d|s|v|w|x|y|X|Y */
67 bool active, center;
68 unsigned int mode[N_PM_ITEMS];
69 unsigned int n_items;
70 } S;
71 struct GMTPMODELER_T { /* -T<fixtime> */
72 bool active;
73 double value;
74 } T;
75 };
76
New_Ctrl(struct GMT_CTRL * GMT)77 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
78 struct GMTPMODELER_CTRL *C;
79
80 C = gmt_M_memory (GMT, NULL, 1, struct GMTPMODELER_CTRL);
81
82 return (C);
83 }
84
Free_Ctrl(struct GMT_CTRL * GMT,struct GMTPMODELER_CTRL * C)85 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GMTPMODELER_CTRL *C) { /* Deallocate control structure */
86 if (!C) return;
87 gmt_M_str_free (C->In.file);
88 gmt_M_str_free (C->E.rot.file);
89 gmt_M_str_free (C->F.file);
90 gmt_M_free (GMT, C);
91 }
92
usage(struct GMTAPI_CTRL * API,int level)93 static int usage (struct GMTAPI_CTRL *API, int level) {
94 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
95 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
96 GMT_Usage (API, 0, "usage: %s <table> %s [-F<polygontable>] [-N<upper_age>] [-SadrswxyXY] [-T<time>] "
97 "[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
98 name, SPOTTER_E_OPT, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_q_OPT, GMT_PAR_OPT);
99
100 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
101
102 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
103 GMT_Usage (API, 1, "\n<table> is a data table with geographic coordinates and optionally crustal ages.");
104 spotter_rot_usage (API);
105 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
106 GMT_Usage (API, 1, "\n-F<polygontable>");
107 GMT_Usage (API, -2, "Specify a multi-segment closed polygon file that describes the area "
108 "of the data table to work on [Default works on the entire table].");
109 GMT_Usage (API, 1, "\n-N<upper_age>");
110 GMT_Usage (API, -2, "Extend earliest stage pole back to <upper_age> [no extension].");
111 GMT_Option (API, "Rg");
112 GMT_Usage (API, 1, "\n-SadrswxyXY");
113 GMT_Usage (API, -2, "Select one or more model predictions as a function of crustal age. Choose from:");
114 GMT_Usage (API, 3, "a: Plate spreading azimuth.");
115 GMT_Usage (API, 3, "d: Distance to origin of crust in km.");
116 GMT_Usage (API, 3, "r: Plate motion rate in mm/yr or km/Myr.");
117 GMT_Usage (API, 3, "s: Plate motion stage ID (1 is youngest).");
118 GMT_Usage (API, 3, "w: Rotation rate in degrees/Myr.");
119 GMT_Usage (API, 3, "x: Change in longitude since formation.");
120 GMT_Usage (API, 3, "y: Change in latitude since formation.");
121 GMT_Usage (API, 3, "X: Longitude at origin of crust.");
122 GMT_Usage (API, 3, "Y: Latitude at origin of crust.");
123 GMT_Usage (API, -2, "Default writes lon,lat,age,<adrswxyXY> to standard output.");
124 GMT_Usage (API, 1, "\n-T<time>");
125 GMT_Usage (API, -2, "Set fixed time of reconstruction to override any input ages.");
126 GMT_Option (API, "bi3,bo,d,e,f,h,i,o,q,s,:,.");
127
128 return (GMT_MODULE_USAGE);
129 }
130
parse(struct GMT_CTRL * GMT,struct GMTPMODELER_CTRL * Ctrl,struct GMT_OPTION * options)131 static int parse (struct GMT_CTRL *GMT, struct GMTPMODELER_CTRL *Ctrl, struct GMT_OPTION *options) {
132 /* This parses the options provided to gmtpmodeler and sets parameters in CTRL.
133 * Any GMT common options will override values set previously by other commands.
134 * It also replaces any file names specified as input or output with the data ID
135 * returned when registering these sources/destinations with the API.
136 */
137
138 unsigned int n_errors = 0, n_files = 0, k;
139 struct GMT_OPTION *opt = NULL;
140 struct GMTAPI_CTRL *API = GMT->parent;
141
142 for (opt = options; opt; opt = opt->next) {
143 switch (opt->option) {
144
145 case '<': /* Input files */
146 if (n_files++ > 0) break;
147 Ctrl->In.active = true;
148 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
149 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
150 break;
151
152 /* Supplemental parameters */
153
154 case 'E': /* File with stage poles */
155 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
156 Ctrl->E.active = true;
157 n_errors += spotter_parse (GMT, opt->option, opt->arg, &(Ctrl->E.rot));
158 break;
159 case 'F':
160 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
161 Ctrl->F.active = true;
162 if (opt->arg[0]) Ctrl->F.file = strdup (opt->arg);
163 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
164 break;
165 case 'N': /* Extend oldest stage back to this time [no extension] */
166 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
167 Ctrl->N.active = true;
168 Ctrl->N.t_upper = atof (opt->arg);
169 break;
170 case 'S':
171 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
172 Ctrl->S.active = true;
173 while (opt->arg[Ctrl->S.n_items]) {
174 switch (opt->arg[Ctrl->S.n_items]) {
175 case 'a': /* Plate spreading azimuth */
176 Ctrl->S.mode[Ctrl->S.n_items] = PM_AZIM; break;
177 case 'd': /* Distance from point to origin at ridge */
178 Ctrl->S.mode[Ctrl->S.n_items] = PM_DIST; break;
179 case 's': /* Plate motion stage ID */
180 Ctrl->S.mode[Ctrl->S.n_items] = PM_STAGE; break;
181 case 'v': case 'r': /* Plate spreading rate [r is backwards compatible] */
182 Ctrl->S.mode[Ctrl->S.n_items] = PM_VEL; break;
183 case 'w': /* Plate rotation omega */
184 Ctrl->S.mode[Ctrl->S.n_items] = PM_OMEGA; break;
185 case 'x': /* Change in longitude since origin */
186 Ctrl->S.mode[Ctrl->S.n_items] = PM_DLON; Ctrl->S.center = true; break;
187 case 'y': /* Change in latitude since origin */
188 Ctrl->S.mode[Ctrl->S.n_items] = PM_DLAT; break;
189 case 'X': /* Plate longitude at crust origin */
190 Ctrl->S.mode[Ctrl->S.n_items] = PM_LON; break;
191 case 'Y': /* Plate latitude at crust origin */
192 Ctrl->S.mode[Ctrl->S.n_items] = PM_LAT; break;
193 default:
194 n_errors++; break;
195 }
196 Ctrl->S.n_items++;
197 }
198 break;
199 case 'T':
200 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
201 Ctrl->T.active = true;
202 Ctrl->T.value = atof (opt->arg);
203 break;
204
205 default: /* Report bad options */
206 n_errors += gmt_default_error (GMT, opt->option);
207 break;
208 }
209 }
210
211 if (Ctrl->S.n_items == 0) { /* Set default which are all the items */
212 Ctrl->S.active = true;
213 Ctrl->S.n_items = N_PM_ITEMS;
214 for (k = 0; k < Ctrl->S.n_items; k++) Ctrl->S.mode[k] = k;
215 }
216
217 n_errors += gmt_M_check_condition (GMT, !Ctrl->E.active, "Must specify -E\n");
218 n_errors += gmt_M_check_condition (GMT, !Ctrl->S.active, "Must specify -S\n");
219 n_errors += gmt_M_check_condition (GMT, Ctrl->S.n_items == 0, "Must specify one or more fields with -S\n");
220 n_errors += gmt_M_check_condition (GMT, Ctrl->T.value < 0.0, "Option -T: Must specify positive age.\n");
221
222 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
223 }
224
225 #define bailout(code) {gmt_M_free_options (mode); return (code);}
226 #define Return(code) {if (p) gmt_M_free (GMT, p); Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
227
GMT_gmtpmodeler(void * V_API,int mode,void * args)228 EXTERN_MSC int GMT_gmtpmodeler (void *V_API, int mode, void *args) {
229 unsigned int inside, stage, n_stages, k;
230 int retval, error = 0, n_fields;
231
232 bool spotted;
233
234 uint64_t seg, n_old = 0, n_outside = 0, n_NaN = 0, n_read = 0;
235
236 double lon = 0, lat = 0, lat_c, d, value = 0.0, age, *in = NULL, *out = NULL;
237
238 char *tag[N_PM_ITEMS] = { "az", "dist", "stage", "vel", "omega", "dlon", "dlat", "lon", "lat" };
239 struct GMT_DATASET *D = NULL;
240 struct GMT_DATATABLE *pol = NULL;
241 struct EULER *p = NULL; /* Pointer to array of stage poles */
242 struct GMT_RECORD *In = NULL, *Out = NULL;
243 struct GMT_OPTION *ptr = NULL;
244 struct GMTPMODELER_CTRL *Ctrl = NULL;
245 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
246 struct GMT_OPTION *options = NULL;
247 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
248
249 /*----------------------- Standard module initialization and parsing ----------------------*/
250
251 if (API == NULL) return (GMT_NOT_A_SESSION);
252 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
253 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
254
255 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
256
257 /* Parse the command-line arguments */
258
259 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 */
260 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
261 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
262 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
263 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
264
265 /*---------------------------- This is the gmtpmodeler main code ----------------------------*/
266
267 if (Ctrl->F.active) { /* Read the user's clip polygon file */
268 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from -C,-F,-L files */
269 if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
270 Return (API->error);
271 }
272 if (D->n_columns < 2) {
273 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least 2 are needed\n", (int)D->n_columns);
274 Return (GMT_DIM_TOO_SMALL);
275 }
276 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
277 pol = D->table[0]; /* Since it is a single file */
278 GMT_Report (API, GMT_MSG_INFORMATION, "Restrict evaluation to within polygons in file %s\n", Ctrl->F.file);
279 gmt_set_inside_mode (GMT, D, GMT_IOO_UNKNOWN);
280 }
281
282 if (Ctrl->E.rot.single) { /* Got a single rotation, no time, create a rotation table with one entry */
283 n_stages = 1;
284 p = gmt_M_memory (GMT, NULL, n_stages, struct EULER);
285 p[0].lon = Ctrl->E.rot.lon; p[0].lat = Ctrl->E.rot.lat; p[0].omega = Ctrl->E.rot.w;
286 if (gmt_M_is_dnan (Ctrl->E.rot.age)) { /* No age, use fake age = 1 everywhere */
287 Ctrl->T.active = true;
288 Ctrl->T.value = Ctrl->N.t_upper = p[0].t_start = 1.0;
289 }
290 spotter_setrot (GMT, &(p[0]));
291 }
292 else { /* Got a file or Gplates plate pair */
293 if ((retval = spotter_init (GMT, Ctrl->E.rot.file, &p, 0, false, Ctrl->E.rot.invert, &Ctrl->N.t_upper)) < 0)
294 Return (-retval);
295 n_stages = (unsigned int)retval;
296 }
297
298 for (stage = 0; stage < n_stages; stage++) {
299 if (p[stage].omega < 0.0) { /* Ensure all stages have positive rotation angles */
300 p[stage].omega = -p[stage].omega;
301 p[stage].lat = -p[stage].lat;
302 p[stage].lon += 180.0;
303 if (p[stage].lon > 360.0) p[stage].lon -= 360.0;
304 }
305 }
306 if (Ctrl->T.active && Ctrl->T.value > Ctrl->N.t_upper) {
307 GMT_Report (API, GMT_MSG_ERROR, "Requested a fixed reconstruction time outside range of rotation table\n");
308 Return (GMT_RUNTIME_ERROR);
309 }
310 /* Set up input */
311 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
312 Return (API->error);
313 }
314 if ((error = GMT_Set_Columns (API, GMT_IN, 2 + !Ctrl->T.active, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
315 Return (error);
316 }
317 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
318 Return (API->error);
319 }
320 /* Set up output */
321 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
322 Return (API->error);
323 }
324 if ((error = GMT_Set_Columns (API, GMT_OUT, Ctrl->S.n_items + 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
325 Return (error);
326 }
327 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
328 Return (API->error);
329 }
330 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
331 Return (API->error);
332 }
333
334 if (gmt_init_distaz (GMT, 'd', GMT_GREATCIRCLE, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE) /* Great circle distances in degrees */
335 Return (GMT_NOT_A_VALID_TYPE);
336 if (Ctrl->S.center) GMT->current.io.geo.range = GMT_IS_M180_TO_P180_RANGE; /* Need +- around 0 here */
337
338 out = gmt_M_memory (GMT, NULL, Ctrl->S.n_items + 3, double);
339 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
340 if (GMT->current.setting.io_header[GMT_OUT]) {
341 char header[GMT_BUFSIZ] = {""};
342 for (k = 0; k < Ctrl->S.n_items; k++) {
343 strncat (header, tag[k], GMT_BUFSIZ-1);
344 if (k < (Ctrl->S.n_items-1)) strncat (header, GMT->current.setting.io_col_separator, GMT_BUFSIZ-1);
345 }
346 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, header); /* Write a header record */
347 }
348
349 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate %d model predictions based on %s\n", Ctrl->S.n_items, Ctrl->E.rot.file);
350
351 /* Read the location data from file or stdin */
352
353 do { /* Keep returning records until we reach EOF */
354 n_read++;
355 if ((In = GMT_Get_Record (API, GMT_READ_DATA, &n_fields)) == NULL) { /* Read next record, get NULL if special case */
356 if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */
357 Return (GMT_RUNTIME_ERROR);
358 if (gmt_M_rec_is_table_header (GMT)) { /* Skip all table headers */
359 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, NULL);
360 continue;
361 }
362 if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
363 break;
364 else if (gmt_M_rec_is_new_segment (GMT)) { /* Parse segment headers */
365 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
366 continue;
367 }
368 }
369
370 if (In->data == NULL) {
371 gmt_quit_bad_record (API, In);
372 Return (API->error);
373 }
374
375 /* Data record to process */
376 in = In->data;
377 Out->text = In->text;
378
379 if (Ctrl->F.active) { /* Use the bounding polygon */
380 for (seg = inside = 0; seg < pol->n_segments && !inside; seg++) { /* Use degrees since function expects it */
381 if (gmt_polygon_is_hole (GMT, pol->segment[seg])) continue; /* Holes are handled within gmt_inonout */
382 inside = (gmt_inonout (GMT, in[GMT_X], in[GMT_Y], pol->segment[seg]) > GMT_OUTSIDE);
383 }
384 if (!inside) {
385 n_outside++; /* Outside the polygon(s), continue */
386 continue;
387 }
388 }
389 /* Get age for this node (or use the constant age) */
390 if (Ctrl->T.active)
391 age = Ctrl->T.value;
392 else if (n_fields == 3) {
393 age = in[GMT_Z];
394 if (gmt_M_is_dnan (age)) { /* No crustal age */
395 n_NaN++;
396 continue;
397 }
398 else if (age > Ctrl->N.t_upper) { /* Outside of model range */
399 n_old++;
400 continue;
401 }
402 }
403 else {
404 GMT_Report (API, GMT_MSG_WARNING, "Point %" PRIu64 " has no age (%g) (skipped)\n", n_read);
405 continue;
406 }
407
408 /* Here we are inside; get the coordinates and rotate back to original grid coordinates */
409 if ((retval = spotter_stage (GMT, age, p, n_stages)) < 0) continue; /* Outside valid stage rotation range */
410 stage = retval; /* Current rotation stage */
411 spotted = false; /* Not yet called spotter_backtrack at this node */
412 out[GMT_X] = in[GMT_X]; out[GMT_Y] = in[GMT_Y]; out[GMT_Z] = age; /* Set the first 3 output coordinates */
413 lat_c = gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O); /* Get concentric latitude */
414
415 for (k = 0; k < Ctrl->S.n_items; k++) { /* Loop over desired output components */
416 switch (Ctrl->S.mode[k]) {
417 case PM_AZIM: /* Compute plate motion direction at this point in time/space */
418 value = gmt_az_backaz (GMT, in[GMT_X], lat_c, p[stage].lon, p[stage].lat, false) + 90.0 * gmt_signum (p[stage].omega);
419 gmt_lon_range_adjust (GMT->current.io.geo.range, &value);
420 break;
421 case PM_DIST: /* Compute great-circle distance between node and point of origin at ridge */
422 if (!spotted) {
423 lon = in[GMT_X] * D2R; lat = lat_c * D2R;
424 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
425 Return (GMT_RUNTIME_ERROR);
426 spotted = true;
427 }
428 value = GMT->current.proj.DIST_KM_PR_DEG * gmt_distance (GMT, in[GMT_X], lat_c, lon * R2D, lat * R2D);
429 break;
430 case PM_STAGE: /* Compute plate rotation stage */
431 value = stage;
432 break;
433 case PM_VEL: /* Compute plate motion speed at this point in time/space */
434 d = gmt_distance (GMT, in[GMT_X], lat_c, p[stage].lon, p[stage].lat);
435 value = sind (d) * p[stage].omega * GMT->current.proj.DIST_KM_PR_DEG; /* km/Myr or mm/yr */
436 break;
437 case PM_OMEGA: /* Compute plate rotation rate omega */
438 value = p[stage].omega; /* degree/Myr */
439 break;
440 case PM_DLON: /* Compute latitude where this point was formed in the model */
441 if (!spotted) {
442 lon = in[GMT_X] * D2R; lat = lat_c * D2R;
443 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
444 Return (GMT_RUNTIME_ERROR);
445 spotted = true;
446 }
447 value = in[GMT_X] - lon * R2D;
448 if (fabs (value) > 180.0) value = copysign (360.0 - fabs (value), -value);
449 break;
450 case PM_DLAT: /* Compute latitude where this point was formed in the model */
451 if (!spotted) {
452 lon = in[GMT_X] * D2R; lat = lat_c * D2R;
453 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
454 Return (GMT_RUNTIME_ERROR);
455 spotted = true;
456 }
457 value = in[GMT_Y] - gmt_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */
458 break;
459 case PM_LON: /* Compute latitude where this point was formed in the model */
460 if (!spotted) {
461 lon = in[GMT_X] * D2R; lat = lat_c * D2R;
462 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
463 Return (GMT_RUNTIME_ERROR);
464 spotted = true;
465 }
466 value = lon * R2D;
467 break;
468 case PM_LAT: /* Compute latitude where this point was formed in the model */
469 if (!spotted) {
470 lon = in[GMT_X] * D2R; lat = lat_c * D2R;
471 if (spotter_backtrack (GMT, &lon, &lat, &age, 1U, p, n_stages, 0.0, 0.0, 0, NULL, NULL) < 0)
472 Return (GMT_RUNTIME_ERROR);
473 spotted = true;
474 }
475 value = gmt_lat_swap (GMT, lat * R2D, GMT_LATSWAP_O2G); /* Convert back to geodetic */
476 break;
477 }
478 out[k+3] = value;
479 }
480 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
481 } while (true);
482
483 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
484 Return (API->error);
485 }
486 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
487 Return (API->error);
488 }
489
490 if (n_outside) GMT_Report (API, GMT_MSG_WARNING, "%" PRIu64 " points fell outside the polygonal boundary\n", n_outside);
491 if (n_old) GMT_Report (API, GMT_MSG_WARNING, "%" PRIu64 " points had ages that exceeded the limit of the rotation model\n", n_old);
492 if (n_NaN) GMT_Report (API, GMT_MSG_WARNING, "%" PRIu64 " points had ages that were NaN\n", n_NaN);
493
494 gmt_M_free (GMT, out);
495 gmt_M_free (GMT, Out);
496
497 Return (GMT_NOERROR);
498 }
499