1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 2016-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 averaging finite rotations, resulting in mean rotations
18 * optionally with a covarience matrix describing the variation.
19 *
20 * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA
21 * Date: 6-AUG-2016
22 * Version: GMT 5 [derived from rotsmooth_age.c in rotsuite]
23 *
24 * AN ASCII total reconstruction file must have the following format:
25 *
26 * 1. Any number of comment lines starting with # in first column
27 * 2. Any number of blank lines (just carriage return, no spaces)
28 * 2. Any number of finite pole records which each have the format:
29 * lon(deg) lat(deg) age(Ma) ccw-angle(deg) [<weight>]
30 * 3. If no age is present the use -A to use angles as pseudo-ages
31 *
32 * Binary data files cannot have header records
33 */
34
35 #include "gmt_dev.h"
36 #include "spotter.h"
37
38 #define THIS_MODULE_CLASSIC_NAME "rotsmoother"
39 #define THIS_MODULE_MODERN_NAME "rotsmoother"
40 #define THIS_MODULE_LIB "spotter"
41 #define THIS_MODULE_PURPOSE "Get mean rotations and covariance matrices from set of finite rotations"
42 #define THIS_MODULE_KEYS "<D{,>D}"
43 #define THIS_MODULE_NEEDS ""
44 #define THIS_MODULE_OPTIONS "-:>Vbdefhios" GMT_OPT("HMm")
45
46 struct ROTSMOOTHER_CTRL { /* All control options for this program (except common args) */
47 /* active is true if the option has been activated */
48 struct ROTSMOOTHER_A { /* -A */
49 bool active;
50 } A;
51 struct ROTSMOOTHER_C { /* -C */
52 bool active;
53 } C;
54 struct ROTSMOOTHER_N { /* -N */
55 bool active;
56 } N;
57 struct ROTSMOOTHER_S { /* -S */
58 bool active;
59 } S;
60 struct ROTSMOOTHER_T { /* -T<time>, -T<start/stop/inc> or -T<tfile.txt> */
61 bool active;
62 unsigned int n_times; /* Number of reconstruction times */
63 double *value; /* Array with one or more reconstruction times */
64 } T;
65 struct ROTSMOOTHER_W { /* -W */
66 bool active;
67 } W;
68 struct ROTSMOOTHER_Z { /* -Z */
69 bool active;
70 } Z;
71 };
72
73 #define K_ANGLE 0
74 #define K_LON 1
75 #define K_LAT 2
76 #define K_AGE 3
77 #define K_WEIGHT 4
78 #define N_ITEMS 5
79
80 struct ROTSMOOTHER_AGEROT {
81 double wxyasn[N_ITEMS];
82 double P[3]; /* For Cartesian components */
83 double Q[4]; /* For quaternion components */
84 };
85
New_Ctrl(struct GMT_CTRL * GMT)86 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
87 struct ROTSMOOTHER_CTRL *C;
88
89 C = gmt_M_memory (GMT, NULL, 1, struct ROTSMOOTHER_CTRL);
90
91 /* Initialize values whose defaults are not 0/false/NULL */
92
93 return (C);
94 }
95
Free_Ctrl(struct GMT_CTRL * GMT,struct ROTSMOOTHER_CTRL * C)96 static void Free_Ctrl (struct GMT_CTRL *GMT, struct ROTSMOOTHER_CTRL *C) { /* Deallocate control structure */
97 if (!C) return;
98 gmt_M_free (GMT, C->T.value);
99 gmt_M_free (GMT, C);
100 }
101
usage(struct GMTAPI_CTRL * API,int level)102 static int usage (struct GMTAPI_CTRL *API, int level) {
103 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
104 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
105 GMT_Usage (API, 0, "usage: %s [<table>] [-A] [-C] [-N] [-S] [-T<time>] [%s] [-W] [-Z] "
106 "[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
107 name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT,
108 GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
109
110 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
111
112 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
113 GMT_Usage (API, 1, "\n<table> (in ASCII, binary, or netCDF) has 3 or more columns. If no file(s) is given, standard input is read. "
114 "First 4 columns must have lon, lat (or lat, lon, see -:), time, and angle (degrees).");
115 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
116 GMT_Usage (API, 1, "\n-A Use opening angles as time. Input is <lon> <lat> <angle> [<weight>] and -T refers to angles "
117 "[Default expects <lon> <lat> <time> <angle> [<weight>] and -T refers to time].");
118 GMT_Usage (API, 1, "\n-C Compute covariance matrix for each mean rotation.");
119 GMT_Usage (API, 1, "\n-N Ensure all poles are in northern hemisphere [Default ensures positive opening angles].");
120 GMT_Usage (API, 1, "\n-S Ensure all poles are in southern hemisphere [Default ensures positive opening angles].");
121 GMT_Usage (API, 1, "\n-T<time>");
122 GMT_Usage (API, -2, "Set the output times when a mean rotation and covariance matrix is desired. "
123 "Append a single time (-T<time>), an equidistant range of times (-T<min>/<max>/<inc>), "
124 "Append +n to t_inc to indicate the number of points instead of an increment. "
125 "Alternatively, give the name of a file with a list of times (-T<tfile>). "
126 "The times indicate bin-boundaries and we output the average rotation time per bin.");
127 GMT_Option (API, "V");
128 GMT_Usage (API, 1, "\n-W Expect weights in last column for a weighted mean rotation [no weights].");
129 GMT_Usage (API, 1, "\n-Z Report negative opening angles [positive].");
130 GMT_Option (API, "bi3,bo,d,e,f,h,i,o,s,:,.");
131
132 return (GMT_MODULE_USAGE);
133 }
134
parse(struct GMT_CTRL * GMT,struct ROTSMOOTHER_CTRL * Ctrl,struct GMT_OPTION * options)135 static int parse (struct GMT_CTRL *GMT, struct ROTSMOOTHER_CTRL *Ctrl, struct GMT_OPTION *options) {
136 /* This parses the options provided to backtracker and sets parameters in CTRL.
137 * Any GMT common options will override values set previously by other commands.
138 * It also replaces any file names specified as input or output with the data ID
139 * returned when registering these sources/destinations with the API.
140 */
141
142 unsigned int n_errors = 0, n_in = 3, t;
143 int k;
144 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, txt_c[GMT_LEN256] = {""};
145 struct GMT_OPTION *opt = NULL;
146 struct GMTAPI_CTRL *API = GMT->parent;
147
148 for (opt = options; opt; opt = opt->next) {
149 switch (opt->option) {
150
151 case '<': /* Input files */
152 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
153 break;
154
155 /* Supplemental parameters */
156
157 case 'A': /* Output only an age-limited segment of the track */
158 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
159 Ctrl->A.active = true;
160 break;
161
162 case 'C':
163 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
164 Ctrl->C.active = true;
165 break;
166
167 case 'N': /* Ensure all poles reported are in northern hemisphere */
168 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
169 Ctrl->N.active = true;
170 break;
171
172 case 'S': /* Ensure all poles reported are in southern hemisphere */
173 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
174 Ctrl->S.active = true;
175 break;
176
177 case 'T': /* New: -Tage, -Tmin/max/inc, -Tmin/max/n+, -Tfile */
178 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
179 Ctrl->T.active = true;
180 if (!gmt_access (GMT, opt->arg, R_OK)) { /* Gave a file with times in first column */
181 uint64_t seg, row;
182 struct GMT_DATASET *T = NULL;
183 struct GMT_DATASEGMENT *S = NULL;
184 if ((T = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_READ_NORMAL, NULL, opt->arg, NULL)) == NULL) {
185 GMT_Report (API, GMT_MSG_ERROR, "Failure while reading file %s\n", opt->arg);
186 n_errors++;
187 continue;
188 }
189 /* Single table, build t array */
190 Ctrl->T.n_times = (unsigned int)T->n_records;
191 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
192 for (seg = t = 0; seg < T->table[0]->n_segments; seg++) {
193 S = T->table[0]->segment[seg]; /* Shorthand to current segment */
194 for (row = 0; row < S->n_rows; seg++, t++)
195 Ctrl->T.value[t] = S->data[GMT_X][row];
196 }
197 if (GMT_Destroy_Data (API, &T) != GMT_NOERROR)
198 n_errors++;
199 break;
200 }
201 /* Not a file */
202 k = sscanf (opt->arg, "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c);
203 if (k == 3) { /* Gave -Ttstart/tstop/tinc[+n] */
204 double min, max, inc;
205 char *c = NULL;
206 min = atof (txt_a); max = atof (txt_b); inc = atof (txt_c);
207 if ((c = strrchr (txt_c, '+')) && (c[1] == 'n' || c[1] == '\0')) /* Gave number of points instead; calculate inc */
208 inc = (max - min) / (inc - 1.0);
209 if (inc <= 0.0) {
210 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Age increment must be positive\n");
211 n_errors++;
212 }
213 else {
214 Ctrl->T.n_times = lrint ((max - min) / inc) + 1;
215 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
216 for (t = 0; t < Ctrl->T.n_times; t++) Ctrl->T.value[t] = (t == (Ctrl->T.n_times-1)) ? max: min + t * inc;
217 }
218 }
219 else { /* Got a single time */
220 Ctrl->T.n_times = 1;
221 Ctrl->T.value = gmt_M_memory (GMT, NULL, Ctrl->T.n_times, double);
222 Ctrl->T.value[0] = atof (txt_a);
223 }
224 break;
225
226 case 'W': /* Use weights */
227 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
228 Ctrl->W.active = true;
229 break;
230 case 'Z': /* Report negative opening angles */
231 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
232 Ctrl->Z.active = true;
233 break;
234
235 default: /* Report bad options */
236 n_errors += gmt_default_error (GMT, opt->option);
237 break;
238 }
239 }
240 if (!Ctrl->A.active) n_in++; /* Got time in input column 3 */
241 if (Ctrl->W.active) n_in++; /* Got weights as extra column */
242 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = n_in;
243 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < n_in, "Binary input data (-bi) must have at least %u columns\n", n_in);
244 n_errors += gmt_M_check_condition (GMT, (Ctrl->N.active + Ctrl->S.active + Ctrl->W.active) > 1, "Only one of -N, -S, -Z can be set.\n");
245
246 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
247 }
248
rotsmoother_compare_ages(const void * point_1v,const void * point_2v)249 GMT_LOCAL int rotsmoother_compare_ages (const void *point_1v, const void *point_2v) {
250 struct ROTSMOOTHER_AGEROT *point_1, *point_2;
251
252 point_1 = (struct ROTSMOOTHER_AGEROT *)point_1v;
253 point_2 = (struct ROTSMOOTHER_AGEROT *)point_2v;
254 if (point_1->wxyasn[K_AGE] < point_2->wxyasn[K_AGE]) return (-1);
255 if (point_1->wxyasn[K_AGE] > point_2->wxyasn[K_AGE]) return (+1);
256 return (0);
257 }
258
259 #define bailout(code) {gmt_M_free_options (mode); return (code);}
260 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
261
GMT_rotsmoother(void * V_API,int mode,void * args)262 EXTERN_MSC int GMT_rotsmoother (void *V_API, int mode, void *args) {
263 bool stop;
264 uint64_t n_read = 0, rot, p, first = 0, last, n_use = 0, n_out = 0, n_total_use = 0, n_minimum, n_alloc = GMT_CHUNK;
265 int error = 0, n_fields;
266 unsigned int n_in = 3, k, j, t_col, w_col, t, n_cols = 4, matrix_dim = 3, nrots;
267 double *in = NULL, min_rot_angle, min_rot_age, max_rot_angle, max_rot_age;
268 double sum_rot_angle, sum_rot_angle2, sum_rot_age, sum_rot_age2, sum_weights;
269 double out[20], khat = 1.0, g = 1.0e-5; /* Common scale factor for all Covariance terms */
270 double sa2, ca2, qlength;
271 double z, mean_rot_angle, std_rot_angle, this_weight, x_comp, y_comp, med_angle;
272 double azimuth, norm, t_lo, t_hi, rot_angle_in_radians, this_rot_angle, lon_mean_pole, lat_mean_pole;
273 double x_in_plane[3], y_in_plane[3], C[9], EigenValue[3], EigenVector[9], work1[3], work2[3];
274 double R[3][3], DR[3][3], Ri[3][3], E[3], this_h[3], xyz_mean_pole[3], xyz_mean_quat[4], z_unit_vector[3];
275 double mean_rot_age, std_rot_age, *H[3], mean_H[3], Ccopy[9], *mangle = NULL, this_lon, this_lat;
276 struct ROTSMOOTHER_AGEROT *D = NULL;
277 struct GMT_RECORD *In = NULL, *Out = NULL;
278 struct ROTSMOOTHER_CTRL *Ctrl = NULL;
279 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
280 struct GMT_OPTION *options = NULL;
281 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
282
283 /*----------------------- Standard module initialization and parsing ----------------------*/
284
285 if (API == NULL) return (GMT_NOT_A_SESSION);
286 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
287 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
288
289 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
290
291 /* Parse the command-line arguments */
292
293 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 */
294 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
295 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
296 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
297
298 /*---------------------------- This is the rotsmoother main code ----------------------------*/
299
300 gmt_set_geographic (GMT, GMT_IN); /* In and out are lon/lat */
301 gmt_set_geographic (GMT, GMT_OUT); /* In and out are lon/lat */
302
303 /* Read the rotation data from file or stdin */
304
305 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
306
307 if (!Ctrl->A.active) n_in++; /* Got time */
308 if (Ctrl->W.active) n_in++; /* Got weights */
309 if (Ctrl->C.active) n_cols = 19; /* Want everything */
310
311 /* Specify input and output expected columns */
312 if ((error = GMT_Set_Columns (API, GMT_IN, n_in, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
313 Return (error);
314 }
315 /* Initialize the i/o for doing record-by-record reading/writing */
316 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
317 Return (API->error);
318 }
319 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
320 Return (API->error);
321 }
322 t_col = (Ctrl->A.active) ? GMT_Z : 3; /* If no time we use angle as proxy for time */
323 w_col = t_col + 1;
324 D = (struct ROTSMOOTHER_AGEROT *) gmt_M_memory (GMT, NULL, n_alloc, struct ROTSMOOTHER_AGEROT);
325 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
326
327 do { /* Keep returning records until we reach EOF */
328 if ((In = GMT_Get_Record (API, GMT_READ_DATA, &n_fields)) == NULL) { /* Read next record, get NULL if special case */
329 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
330 gmt_M_free (GMT, D);
331 gmt_M_free (GMT, Out);
332 Return (GMT_RUNTIME_ERROR);
333 }
334 if (gmt_M_rec_is_table_header (GMT)) { /* Skip all table headers */
335 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, NULL);
336 continue;
337 }
338 if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
339 break;
340 else if (gmt_M_rec_is_new_segment (GMT)) { /* Parse segment headers */
341 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
342 continue;
343 }
344 }
345 if (In->data == NULL) {
346 gmt_quit_bad_record (API, In);
347 Return (API->error);
348 }
349
350 in = In->data; /* Only need to process numerical part here */
351
352 /* Convert to geocentric, load parameters */
353 D[n_read].wxyasn[K_LON] = in[GMT_X];
354 D[n_read].wxyasn[K_LAT] = gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O);
355 D[n_read].wxyasn[K_ANGLE] = in[GMT_Z];
356 D[n_read].wxyasn[K_AGE] = in[t_col];
357 D[n_read].wxyasn[K_WEIGHT] = (Ctrl->W.active) ? in[w_col] : 1.0; /* Optionally use weights */
358 n_read++;
359 if (n_read == n_alloc) { /* Need larger arrays */
360 n_alloc <<= 1;
361 D = gmt_M_memory (GMT, D, n_alloc, struct ROTSMOOTHER_AGEROT);
362 }
363 } while (true);
364
365 if (n_read < n_alloc) D = gmt_M_memory (GMT, D, n_read, struct ROTSMOOTHER_AGEROT);
366
367 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
368 gmt_M_free (GMT, D);
369 gmt_M_free (GMT, Out);
370 Return (API->error);
371 }
372
373 min_rot_angle = min_rot_age = DBL_MAX; max_rot_angle = max_rot_age = -DBL_MAX;
374 for (rot = 0; rot < n_read; rot++) { /* Process rotations and make pseudo-vectors and quaternions */
375 if (D[rot].wxyasn[K_ANGLE] < min_rot_angle) min_rot_angle = D[rot].wxyasn[K_ANGLE];
376 if (D[rot].wxyasn[K_ANGLE] > max_rot_angle) max_rot_angle = D[rot].wxyasn[K_ANGLE];
377 if (D[rot].wxyasn[K_AGE] < min_rot_age) min_rot_age = D[rot].wxyasn[K_AGE];
378 if (D[rot].wxyasn[K_AGE] > max_rot_age) max_rot_age = D[rot].wxyasn[K_AGE];
379 gmt_geo_to_cart (GMT, D[rot].wxyasn[K_LAT], D[rot].wxyasn[K_LON], D[rot].P, true); /* Get unit vector on sphere */
380 rot_angle_in_radians = D[rot].wxyasn[K_ANGLE] * D2R; /* pseudo-vector length in radians */
381 /* Get quarternion by scaling the unit vector by sine of half the opening angle */
382 sincos (0.5 * rot_angle_in_radians, &sa2, &ca2);
383 D[rot].Q[0] = ca2; /* Store quarternion constant in first slot. */
384 for (k = 0; k < 3; k++) { /* Compute both pseudo-vector and quarternion components */
385 D[rot].Q[k+1] = D[rot].P[k] * sa2; /* Quarternion representation */
386 D[rot].P[k] *= rot_angle_in_radians; /* Pseudo-vector representation */
387 }
388 }
389
390 if (!Ctrl->A.active) GMT_Report (API, GMT_MSG_INFORMATION, "Range of input ages = %g/%g\n", min_rot_age, max_rot_age);
391 GMT_Report (API, GMT_MSG_INFORMATION, "Range of input angles = %g/%g\n", min_rot_angle, max_rot_angle);
392
393 if ((error = GMT_Set_Columns (API, GMT_OUT, n_cols, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
394 gmt_M_free (GMT, D);
395 gmt_M_free (GMT, Out);
396 Return (error);
397 }
398 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
399 gmt_M_free (GMT, D);
400 gmt_M_free (GMT, Out);
401 Return (API->error);
402 }
403
404 /* Sort the entire dataset on increasing ages */
405
406 qsort (D, n_read, sizeof (struct ROTSMOOTHER_AGEROT), rotsmoother_compare_ages);
407
408 if (GMT->common.h.add_colnames) { /* Create meaningful column header */
409 static char *short_header = "lon\tlat\ttime\tangle";
410 static char *long_header = "lon\tlat\ttime\tangle\tk_hat\ta\tb\tc\td\te\tf\tg\tdf\tstd_t\tstd_w\taz\tS1\tS2\tS3";
411 char *header = (Ctrl->C.active) ? long_header : short_header;
412 if (GMT_Set_Comment (API, GMT_IS_DATASET, GMT_COMMENT_IS_COLNAMES, header, NULL)) {
413 gmt_M_free (GMT, D);
414 gmt_M_free (GMT, Out);
415 Return (API->error);
416 }
417 }
418 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
419 gmt_M_free (GMT, D);
420 gmt_M_free (GMT, Out);
421 Return (API->error);
422 }
423 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
424 gmt_M_free (GMT, D);
425 gmt_M_free (GMT, Out);
426 Return (API->error);
427 }
428
429 z_unit_vector[0] = z_unit_vector[1] = 0.0; z_unit_vector[2] = 1.0; /* The local z unit vector */
430 n_minimum = (Ctrl->C.active) ? 2 : 1; /* Need at least two rotations to compute covariance, at least one to report the mean */
431
432 for (t = 1; t < Ctrl->T.n_times; t++) { /* For each desired output time interval */
433 t_lo = Ctrl->T.value[t-1]; t_hi = Ctrl->T.value[t]; /* The current interval */
434 for (rot = first, stop = false; !stop && rot < n_read; rot++) /* Determine index of first rotation inside this age window */
435 if (D[rot].wxyasn[K_AGE] >= t_lo) stop = true;
436 first = rot - 1; /* Index to first rotation inside this time interval */
437 for (stop = false; !stop && rot < n_read; rot++) /* Determine index of last rotation inside this age window */
438 if (D[rot].wxyasn[K_AGE] > t_hi) stop = true;
439 last = rot - 1; /* Index to first rotation outside this time interval */
440 n_use = last - first; /* Number of rotations in the interval */
441 GMT_Report (API, GMT_MSG_INFORMATION, "Found %d rots for the time interval %g <= t < %g\n", n_use, t_lo, t_hi);
442 if (n_use < n_minimum) continue; /* Need at least 1 or 2 poles to do anything useful */
443
444 /* Now estimate the average rotation */
445
446 gmt_M_memset (xyz_mean_pole, 3, double); /* Reset sum of mean components and weight sums */
447 gmt_M_memset (xyz_mean_quat, 4, double); /* Reset sum of mean components and weight sums */
448 sum_rot_angle = sum_rot_angle2 = sum_rot_age = sum_rot_age2 = sum_weights = 0.0;
449 for (rot = first; rot < last; rot++) { /* Loop over every rotation in this time interval... */
450 this_weight = D[rot].wxyasn[K_WEIGHT]; /* Either use weights or it is already set to 1 */
451 for (k = 0; k < 4; k++) xyz_mean_quat[k] += (this_weight * D[rot].Q[k]); /* Sum up the [weighted] quarternion components */
452 for (k = 0; k < 3; k++) xyz_mean_pole[k] += (this_weight * D[rot].P[k]); /* Sum up the [weighted] Cartesian components */
453 z = this_weight * D[rot].wxyasn[K_ANGLE]; /* Separately sum up opening [weighted] angle stats */
454 sum_rot_angle += z;
455 sum_rot_angle2 += (z * D[rot].wxyasn[K_ANGLE]);
456 sum_weights += this_weight;
457 z = this_weight * D[rot].wxyasn[K_AGE]; /* Separately sum up opening [weighted] age stats */
458 sum_rot_age += z;
459 sum_rot_age2 += (z * D[rot].wxyasn[K_AGE]);
460 }
461
462 /* Get [weighted] mean age */
463 mean_rot_age = sum_rot_age / sum_weights;
464 std_rot_age = sqrt ((sum_weights * sum_rot_age2 - sum_rot_age * sum_rot_age) / (sum_weights * sum_weights * ((n_use - 1.0) / n_use)));
465 /* Get [weighted] mean angle */
466 mean_rot_angle = sum_rot_angle / sum_weights;
467 std_rot_angle = sqrt ((sum_weights * sum_rot_angle2 - sum_rot_angle * sum_rot_angle) / (sum_weights * sum_weights * ((n_use - 1.0) / n_use)));
468
469 /* Get Euclidean [weighted] mean pseudo-vector location */
470 for (k = 0; k < 3; k++) xyz_mean_pole[k] /= sum_weights; /* Components of [weighted] mean pseudovector */
471 gmt_normalize3v (GMT, xyz_mean_pole); /* Make it a unit vector ... */
472
473 /* Get [weighted] mean quaternion */
474 for (k = 0; k < 4; k++) xyz_mean_quat[k] /= sum_weights; /* Components of mean quaternion */
475 for (k = 0, qlength = 0.0; k < 4; k++) qlength += xyz_mean_quat[k] * xyz_mean_quat[k]; /* Sum squared components of total quaternion */
476 qlength = sqrt (qlength); /* Length of the quaternion */
477 for (k = 0; k < 4; k++) xyz_mean_quat[k] /= qlength; /* Components of mean unit beast */
478 mean_rot_angle = 2.0 * d_acos (xyz_mean_quat[0]) * R2D; /* Rotation (in degrees) represented by mean quaternion */
479 norm = d_sqrt (1.0 - xyz_mean_quat[0] * xyz_mean_quat[0]); /* This is sin (angle/2) */
480 for (k = 1; k < 4; k++) xyz_mean_quat[k] /= norm; /* Unit vector of mean quaternion */
481 gmt_cart_to_geo (GMT, &lat_mean_pole, &lon_mean_pole, &xyz_mean_quat[1], true); /* Convert it to lon/lat pole */
482 mangle = (double *) gmt_M_memory (GMT, NULL, n_use, double); /* For median angle calculation we need to make an array */
483 /* Now compute the median of the opening angle */
484 for (rot = first, p = 0; rot < last; rot++, p++) /* For each angle of rotation in this rotation interval... */
485 mangle[p] = D[rot].wxyasn[K_ANGLE]; /* Place angles in temp array for computing median angle */
486 gmt_median (GMT, mangle, n_use, min_rot_angle, max_rot_angle, 0.5 * (min_rot_angle + max_rot_angle), &med_angle);
487 gmt_M_free (GMT, mangle);
488
489 if ((Ctrl->N.active && mean_rot_angle > 0.0) || (Ctrl->N.active && lat_mean_pole < 0.0) || (!Ctrl->N.active && lat_mean_pole > 0.0)) {
490 /* Depending on -S, -N, -W. flip pole to right hemisphere or sign of angle */
491 lat_mean_pole = -lat_mean_pole;
492 lon_mean_pole += 180.0;
493 mean_rot_angle = -mean_rot_angle;
494 while (lon_mean_pole > 180.0) lon_mean_pole -= 360.0;
495 }
496 /* Prepare output array for the mean rotation */
497 out[GMT_X] = lon_mean_pole;
498 out[GMT_Y] = gmt_lat_swap (GMT, lat_mean_pole, GMT_LATSWAP_O2G); /* Convert back from geocentric to geodetic */
499 out[GMT_Z] = mean_rot_age; out[3] = mean_rot_angle;
500 GMT_Report (API, GMT_MSG_INFORMATION, "Mean opening angle %8.4f vs median opening angle %8.4f\n", mean_rot_angle, med_angle); /* Report mean and median angle for testing */
501 n_total_use += n_use;
502 n_out++;
503
504 if (!Ctrl->C.active) { /* No covariance requested, print this rotation and continue */
505 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
506 continue;
507 }
508
509 /* Now get the covariance matrix */
510
511 /* Here we also want to get the covariance matrix. To do so we want to parametrize all
512 * the n_use rotations as a small incremental rotation DR_i followed by the mean rotation
513 * R, i.e., R_i = R * DR_i. We must therefore first solve for the incremental rotation
514 * DR_i (a 3x3 matrix) which becomes R^T * R_i and then convert it to a pseudo-vector h
515 * so that we can calculate covariances. Because the incremental rotation angles are
516 * very small (|pseudo-vector| << 1 the pseudo-vector approach is warranted. */
517
518 for (k = 0; k < 3; k++) H[k] = gmt_M_memory (GMT, NULL, n_use, double);
519 gmt_M_memset (mean_H, 3, double); /* Zero out the mean sum matrix for pseudovectors */
520 gmt_make_rot_matrix (GMT, lon_mean_pole, lat_mean_pole, -mean_rot_angle, R); /* Mean rotation R; use -ve angle since we actually need R^T */
521 for (rot = first, p = 0; rot < last; rot++, p++) { /* For each individually determined pole of rotation in this rotation interval... */
522 gmt_M_memcpy (E, D[rot].P, 3, double); /* Set up the rotation pseudo-vector E */
523 this_rot_angle = gmt_mag3v (GMT, E) * R2D; /* Magnitude of this vector -> opening in degrees */
524 gmt_normalize3v (GMT, E); /* Make rotation pole a unit vector */
525 gmt_make_rot_matrix2 (GMT, E, this_rot_angle, Ri); /* Compute this individual rotation matrix R_i */
526 spotter_matrix_mult (GMT, R, Ri, DR); /* The incremental rotation DR_i = R^T * R_i */
527 spotter_matrix_to_pole (GMT, DR, &this_lon, &this_lat, &this_rot_angle); /* Get (lon, lat, angle) of the incremental rotation */
528 gmt_geo_to_cart (GMT, this_lat, this_lon, this_h, true); /* Make unit vector this_h for incremental rotation pole */
529 rot_angle_in_radians = this_rot_angle * D2R; /* Scale to get pseudo-vector with length in radians */
530 for (k = 0; k < 3; k++) { /* For each x,y,z component of pseudo-vector this_h */
531 H[k][p] = this_h[k] * rot_angle_in_radians; /* Final pseudo-vector for the rotation increment */
532 mean_H[k] += H[k][p]; /* Sum of x,y,z components of all such pseudo-vectors */
533 }
534 }
535
536 gmt_M_memset (C, 9, double); /* Zero out the covariance matrix */
537 for (k = 0; k < 3; k++) mean_H[k] /= n_use; /* The mean pseudo-vector */
538 for (p = 0; p < n_use; p++) { /* Calculate the covariance matrix */
539 for (k = 0; k < 3; k++) H[k][p] -= mean_H[k]; /* Remove mean component from each vector */
540 for (j = 0; j < 3; j++) for (k = 0; k < 3; k++) {
541 C[3*j+k] += (H[j][p]) * (H[k][p]); /* Products of deviations from each mean */
542 }
543 }
544 for (k = 0; k < 9; k++) C[k] /= (n_use - 1.0); /* Normalize to get final var-covar matrix, C */
545 for (k = 0; k < 3; k++) gmt_M_free (GMT, H[k]); /* Release temporary memory */
546
547 /* We also want to get a 95% ellipse projected onto a tangent plane to the surface of Earth at the pole */
548
549 gmt_M_memcpy (Ccopy, C, 9, double); /* Duplicate since it will get trodden on */
550 if (gmt_jacobi (GMT, Ccopy, matrix_dim, matrix_dim, EigenValue, EigenVector, work1, work2, &nrots)) { /* Solve eigen-system C = EigenVector * EigenValue * EigenVector^T */
551 GMT_Report (API, GMT_MSG_ERROR, "Eigenvalue routine gmt_jacobi failed to converge in 50 sweeps.\n");
552 }
553
554 /* In addition to reporting the covariance, we will report the azimuth of the major axis and the
555 * lengths of the major and minor axes as they appear when projected onto the local plane tangent
556 * to the surface at the mean pole. Note that EigenVector[0], EigenVector[1], EigenVector[2] are
557 * the components of the 1st eigenvector. */
558
559 gmt_cross3v (GMT, z_unit_vector, xyz_mean_pole, x_in_plane); /* Local x-axis in plane normal to mean pole */
560 gmt_cross3v (GMT, xyz_mean_pole, x_in_plane, y_in_plane); /* Local y-axis in plane normal to mean pole */
561 x_comp = gmt_dot3v (GMT, EigenVector, x_in_plane); /* x-component of major axis in tangent plane */
562 y_comp = gmt_dot3v (GMT, EigenVector, y_in_plane); /* y-component of major axis in tangent plane */
563 azimuth = fmod (360.0 + (90.0 - atan2 (y_comp, x_comp) * R2D), 360.0); /* Azimuth of major axis */
564 if (azimuth > 180.0) azimuth -= 180.0;
565 /* Get the three axes lengths from the eigenvalues scaled from radians to km */
566 EigenValue[0] = sqrt (EigenValue[0]/khat) * EQ_RAD * SQRT_CHI2;
567 EigenValue[1] = sqrt (EigenValue[1]/khat) * EQ_RAD * SQRT_CHI2;
568 EigenValue[2] = sqrt (EigenValue[2]/khat) * EQ_RAD * SQRT_CHI2;
569
570 /* Report the results. It is conventional to define the covariance matrix as follows:
571 *
572 * [ a b d ]
573 * C = cov (u) = k^(-1) * | b c e | * g
574 * [ d e f ]
575 *
576 * where k is a scale factor that relates true estimate of uncertainty in the data to the assigned
577 * estimate. For this work we have no uncertainty estimate so k == 1. a-f are six unique entries in the
578 * symmetric C matrix, and g is just a common scale factor to avoid too many decimal places in the
579 * tables; typically g is 1e-5 ,which we adopt here. The output is thus written in the order k, a-g,
580 * followed by d.f., the degrees of freedom.
581 */
582
583 for (k = 0; k < 9; k++) C[k] /= g; /* Take out the common factor */
584 /* Place khat a b c d e f g df */
585 out[4] = khat;
586 out[5] = C[0];
587 out[6] = C[1];
588 out[7] = C[4];
589 out[8] = C[2];
590 out[9] = C[5];
591 out[10] = C[8];
592 out[11] = g;
593 out[12] = (double)n_use - 1;
594 out[13] = std_rot_age; /* Add std_age and std_angle */
595 out[14] = std_rot_angle;
596 out[15] = azimuth; /* Axes and orientation for surface 95% error ellipse */
597 out[16] = EigenValue[0];
598 out[17] = EigenValue[1];
599 out[18] = EigenValue[2];
600 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
601 }
602
603 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
604 Return (API->error);
605 }
606
607 if (n_out == 0) {
608 GMT_Report (API, GMT_MSG_ERROR, "No smoothed poles created - is filterwidth too small?\n");
609 if (n_total_use != n_read) GMT_Report (API, GMT_MSG_ERROR, "Read %ld poles, used %7ld\n", n_read, n_total_use);
610 }
611
612 /* Free up memory used for the array */
613
614 gmt_M_free (GMT, D);
615 gmt_M_free (GMT, Out);
616
617 Return (GMT_NOERROR);
618 }
619