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 * Program for converting between total reconstruction and stage poles or to add rotations.
18 *
19 * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
20 * Date: 24-OCT-2001
21 * Version: 1.0
22 * 31-MAR-2006: Changed -H to -C to avoid clash with GMT -H
23 *
24 *-------------------------------------------------------------------------
25 * An ASCII stage pole (Euler) file must have following format:
26 *
27 * 1. Any number of comment lines starting with # in first column
28 * 2. Any number of blank lines (just carriage return, no spaces)
29 * 2. Any number of stage pole records which each have the format:
30 * lon(deg) lat(deg) tstart(Ma) tstop(Ma) ccw-angle(deg)
31 * 3. stage records must go from oldest to youngest rotation
32 * 4. Note tstart is larger (older) that tstop for each record
33 * 5. No gaps allowed: tstart must equal the previous records tstop
34 *
35 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
36 *
37 * # Time in Ma, angles in degrees
38 * # lon lat tstart tend ccw-angle
39 * 165 85 150 100 24.0
40 * 284 36 100 74 15.0
41 * 265 22 74 65 7.5
42 * 253 17 65 42 14.0
43 * 285 68 42 0 34.0
44 *
45 * AN ASCII total reconstruction file must have the following format:
46 *
47 * 1. Any number of comment lines starting with # in first column
48 * 2. Any number of blank lines (just carriage return, no spaces)
49 * 2. Any number of finite pole records which each have the format:
50 * lon(deg) lat(deg) tstop(Ma) ccw-angle(deg)
51 * 3. total reconstruction rotations must go from youngest to oldest
52 *
53 * Example: Duncan & Clague [1985] Pacific-Hotspot rotations:
54 *
55 * # Time in Ma, angles in degrees
56 * #longitude latitude time(My) angle(deg)
57 * 285.00000 68.00000 42.0000 34.0000
58 * 275.66205 53.05082 65.0000 43.5361
59 * 276.02501 48.34232 74.0000 50.0405
60 * 279.86436 46.30610 100.0000 64.7066
61 * 265.37800 55.69932 150.0000 82.9957
62 */
63
64 #include "gmt_dev.h"
65 #include "spotter.h"
66
67 #define THIS_MODULE_CLASSIC_NAME "rotconverter"
68 #define THIS_MODULE_MODERN_NAME "rotconverter"
69 #define THIS_MODULE_LIB "spotter"
70 #define THIS_MODULE_PURPOSE "Manipulate total reconstruction and stage rotations"
71 #define THIS_MODULE_KEYS ">D}"
72 #define THIS_MODULE_NEEDS ""
73 #define THIS_MODULE_OPTIONS "-:>Vh"
74
75 struct ROTCONVERTER_CTRL { /* All control options for this program (except common args) */
76 /* active is true if the option has been activated */
77 struct ROTCONVERTER_A { /* -A */
78 bool active;
79 } A;
80 struct ROTCONVERTER_D { /* -D */
81 bool active;
82 } D;
83 struct ROTCONVERTER_F { /* -F */
84 bool active;
85 bool mode; /* out mode (true if total reconstruction rotations) */
86 } F;
87 struct ROTCONVERTER_G { /* -G */
88 bool active;
89 } G;
90 struct ROTCONVERTER_M { /* -M[<value>] */
91 bool active;
92 double value;
93 } M;
94 struct ROTCONVERTER_N { /* -N */
95 bool active;
96 } N;
97 struct ROTCONVERTER_S { /* -S */
98 bool active;
99 } S;
100 struct ROTCONVERTER_T { /* -T */
101 bool active;
102 } T;
103 struct ROTCONVERTER_W { /* -W */
104 bool active;
105 } W;
106 };
107
New_Ctrl(struct GMT_CTRL * GMT)108 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
109 struct ROTCONVERTER_CTRL *C;
110
111 C = gmt_M_memory (GMT, NULL, 1, struct ROTCONVERTER_CTRL);
112
113 /* Initialize values whose defaults are not 0/false/NULL */
114
115 C->M.value = 0.5; /* To get half-angles */
116 C->F.mode = true; /* Default format is total reconstruction rotation poles for both input and output */
117 return (C);
118 }
119
Free_Ctrl(struct GMT_CTRL * GMT,struct ROTCONVERTER_CTRL * C)120 static void Free_Ctrl (struct GMT_CTRL *GMT, struct ROTCONVERTER_CTRL *C) { /* Deallocate control structure */
121 if (!C) return;
122 gmt_M_free (GMT, C);
123 }
124
usage(struct GMTAPI_CTRL * API,int level)125 static int usage (struct GMTAPI_CTRL *API, int level) {
126 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
127 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
128 GMT_Usage (API, 0, "usage: %s [+][-] <rotA> [[+][-] <rotB>] [[+][-] <rotC>] ... [-A] [-D] "
129 "[-Fs|t] [-G] [-M[<factor>]] [-N] [-S] [-T] [%s] [-W] [%s] [%s]\n", name, GMT_V_OPT, GMT_h_OPT, GMT_PAR_OPT);
130
131 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
132
133 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
134 GMT_Usage (API, 1, "\n<rotA>, ... are total reconstruction or stage rotation pole files. "
135 "Alternatively, give two plate IDs separated by a hyphen (e.g., PAC-MBL) "
136 "to extract that rotation from the GPlates rotation database (if installed). "
137 "Or, they can be a single rotation in lon/lat[/tstart[/tstop]]/angle format. "
138 "All rotation poles are assumed to be in geocentric coordinates. "
139 "Rotations will be added/subtracted in the order given.");
140 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
141 GMT_Usage (API, 1, "\n-A Report angles as time [Default uses time].");
142 GMT_Usage (API, 1, "\n-D Report all longitudes in -180/+180 range [Default is 0-360].");
143 GMT_Usage (API, 1, "\n-Fs|t");
144 GMT_Usage (API, -2, "Set output file directive:");
145 GMT_Usage (API, 3, "t: Total reconstruction [Default].");
146 GMT_Usage (API, 3, "s: Stage rotations.");
147 GMT_Usage (API, 1, "\n-G Write rotations using GPlates format [Default is spotter format].");
148 GMT_Usage (API, 1, "\n-M[<factor>]");
149 GMT_Usage (API, -2, "Reduce opening angles for stage rotations by <factor> [0.5]. "
150 "Typically used to get half-rates needed for flowlines.");
151 GMT_Usage (API, 1, "\n-N Ensure all poles are in northern hemisphere [Default ensures positive opening angles/rates].");
152 GMT_Usage (API, 1, "\n-S Ensure all poles are in southern hemisphere [Default ensures positive opening angles/rates].");
153 GMT_Usage (API, 1, "\n-T Transpose the result (i.e., change sign of final rotation angle).");
154 GMT_Option (API, "V");
155 GMT_Usage (API, 1, "\n-W Ensure all rotations have negative opening angles/rates [Default ensures positive opening angles/rates].");
156 GMT_Usage (API, 1, "\nOnly one of -N, -S, -W may be used at the same time.");
157 GMT_Option (API, "h,.");
158
159 return (GMT_MODULE_USAGE);
160 }
161
parse(struct GMT_CTRL * GMT,struct ROTCONVERTER_CTRL * Ctrl,struct GMT_OPTION * options)162 static int parse (struct GMT_CTRL *GMT, struct ROTCONVERTER_CTRL *Ctrl, struct GMT_OPTION *options) {
163 /* This parses the options provided to rotconverter and sets parameters in CTRL.
164 * Any GMT common options will override values set previously by other commands.
165 * It also replaces any file names specified as input or output with the data ID
166 * returned when registering these sources/destinations with the API.
167 */
168
169 unsigned int n_errors = 0;
170 struct GMT_OPTION *opt = NULL;
171 struct GMTAPI_CTRL *API = GMT->parent;
172
173 for (opt = options; opt; opt = opt->next) {
174 switch (opt->option) {
175
176 case '<': /* Skip input files */
177 break;
178
179 /* Supplemental parameters */
180
181 case 'A': /* Angle, not time */
182 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
183 Ctrl->A.active = true;
184 break;
185 case 'D': /* Dateline */
186 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
187 Ctrl->D.active = true;
188 break;
189
190 case 'E': /* Convert to half-spreading stage rotations [NOW -M] */
191 if (gmt_M_compat_check (GMT, 5)) { /* Warn and fall through */
192 GMT_Report (API, GMT_MSG_COMPAT, "-E is deprecated; use -M instead.\n");
193 Ctrl->M.active = true;
194 if (opt->arg[0]) Ctrl->M.value = atof (opt->arg);
195 }
196 else {
197 GMT_Report (API, GMT_MSG_ERROR, "No such option -E\n");
198 n_errors++;
199 }
200 break;
201
202 case 'F':
203 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
204 Ctrl->F.active = true;
205 if (strlen (opt->arg) != 1) {
206 GMT_Report (API, GMT_MSG_ERROR, "Must specify -F<out>\n");
207 n_errors++;
208 continue;
209 }
210 switch (opt->arg[0]) { /* Output format */
211 case 'f':
212 if (gmt_M_compat_check (GMT, 4)) /* Warn and fall through */
213 GMT_Report (API, GMT_MSG_COMPAT, "-Ff is deprecated; use -Ft instead.\n");
214 /* Intentionally fall through - to 't' */
215 else {
216 GMT_Report (API, GMT_MSG_ERROR, "Must specify t|s\n");
217 n_errors++;
218 break;
219 }
220 /* Intentionally fall through */
221 case 't':
222 Ctrl->F.mode = true;
223 break;
224 case 's':
225 Ctrl->F.mode = false;
226 break;
227 default:
228 GMT_Report (API, GMT_MSG_ERROR, "Must specify t|s\n");
229 n_errors++;
230 break;
231 }
232 break;
233
234 case 'G': /* GPlates output format */
235 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
236 Ctrl->G.active = true;
237 break;
238
239 case 'M': /* Convert to total reconstruction rotation poles instead */
240 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
241 Ctrl->M.active = true;
242 if (opt->arg[0]) Ctrl->M.value = atof (opt->arg);
243 break;
244
245 case 'N': /* Ensure all poles reported are in northern hemisphere */
246 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
247 Ctrl->N.active = true;
248 break;
249
250 case 'S': /* Ensure all poles reported are in southern hemisphere */
251 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
252 Ctrl->S.active = true;
253 break;
254
255 case 'T': /* Transpose the final result (i.e., change sign of rotation) */
256 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
257 Ctrl->T.active = true;
258 break;
259
260 case 'W': /* Ensure all poles reported have negative opening angles */
261 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
262 Ctrl->W.active = true;
263 break;
264
265 case '0': case '1': case '2': case '3': case '4': case '5': case '6':
266 case '7': case '8': case '9': case '.':
267 break; /* Probably a rotation lon/lat/angle with negative longitude */
268
269 default: /* Report bad options */
270 n_errors += gmt_default_error (GMT, opt->option);
271 break;
272 }
273 }
274
275 n_errors += gmt_M_check_condition (GMT, (Ctrl->S.active + Ctrl->N.active + Ctrl->W.active) > 1, "Specify only one of -N, -S, and -W!\n");
276 n_errors += gmt_M_check_condition (GMT, Ctrl->M.active && Ctrl->F.mode, "Option -M requires stage rotations on output. Please add -Fs\n");
277 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && !Ctrl->F.mode, "Option -G requires total reconstruction rotations on output\n");
278
279 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
280 }
281
282 #define bailout(code) {gmt_M_free_options (mode); return (code);}
283 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
284
GMT_rotconverter(void * V_API,int mode,void * args)285 EXTERN_MSC int GMT_rotconverter (void *V_API, int mode, void *args) {
286 struct EULER *p = NULL; /* Pointer to array of stage poles */
287 struct EULER *a = NULL, *b = NULL; /* Pointer to arrays of stage poles */
288
289 unsigned int col, k, stage; /* Misc. counters */
290 unsigned int n_slash, n_out = 0, n_opt = 0, n_p, n_a = 1, n_b;
291 int last_sign;
292 bool confusion = false, online_stage = false;
293 int error = 0; /* nonzero if arguments are inconsistent */
294 bool first = true; /* true for first input file */
295 bool online_rot = false; /* true if we gave a rotation on the commandline rather than file name */
296 bool no_time = false; /* true if we gave a rotation on the commandline as lon/lat/angle only */
297
298 double zero = 0.0; /* Needed to pass to spotter_init */
299 double lon = 0.0, lat = 0.0; /* Pole location for online rotations */
300 double t0 = 0.0, t1 = 0.0; /* Start, stop times for online rotations */
301 double angle = 0.0; /* Rotation angle for online rotations */
302 double out[20];
303
304 char *start_text[2] = {"tstart(My)", "astart(deg)"}; /* Misc. column titles for rates or angles */
305 char *end_text[2] = {"tend(My)", "aend(deg)"};
306 char *time_text[2] = {"ttime(My)", "tangle(deg)"};
307 char record[GMT_BUFSIZ] = {""};
308
309 struct GMT_RECORD *Out = NULL;
310 struct GMT_OPTION *ptr = NULL, *opt = NULL;
311 struct ROTCONVERTER_CTRL *Ctrl = NULL;
312 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
313 struct GMT_OPTION *options = NULL;
314 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
315
316 /*----------------------- Standard module initialization and parsing ----------------------*/
317
318 if (API == NULL) return (GMT_NOT_A_SESSION);
319 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
320 options = GMT_Create_Options (API, mode, args);
321 if (API->error) bailout (API->error); /* Set or get option list */
322
323 /* Special preprocessing since online rotations like -144/34/-9 and -.55/33/2 will
324 * have been decoded as options -4 and option -., respectively. Here we simply
325 * undo these and make them all "file" options -<.
326 * Also, a single - sign would have been decoded as the synopsis option. */
327
328 for (opt = options; opt; opt = opt->next, n_opt++) {
329 switch (opt->option) {
330 case '0': case '1': case '2': case '3': case '4': case '5':
331 case '6': case '7': case '8': case '9': case '.':
332 sprintf (record, "-%c%s", opt->option, opt->arg);
333 gmt_M_str_free (opt->arg);
334 opt->arg = strdup (record);
335 opt->option = GMT_OPT_INFILE;
336 break;
337 case GMT_OPT_SYNOPSIS:
338 gmt_M_str_free (opt->arg);
339 opt->arg = strdup ("-");
340 opt->option = GMT_OPT_INFILE;
341 confusion = true; /* Since we don't know if just a single - was given */
342 break;
343 default: /* Do nothing */
344 break;
345 }
346 }
347 if (!options || options->option == GMT_OPT_USAGE) bailout (usage (API, GMT_USAGE)); /* Return the usage message */
348 if (n_opt == 1 && confusion) bailout (usage (API, GMT_SYNOPSIS)); /* Return the synopsis */
349
350 /* Parse the command-line arguments */
351
352 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 */
353 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
354 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
355 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
356 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
357
358 /*---------------------------- This is the rotconverter main code ----------------------------*/
359
360 gmt_M_memset (out, 20, double);
361 if (Ctrl->G.active) {
362 gmt_set_column_type (GMT, GMT_OUT, 0, GMT_IS_FLOAT);
363 gmt_set_column_type (GMT, GMT_OUT, 1, GMT_IS_FLOAT);
364 gmt_set_column_type (GMT, GMT_OUT, 2, GMT_IS_LAT);
365 gmt_set_column_type (GMT, GMT_OUT, 3, GMT_IS_LON);
366 strcpy (GMT->current.setting.format_float_out, "%g");
367 }
368
369 last_sign = +1;
370 for (opt = options; opt; opt = opt->next) {
371 if (opt->option != GMT_OPT_INFILE) continue; /* Only consider files (or rotations) from here on */
372
373 /* Here a single + would have been parsed as a file with name "+"; change to a plus sign */
374 if (opt->option == GMT_OPT_INFILE && opt->arg[0] == '+' && opt->arg[1] == '\0') {
375 last_sign = +1;
376 continue;
377 }
378 /* Here a single - would have been parsed as a file with name "-"; change to a minus sign */
379 if (opt->option == GMT_OPT_INFILE && opt->arg[0] == '-' && opt->arg[1] == '\0') {
380 last_sign = -1;
381 continue;
382 }
383
384 if (spotter_GPlates_pair (opt->arg)) { /* A GPlates plate pair to look up in the rotation table */
385 online_rot = false;
386 }
387 else if (gmt_access (GMT, opt->arg, R_OK)) { /* Not a readable file, is it a lon/lat/t0[/t1]/omega specification? */
388 for (k = n_slash = 0; opt->arg[k]; k++) if (opt->arg[k] == '/') n_slash++;
389 if (n_slash < 2 || n_slash > 4) { /* No way it can be a online rotation, cry foul */
390 GMT_Report (API, GMT_MSG_ERROR, "Cannot read file %s\n", opt->arg);
391 gmt_M_free (GMT, a);
392 Return (GMT_RUNTIME_ERROR);
393 }
394 else { /* Try to decode as a single rotation */
395
396 k = sscanf (opt->arg, "%lf/%lf/%lf/%lf/%lf", &lon, &lat, &t0, &t1, &angle);
397 if (k == 4) angle = t1, t1 = 0.0; /* Only 4 input values */
398 if (n_slash == 2) angle = t0, t0 = 1.0, t1 = 0.0, no_time = true; /* Quick lon/lat/angle total reconstruction rotation, no time */
399 if (t0 < t1) {
400 GMT_Report (API, GMT_MSG_ERROR, "Online rotation has t_start (%g) younger than t_stop (%g)\n", t0, t1);
401 Return (GMT_RUNTIME_ERROR);
402 }
403 if (angle == 0.0) {
404 GMT_Report (API, GMT_MSG_ERROR, "Online rotation has zero opening angle\n");
405 Return (GMT_RUNTIME_ERROR);
406 }
407 online_rot = true;
408 online_stage = (k == 5);
409 }
410 }
411 else
412 online_rot = false;
413
414 if (first) { /* First time loading a rotation model */
415 if (online_rot) {
416 n_a = 1;
417 a = gmt_M_memory (GMT, NULL, 1, struct EULER);
418 a[0].lon = lon; a[0].lat = lat;
419 a[0].t_start = t0; a[0].t_stop = t1;
420 a[0].duration = t0 - t1;
421 a[0].omega = angle / a[0].duration;
422 if (online_stage) spotter_stages_to_total (GMT, a, n_a, true, true);
423 }
424 else {
425 if ((error = spotter_init (GMT, opt->arg, &a, 0, true, false, &zero)) < 0) /* Return total reconstruction rotations */
426 Return (-error);
427 n_a = (unsigned int)error;
428 }
429 zero = 0.0;
430 if (last_sign == -1) { /* Leading - sign, simply reverse the rotation angles */
431 for (stage = 0; stage < n_a; stage++) {
432 a[stage].omega = -a[stage].omega;
433 spotter_cov_of_inverse (GMT, &a[stage], a[stage].C);
434 }
435 last_sign = +1;
436 }
437 first = false;
438 }
439 else { /* For additional times, load a second model and add/subtract them */
440 if (online_rot) {
441 n_b = 1;
442 b = gmt_M_memory (GMT, NULL, 1, struct EULER);
443 b[0].lon = lon; b[0].lat = lat;
444 b[0].t_start = t0; b[0].t_stop = t1;
445 b[0].duration = t0 - t1;
446 b[0].omega = angle / b[0].duration;
447 if (online_stage) spotter_stages_to_total (GMT, b, n_b, true, true);
448 }
449 else {
450 if ((error = spotter_init (GMT, opt->arg, &b, 0, true, false, &zero)) < 0) /* Return total reconstruction rotations */
451 Return (-error);
452 n_b = (unsigned int)error;
453 }
454 zero = 0.0;
455 spotter_add_rotations (GMT, a, n_a, b, last_sign * n_b, &p, &n_p); /* Add the two total reconstruction rotations sets, returns total reconstruction rotations in p */
456 gmt_M_free (GMT, a);
457 gmt_M_free (GMT, b);
458 a = p;
459 n_a = n_p;
460 }
461 }
462
463 if (a == NULL) {
464 GMT_Report (API, GMT_MSG_ERROR, "No rotation resulting from operation\n");
465 Return (GMT_RUNTIME_ERROR);
466 }
467 n_out = 3 + ((Ctrl->F.mode) ? 1 - no_time : 2);
468 if (online_rot && n_out > 3) n_out--;
469 if (a[0].has_cov) n_out += 9;
470 if (Ctrl->G.active) n_out = 6;
471
472 if ((error = GMT_Set_Columns (API, GMT_OUT, n_out, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
473 gmt_M_free (GMT, a);
474 Return (error);
475 }
476 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
477 gmt_M_free (GMT, a);
478 Return (API->error);
479 }
480 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
481 gmt_M_free (GMT, a);
482 Return (API->error);
483 }
484 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_NONE) != GMT_NOERROR) { /* Sets output geometry */
485 gmt_M_free (GMT, a);
486 Return (API->error);
487 }
488
489 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
490 if (Ctrl->G.active) /* GPlates header */
491 sprintf (record, "plateid%stime%slatitude%slongitude%sangle%sfixedplateid\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, \
492 GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator);
493 else if (Ctrl->F.mode && no_time)
494 sprintf (record, "longitude%slatitude%sangle(deg)\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator);
495 else if (Ctrl->F.mode) /* Easy, simply output what we've got following a header*/
496 sprintf (record, "longitude%slatitude%s%s%sangle(deg)\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, time_text[Ctrl->A.active], GMT->current.setting.io_col_separator);
497 else { /* Convert total reconstruction to stages before output */
498 spotter_total_to_stages (GMT, a, n_a, true, true); /* To ensure we have the right kind of poles for output */
499 printf (record, "longitude%slatitude%s%s%s%s%sangle(deg)\n", GMT->current.setting.io_col_separator, GMT->current.setting.io_col_separator, start_text[Ctrl->A.active], GMT->current.setting.io_col_separator, end_text[Ctrl->A.active], GMT->current.setting.io_col_separator);
500 }
501 if (GMT->current.setting.io_header[GMT_OUT]) GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, record);
502
503 for (stage = 0; stage < n_a; stage++) {
504 if (Ctrl->T.active) a[stage].omega = -a[stage].omega;
505 if ((Ctrl->W.active && a[stage].omega > 0.0) || (Ctrl->S.active && a[stage].lat > 0.0) || (Ctrl->N.active && a[stage].lat < 0.0) || (!(Ctrl->S.active || Ctrl->N.active || Ctrl->W.active) && a[stage].omega < 0.0)) /* flip to antipole */
506 a[stage].lat = -a[stage].lat, a[stage].lon += 180.0, a[stage].omega = -a[stage].omega;
507 while (a[stage].lon >= 360.0) a[stage].lon -= 360.0; /* Force geodetic longitude range */
508 while (a[stage].lon <= -180.0) a[stage].lon += 360.0; /* Force geodetic longitude range */
509 while (!Ctrl->D.active && a[stage].lon < 0.0) a[stage].lon += 360.0; /* Force geodetic longitude range */
510 while (Ctrl->D.active && a[stage].lon > 180.0) a[stage].lon -= 360.0; /* Force a geographic longitude range */
511 if (Ctrl->M.active) a[stage].omega *= Ctrl->M.value;
512 out[GMT_X] = a[stage].lon; out[GMT_Y] = a[stage].lat;
513 col = 2;
514 if (Ctrl->G.active) {
515 out[0] = (double)a[stage].id[0]; out[5] = (double)a[stage].id[1];
516 out[1] = a[stage].t_start;
517 out[2] = a[stage].lat;
518 out[3] = a[stage].lon;
519 out[4] = a[stage].omega * a[stage].duration;
520 }
521 else if (Ctrl->F.mode && no_time) {
522 out[col++] = a[stage].omega * a[stage].duration;
523 }
524 else if (Ctrl->F.mode) {
525 out[col++] = a[stage].t_start;
526 out[col++] = a[stage].omega * a[stage].duration;
527 }
528 else {
529 out[col++] = a[stage].t_start;
530 out[col++] = a[stage].t_stop;
531 out[col++] = a[stage].omega * a[stage].duration;
532 }
533 if (a[stage].has_cov) {
534 double K[9];
535 spotter_covar_to_record (GMT, &a[stage], K);
536 for (k = 0; k < 9; k++) out[col++] = K[k];
537 }
538 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
539 }
540
541 gmt_M_free (GMT, a);
542 gmt_M_free (GMT, Out);
543
544 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
545 Return (API->error);
546 }
547
548 Return (GMT_NOERROR);
549 }
550