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 */
18
19 #include "gmt_dev.h"
20 #include "spotter.h"
21
22 #define THIS_MODULE_CLASSIC_NAME "polespotter"
23 #define THIS_MODULE_MODERN_NAME "polespotter"
24 #define THIS_MODULE_LIB "spotter"
25 #define THIS_MODULE_PURPOSE "Find stage poles given fracture zones and abyssal hills"
26 #define THIS_MODULE_KEYS "AD(,CD),FD(,GG},LD)"
27 #define THIS_MODULE_NEEDS ""
28 #define THIS_MODULE_OPTIONS "-:>RVbdefghios" GMT_OPT("HMm")
29
30 #define RADIAN2KM (R2D * GMT->current.proj.DIST_KM_PR_DEG)
31
32 struct POLESPOTTER_CTRL { /* All control options for this program (except common args) */
33 /* active is true if the option has been activated */
34 struct POLESPOTTER_A { /* -A<abyssalhilefile> */
35 bool active;
36 char *file;
37 double weight;
38 } A;
39 struct POLESPOTTER_D { /* -D<spacing> */
40 bool active;
41 double length;
42 } D;
43 struct POLESPOTTER_E { /* -Ea|f<sigma> */
44 bool active[2];
45 } E;
46 struct POLESPOTTER_F { /* -F<fzfile> */
47 bool active;
48 char *file;
49 double weight;
50 } F;
51 struct POLESPOTTER_G { /* -Goutfile */
52 bool active;
53 char *file;
54 } G;
55 struct POLESPOTTER_I { /* -I (for checking only) */
56 bool active;
57 } I;
58 struct POLESPOTTER_N { /* -N */
59 bool active;
60 } N;
61 struct POLESPOTTER_S { /* -Ss|l|p[<modifiers>] */
62 bool active;
63 bool dump_lines;
64 bool dump_crossings;
65 bool midpoint;
66 unsigned int mode;
67 char *file;
68 double plon, plat;
69 } S;
70 };
71
72 enum spotter_modes {
73 SPOTTER_SCAN_SPOTS = 0,
74 SPOTTER_SCAN_LINES = 1,
75 SPOTTER_SCAN_POLES = 2};
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 POLESPOTTER_CTRL *C;
79
80 C = gmt_M_memory (GMT, NULL, 1, struct POLESPOTTER_CTRL);
81
82 /* Initialize values whose defaults are not 0/false/NULL */
83 C->A.weight = C->F.weight = 1.0;
84 C->D.length = 5.0; /* 5 km spacing */
85 return (C);
86 }
87
Free_Ctrl(struct GMT_CTRL * GMT,struct POLESPOTTER_CTRL * C)88 static void Free_Ctrl (struct GMT_CTRL *GMT, struct POLESPOTTER_CTRL *C) { /* Deallocate control structure */
89 if (!C) return;
90 gmt_M_str_free (C->A.file);
91 gmt_M_str_free (C->F.file);
92 gmt_M_str_free (C->S.file);
93 gmt_M_str_free (C->G.file);
94 gmt_M_free (GMT, C);
95 }
96
usage(struct GMTAPI_CTRL * API,int level)97 static int usage (struct GMTAPI_CTRL *API, int level) {
98 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
99 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
100 GMT_Usage (API, 0, "usage: %s [-A<abyssalhills>] [-D<step>] [-Ea|f<sigma>] [-F<FZfile>] [-G%s] "
101 "[%s] [-N] [%s] [-Ss|p|l[<modifiers>]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
102 name, GMT_OUTGRID, GMT_Id_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_d_OPT, GMT_e_OPT, GMT_h_OPT,
103 GMT_i_OPT, GMT_r_OPT, GMT_o_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
104
105 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
106
107 GMT_Message (API, GMT_TIME_NONE, " OPTIONAL ARGUMENTS:\n");
108 GMT_Usage (API, 1, "\n-A<abyssalhills>");
109 GMT_Usage (API, -2, "Give multisegment file with abyssal hill lineaments [none].");
110 GMT_Usage (API, 1, "\n-D<step>");
111 GMT_Usage (API, -2, "Give step-length along great circles in km [5].");
112 GMT_Usage (API, 1, "\n-Ea|f<sigma>");
113 GMT_Usage (API, -2, "Specify typical angular error (in degrees) for (a)byssal hills or (f)racture zones [1]. Repeatable");
114 GMT_Usage (API, 1, "\n-F<FZfile>");
115 GMT_Usage (API, -2, "Give multisegment file with fracture zone lineaments [none].");
116 gmt_outgrid_syntax (API, 'G', "Specify file name for output pole-density grid [no grid]. Requires -R -I [-r]. "
117 "Accumulates weighted great-circle length density on the grid");
118 GMT_Usage (API, 1, "\n%s", GMT_Id_OPT);
119 GMT_Usage (API, -2, "Specify grid interval(s); Append m [or s] to <dx> and/or <dy> for minutes [or seconds].");
120 GMT_Usage (API, 1, "\n-N Normalize grid so maximum is 1 [no normalization].");
121 GMT_Option (API, "Rg");
122 GMT_Usage (API, 1, "\n-Ss|p|l[<modifiers>]");
123 GMT_Usage (API, -2, "Set the spotter directive:");
124 GMT_Usage (API, 3, "s: Scan for spots [Default]. This mode offers two optional modifiers:");
125 GMT_Usage (API, 4, "+l Dump all great circles to standard output [none].");
126 GMT_Usage (API, 4, "+c Save all great circle intersections to appended file <xfile> [no crossings].");
127 GMT_Usage (API, 3, "p: Scan for poles. Writes a misfit grid to <outgrid>.");
128 GMT_Usage (API, 3, "l: Scan for compatible lines given appended <plon>/<plat> trial pole. "
129 "Append +m to report misfit for each midpoint.");
130 GMT_Option (API, "V");
131 GMT_Option (API, "bi2,d,e,h,i,o,r,s,:,.");
132
133 return (GMT_MODULE_USAGE);
134 }
135
parse(struct GMT_CTRL * GMT,struct POLESPOTTER_CTRL * Ctrl,struct GMT_OPTION * options)136 static int parse (struct GMT_CTRL *GMT, struct POLESPOTTER_CTRL *Ctrl, struct GMT_OPTION *options) {
137 /* This parses the options provided to polespotter and sets parameters in CTRL.
138 * Any GMT common options will override values set previously by other commands.
139 * It also replaces any file names specified as input or output with the data ID
140 * returned when registering these sources/destinations with the API.
141 */
142
143 unsigned int n_errors = 0;
144 int n;
145 char txt_a[GMT_LEN256] = {""}, txt_b[GMT_LEN256] = {""}, *c = NULL;
146 struct GMT_OPTION *opt = NULL;
147 struct GMTAPI_CTRL *API = GMT->parent;
148
149 for (opt = options; opt; opt = opt->next) {
150 switch (opt->option) {
151
152 /* Supplemental parameters */
153
154 case 'A': /* File with abyssal hill traces */
155 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
156 Ctrl->A.active = true;
157 if (opt->arg[0]) Ctrl->A.file = strdup (opt->arg);
158 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->A.file))) n_errors++;
159 break;
160 case 'D': /* Step length */
161 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
162 Ctrl->D.active = true;
163 Ctrl->D.length = atof (opt->arg);
164 break;
165 case 'E': /* Sigma of lines (store 1/sigma here) */
166 switch (opt->arg[0]) {
167 case 'a':
168 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active[0]);
169 Ctrl->E.active[0] = true;
170 Ctrl->A.weight = 1.0 / atof (&opt->arg[1]);
171 break;
172 case 'f':
173 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active[1]);
174 Ctrl->E.active[1] = true;
175 Ctrl->F.weight = 1.0 / atof (&opt->arg[1]);
176 break;
177 }
178 break;
179 case 'F': /* File with fracture zone traces */
180 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
181 Ctrl->F.active = true;
182 if (opt->arg[0]) Ctrl->F.file = strdup (opt->arg);
183 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->F.file))) n_errors++;
184 break;
185 case 'G':
186 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
187 Ctrl->G.active = true;
188 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
189 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_REMOTE, &(Ctrl->G.file))) n_errors++;
190 break;
191 case 'I':
192 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
193 Ctrl->I.active = true;
194 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
195 break;
196 case 'N': /* Normalize grid */
197 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
198 Ctrl->N.active = true;
199 break;
200 case 'S': /* modes */
201 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
202 switch (opt->arg[0]) {
203 case 'p': Ctrl->S.mode = SPOTTER_SCAN_POLES;
204 break;
205 case 'l': Ctrl->S.mode = SPOTTER_SCAN_LINES;
206 if (gmt_validate_modifiers (GMT, opt->arg, 'S', "m", GMT_MSG_ERROR)) n_errors++;
207 if ((c = strstr (opt->arg, "+m"))) { /* Do midpoint analysis instead */
208 Ctrl->S.midpoint = true;
209 c[0] = '\0'; /* Chop off modifier */
210 }
211 if ((n = sscanf (&opt->arg[1], "%[^/]/%s", txt_a, txt_b)) != 2) {
212 GMT_Report (API, GMT_MSG_ERROR, "Option -Sp: No pole given\n");
213 n_errors++;
214 }
215 else {
216 n_errors += gmt_verify_expectations (GMT, GMT_IS_LON, gmt_scanf_arg (GMT, txt_a, GMT_IS_LON, false, &Ctrl->S.plon), txt_a);
217 n_errors += gmt_verify_expectations (GMT, GMT_IS_LAT, gmt_scanf_arg (GMT, txt_b, GMT_IS_LAT, false, &Ctrl->S.plat), txt_b);
218 }
219 if (c) c[0] = '+'; /* Restore modifier */
220 break;
221 case 's': Ctrl->S.mode = SPOTTER_SCAN_SPOTS;
222 if (gmt_validate_modifiers (GMT, opt->arg, 'S', "cl", GMT_MSG_ERROR)) n_errors++;
223 if (gmt_get_modifier (opt->arg, 'l', txt_a)) /* Dump lines to stdout */
224 Ctrl->S.dump_lines = true;
225 if (gmt_get_modifier (opt->arg, 'c', txt_a) && txt_a[0]) { /* Crossing output file */
226 Ctrl->S.dump_crossings = true;
227 Ctrl->S.file = strdup (txt_a);
228
229 }
230 break;
231 default:
232 n_errors++;
233 }
234 Ctrl->S.active = true;
235 break;
236
237 default: /* Report bad options */
238 n_errors += gmt_default_error (GMT, opt->option);
239 break;
240 }
241 }
242
243 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0) GMT->common.b.ncol[GMT_IN] = 2;
244 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] < 2, "Binary input data (-bi) must have at least 3 columns\n");
245 n_errors += gmt_M_check_condition (GMT, !Ctrl->A.active && !Ctrl->F.active, "At least one of -A or -F is required.\n");
246 if (Ctrl->G.active) {
247 n_errors += gmt_M_check_condition (GMT, Ctrl->G.file == NULL, "Option -G: Must specify output file\n");
248 n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0, "Option -I: Must specify positive increment(s)\n");
249 }
250 n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->D.length <= 0.0, "Option -D: Must specify a positive length step.\n");
251 n_errors += gmt_M_check_condition (GMT, Ctrl->S.mode == SPOTTER_SCAN_SPOTS && !Ctrl->S.dump_lines && !Ctrl->G.active, "Option -Ss: Must specify at least one of -G, -L.\n");
252 n_errors += gmt_M_check_condition (GMT, Ctrl->S.dump_crossings && !Ctrl->S.file, "Option -Ss: Must specify a file name if +c is used.\n");
253
254 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
255 }
256
257 #define POLESPOTTER_AH 0
258 #define POLESPOTTER_FZ 1
259
260 EXTERN_MSC void gmtlib_load_rot_matrix (double w, double R[3][3], double E[]);
261 EXTERN_MSC void gmtlib_init_rot_matrix (double R[3][3], double E[]);
262
polespotter_get_cross_normalized(struct GMT_CTRL * GMT,double P1[],double P2[],double G[])263 GMT_LOCAL void polespotter_get_cross_normalized (struct GMT_CTRL *GMT, double P1[], double P2[], double G[]) {
264 gmt_cross3v (GMT, P1, P2, G); /* G is the pole of the great circle that passes through P1 & P2 */
265 gmt_normalize3v (GMT, G); /* But we need to normalize it */
266 }
267
polespotter_get_great_circle_pole(struct GMT_CTRL * GMT,double P1[],double P2[],unsigned int type,double M[],double G[])268 GMT_LOCAL void polespotter_get_great_circle_pole (struct GMT_CTRL *GMT, double P1[], double P2[], unsigned int type, double M[], double G[]) {
269 /* Input is P1 and P2, the two cartesian points defining a small great circle segment.
270 * Output is M (the mid point of the segment) and G, the great circle defining the bisector (FZ) or segment itself (AH) */
271 unsigned int k;
272 for (k = 0; k < 3; k++) M[k] = 0.5 * (P1[k] + P2[k]); /* Mid-point M */
273 gmt_normalize3v (GMT, M); /* To get a unit length vector M */
274 polespotter_get_cross_normalized (GMT, P1, P2, G); /* G is the pole of the great circle that passes through P1 & P2 */
275 if (type == POLESPOTTER_FZ) { /* Must get the bisector pole instead, so cross it with M */
276 double B[3]; /* Temp vector */
277 polespotter_get_cross_normalized (GMT, M, G, B); /* This gives the normalied bisector pole instead */
278 gmt_M_memcpy (G, B, 3, double); /* Put bisector pole into G which is what we return */
279 }
280 }
281
polespotter_get_angle_between_trends(struct GMT_CTRL * GMT,double P1[],double P2[],unsigned int type,double X[])282 GMT_LOCAL double polespotter_get_angle_between_trends (struct GMT_CTRL *GMT, double P1[], double P2[], unsigned int type, double X[]) {
283 /* P1 and P2 are two points on a FZ or AH. Midpoint of P1 and P2 is M. X is a trial pole.
284 * If type = FZ: Find difference in orientations between bisector to P1-P2 and great circle from X through M.
285 * If type = AH: Find difference in orientations between great circle through P1-P2 and great circle from X through M.
286 */
287 double M[3], G[3], B[3], cos_del_angle, del_angle;
288 polespotter_get_great_circle_pole (GMT, P1, P2, type, M, G); /* Obtain great circle pole to segment (or bisector if FZ) as well as mid-point M */
289 polespotter_get_cross_normalized (GMT, M, X, B); /* B is pole for great circle through selected pole and M */
290 cos_del_angle = gmt_dot3v (GMT, G, B); /* Cos of angle between great circles is our other weight */
291 del_angle = fabs (d_acos (cos_del_angle)); /* Get |angle| between the two great circles */
292 if (del_angle > M_PI_2) del_angle = M_PI - del_angle; /* Since angles are actually orientation differences */
293 return (del_angle);
294 }
295
296 #define bailout(code) {gmt_M_free_options (mode); return (code);}
297 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
298
GMT_polespotter(void * V_API,int mode,void * args)299 EXTERN_MSC int GMT_polespotter (void *V_API, int mode, void *args) {
300 bool create_great_circles;
301 int error;
302 unsigned int d, n_steps, grow, gcol, k;
303 uint64_t node, tbl, seg, row, ng = 0;
304 size_t n_alloc = 0;
305 char header[GMT_LEN128] = {""}, *code = NULL;
306 const char *label[2] = {"AH", "FZ"};
307 gmt_grdfloat *layer = NULL;
308 double weight, seg_weight, angle_radians, d_angle_radians, mlon, mlat, glon, glat, L, in[2];
309 double P1[3], P2[3], M[3], G[3], X[3];
310
311 struct GMT_OPTION *ptr = NULL;
312 struct GMT_GRID *Grid = NULL;
313 struct GMT_DATASET *In[2] = {NULL, NULL};
314 struct GMT_DATASEGMENT *S = NULL;
315 struct GMT_RECORD *Out = NULL;
316 struct POLESPOTTER_CTRL *Ctrl = NULL;
317 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
318 struct GMT_OPTION *options = NULL;
319 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
320
321 /*----------------------- Standard module initialization and parsing ----------------------*/
322
323 if (API == NULL) return (GMT_NOT_A_SESSION);
324 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
325 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
326
327 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
328
329 /* Parse the command-line arguments */
330
331 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 */
332 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
333 if ((ptr = GMT_Find_Option (API, 'f', options)) == NULL) gmt_parse_common_options (GMT, "f", 'f', "g"); /* Did not set -f, implicitly set -fg */
334 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
335 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
336
337 /*---------------------------- This is the polespotter main code ----------------------------*/
338
339 if (Ctrl->A.active && (In[POLESPOTTER_AH] = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->A.file, NULL)) == NULL) {
340 GMT_Report (API, GMT_MSG_ERROR, "Unable to open file with abyssal hill lineaments: %s", Ctrl->A.file);
341 Return (API->error);
342 }
343 if (Ctrl->F.active && (In[POLESPOTTER_FZ] = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->F.file, NULL)) == NULL) {
344 GMT_Report (API, GMT_MSG_ERROR, "Unable to open file with fracture zone lineaments: %s", Ctrl->F.file);
345 Return (API->error);
346 }
347
348 /* Initialize the pole search grid and structure */
349
350 if (Ctrl->G.active) {
351 if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
352 GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
353 if (Ctrl->S.mode == SPOTTER_SCAN_SPOTS)
354 layer = gmt_M_memory_aligned (GMT, NULL, Grid->header->size, gmt_grdfloat);
355 }
356
357 /* Determine how often to sample the great circle given -D. Since all great cirlces are sampled this way
358 * we only do a cos(lat) weighting for the line density if -G is used. */
359 n_steps = urint (360.0 * GMT->current.proj.DIST_KM_PR_DEG / Ctrl->D.length);
360 d_angle_radians = TWO_PI / (n_steps - 1);
361
362 if (Ctrl->S.dump_lines) {
363 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
364 Return (API->error);
365 }
366 if ((error = GMT_Set_Columns (API, GMT_OUT, 2, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
367 Return (error);
368 }
369 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
370 Return (API->error);
371 }
372 Out = gmt_new_record (GMT, in, NULL); /* Since we only need to worry about numerics in this module */
373 }
374
375 if (Ctrl->S.mode == SPOTTER_SCAN_SPOTS) {
376 double Rot0[3][3], Rot[3][3], *GG = NULL;
377 /* Loop over all abyssal hill and fracture zone lines and consider each subsecutive pair of points to define a great circle that intersects the pole */
378
379 GMT_Report (API, GMT_MSG_INFORMATION, "Entering scan mode: spots\n");
380 create_great_circles = (Ctrl->G.active || Ctrl->S.dump_lines);
381 if (Ctrl->S.dump_crossings) { /* Need temporary storage for all great circle poles and their weights and type */
382 n_alloc = GMT_BIG_CHUNK;
383 GG = gmt_M_memory (GMT, NULL, n_alloc*4, double);
384 code = gmt_M_memory (GMT, NULL, n_alloc, char);
385 }
386 for (d = POLESPOTTER_AH; d <= POLESPOTTER_FZ; d++) {
387 if (In[d] == NULL) continue; /* Don't have this data set */
388 weight = (d == POLESPOTTER_AH) ? Ctrl->A.weight : Ctrl->F.weight;
389 for (tbl = 0; tbl < In[d]->n_tables; tbl++) {
390 for (seg = 0; seg < In[d]->table[tbl]->n_segments; seg++) { /* For each segment in the table */
391 S = In[d]->table[tbl]->segment[seg]; /* Set shortcut to current segment */
392 if (gmt_parse_segment_item (GMT, S->header, "-D", header)) /* Found -D<sigma> */
393 seg_weight = 1.0 / atof (header);
394 else
395 seg_weight = weight; /* Already got 1/sigma via -E, actually */
396 /* Convert the entire segment to geocentric latitude */
397 for (row = 0; row < S->n_rows; row++) S->data[GMT_Y][row] = gmt_lat_swap (GMT, S->data[GMT_Y][row], GMT_LATSWAP_G2O);
398 gmt_geo_to_cart (GMT, S->data[GMT_Y][0], S->data[GMT_X][0], P1, true); /* get x/y/z of first point P1 */
399 for (row = 1; row < S->n_rows; row++) {
400 if (Ctrl->G.active) gmt_M_memset (layer, Grid->header->size, gmt_grdfloat);
401 gmt_geo_to_cart (GMT, S->data[GMT_Y][row], S->data[GMT_X][row], P2, true); /* get x/y/z of 2nd point P2 */
402 L = d_acos (gmt_dot3v (GMT, P1, P2)) * RADIAN2KM * seg_weight; /* Weighted length of this segment */
403 polespotter_get_great_circle_pole (GMT, P1, P2, d, M, G); /* Obtain great circle pole to segment (or bisector if FZ) */
404 if (Ctrl->S.dump_crossings) { /* Keep track of great circles to each line */
405 gmt_M_memcpy (&GG[4*ng], G, 3, double);
406 GG[4*ng+3] = L;
407 code[ng++] = (char)d;
408 if (ng == n_alloc) {
409 n_alloc <<= 1;
410 GG = gmt_M_memory (GMT, GG, n_alloc*4, double);
411 code = gmt_M_memory (GMT, code, n_alloc, char);
412 }
413 }
414 if (create_great_circles) {
415 gmt_cart_to_geo (GMT, &mlat, &mlon, M, true); /* Get lon/lat of the mid point */
416 gmt_cart_to_geo (GMT, &glat, &glon, G, true); /* Get lon/lat of the mid point */
417 gmtlib_init_rot_matrix (Rot0, G); /* Get partial rotation matrix since no actual angle is applied yet */
418 sprintf (header, "Great circle: Center = %g/%g and pole is %g/%g -I%s", mlon, mlat, glon, glat, label[d]);
419 if (Ctrl->S.dump_lines) GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, header);
420 for (k = 0; k < n_steps; k++) {
421 angle_radians = k * d_angle_radians; /* The angle around the great circle (0-360) */
422 gmt_M_memcpy (Rot, Rot0, 9, double); /* Get a copy of the "0-angle" rotation matrix */
423 gmtlib_load_rot_matrix (angle_radians, Rot, G); /* Build the actual rotation matrix for this angle */
424 gmt_matrix_vect_mult (GMT, 3U, Rot, M, X); /* Rotate the mid point along the great circle */
425 gmt_cart_to_geo (GMT, &in[GMT_Y], &in[GMT_X], X, true); /* Get lon/lat of this point along crossing profile */
426 in[GMT_Y] = gmt_lat_swap (GMT, in[GMT_Y], GMT_LATSWAP_G2O + 1); /* Convert back to geodetic */
427 if (Ctrl->S.dump_lines) GMT_Put_Record (API, GMT_WRITE_DATA, Out); /* We are writing these circles to output */
428 if (!Ctrl->G.active) continue; /* Not doing density grid here */
429 if (gmt_M_y_is_outside (GMT, in[GMT_Y], Grid->header->wesn[YLO], Grid->header->wesn[YHI])) continue; /* Outside y-range */
430 if (gmt_x_is_outside (GMT, &in[GMT_X], Grid->header->wesn[XLO], Grid->header->wesn[XHI])) continue; /* Outside x-range (or periodic longitude) */
431 if (gmt_row_col_out_of_bounds (GMT, in, Grid->header, &grow, &gcol)) continue; /* Sorry, outside after all */
432 node = gmt_M_ijp (Grid->header, grow, gcol); /* Bin index */
433 layer[node] = (gmt_grdfloat)(cosd (in[GMT_Y]) * L); /* Any bin intersected will have this single value despite perhaps many intersections */
434 }
435 if (Ctrl->G.active) { /* Add density layer of this great circle to the total density grid */
436 for (node = 0; node < Grid->header->size; node++) Grid->data[node] += layer[node];
437 }
438 }
439 gmt_M_memcpy (P1, P2, 3, double); /* Let old P2 be next P1 */
440 }
441 }
442 }
443 }
444 if (Ctrl->G.active) gmt_M_free_aligned (GMT, layer);
445
446 if (Ctrl->S.dump_lines) { /* Disables further line output */
447 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further line output */
448 Return (API->error);
449 }
450 gmt_M_free (GMT, Out);
451 }
452
453 if (Ctrl->S.dump_crossings) { /* Generate crossings and write them to the crossings file */
454 uint64_t dim[4] = {1, 1, 0, 5}, n_cross, g1, g2;
455 struct GMT_DATASET *C = NULL;
456 struct GMT_DATASEGMENT *S = NULL;
457 n_cross = ng * (ng - 1);
458 dim[GMT_ROW] = n_cross;
459 if ((C = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL)
460 Return (API->error);
461 S = C->table[0]->segment[0]; /* Only have a single segment here*/
462 for (g1 = k = 0; g1 < ng; g1++) {
463 for (g2 = g1+1; g2 < ng; g2++) { /* Get circle intersections (2) */
464 polespotter_get_cross_normalized (GMT, &GG[4*g1], &GG[4*g2], X); /* X is great circle intersection */
465 gmt_cart_to_geo (GMT, &S->data[GMT_Y][k], &S->data[GMT_X][k], X, true); /* Get lon/lat of this point along crossing profile */
466 S->data[GMT_Y][k] = gmt_lat_swap (GMT, S->data[GMT_Y][k], GMT_LATSWAP_G2O + 1); /* Convert to geodetic */
467 S->data[GMT_Z][k] = S->data[GMT_Z][k+1] = hypot (GG[4*g1+3], GG[4*g2+3]); /* Combined length in quadrature */
468 S->data[3][k] = S->data[3][k+1] = gmt_dot3v (GMT, &GG[4*g1], &GG[4*g2]); /* Cos of angle between great circles is our other weight */
469 S->data[4][k] = S->data[4][k+1] = code[g1] + code[g2]; /* 0 = AH&AH, 1 = AH&FZ, 2 = FZ&FZ */
470 S->data[GMT_X][k+1] = S->data[GMT_X][k] + 180.0;
471 S->data[GMT_Y][k+1] = -S->data[GMT_Y][k];
472 k += 2;
473 }
474 }
475 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->S.file, C) != GMT_NOERROR) {
476 Return (API->error);
477 }
478 gmt_M_free (GMT, GG);
479 gmt_M_free (GMT, code);
480 }
481 }
482 else if (Ctrl->S.mode == SPOTTER_SCAN_LINES) { /* Determine which lines are compatible with the selected test pole */
483 double out[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, plat, sum_L = 0.0, del_angle, chi2, this_chi2;
484 unsigned int n_out = (Ctrl->S.midpoint) ? 6 : 3;
485 GMT_Report (API, GMT_MSG_INFORMATION, "Entering scan mode: lines [EXPERIMENTAL]\n");
486 gmt_set_cartesian (GMT, GMT_OUT); /* Since x here will be table number and y is segment number */
487 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data output */
488 Return (API->error);
489 }
490 if ((error = GMT_Set_Columns (API, GMT_OUT, n_out, GMT_COL_FIX)) != GMT_NOERROR) {
491 Return (error);
492 }
493 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
494 Return (API->error);
495 }
496 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
497
498 plat = gmt_lat_swap (GMT, Ctrl->S.plat, GMT_LATSWAP_G2O); /* Convert latitude to geodetic */
499 gmt_geo_to_cart (GMT, plat, Ctrl->S.plon, X, true); /* Get x/y/z of selected pole X */
500
501 /* Now visit all our segments */
502 for (d = POLESPOTTER_AH; d <= POLESPOTTER_FZ; d++) {
503 if (In[d] == NULL) continue; /* Don't have this data set */
504 weight = (d == POLESPOTTER_AH) ? Ctrl->A.weight : Ctrl->F.weight;
505 Out->text = (char *)label[d];
506 for (tbl = 0; tbl < In[d]->n_tables; tbl++) {
507 for (seg = 0; seg < In[d]->table[tbl]->n_segments; seg++) { /* For each segment in the table */
508 S = In[d]->table[tbl]->segment[seg]; /* Set shortcut to current segment */
509 if (gmt_parse_segment_item (GMT, S->header, "-D", header)) /* Found -D<val> */
510 seg_weight = 1.0 / atof (header);
511 else
512 seg_weight = weight; /* Already 1/sigma, actually */
513 /* Reminder, latitudes in segments are now geocentric latitudes */
514 chi2 = sum_L = 0.0;
515 if (Ctrl->S.midpoint) GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); /* Output the segment header first */
516 gmt_geo_to_cart (GMT, S->data[GMT_Y][0], S->data[GMT_X][0], P1, true); /* Get x/y/z of first point P1 */
517 for (row = 1; row < S->n_rows; row++) {
518 gmt_geo_to_cart (GMT, S->data[GMT_Y][row], S->data[GMT_X][row], P2, true); /* get x/y/z of 2nd point P2 */
519 L = d_acos (gmt_dot3v (GMT, P1, P2)) * RADIAN2KM; /* Length of this segment */
520 del_angle = polespotter_get_angle_between_trends (GMT, P1, P2, d, X);
521 this_chi2 = pow (del_angle * seg_weight, 2.0); /* The chi2 increment from the P1-P2 line */
522 gmt_M_memcpy (P1, P2, 3, double); /* Let old P2 be next P1 */
523 if (Ctrl->S.midpoint) { /* Report for this mid-point */
524 gmt_cart_to_geo (GMT, &mlat, &mlon, M, true); /* Get lon/lat of the mid point */
525 mlat = gmt_lat_swap (GMT, mlat, GMT_LATSWAP_G2O + 1); /* Convert back to geodetic */
526 out[GMT_X] = mlon; out[GMT_Y] = mlat;
527 out[GMT_Z] = del_angle;
528 out[3] = this_chi2;
529 out[4] = (double)tbl;
530 out[5] = (double)seg;
531 GMT_Put_Record (API, GMT_WRITE_DATA, Out); /* Report <mlon mlat del_angle chi2 tbl seg type> for this M */
532 }
533 else { /* Summarize for this line instead */
534 chi2 += L * this_chi2; /* The weighted chi2 sum from this line */
535 sum_L += L; /* Add up total weight sum */
536 }
537 }
538 if (!Ctrl->S.midpoint) { /* Write <tbl seg chi2 type> for this segment */
539 Out->data[GMT_X] = (double)tbl; Out->data[GMT_Y] = (double)seg; Out->data[GMT_Z] = chi2 / sum_L;
540 GMT_Put_Record (API, GMT_WRITE_DATA, Out); /* We are writing these circles to output */
541 }
542 }
543 }
544 }
545 /* Disables further line output */
546 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further line output */
547 Return (API->error);
548 }
549 gmt_M_free (GMT, Out);
550 }
551 else { /* SPOTTER_SCAN_POLES */
552 double *plon = NULL, *plat = NULL, sum_L = 0.0, del_angle, chi2;
553
554 /* Now visit all our segments to convert to geocentric and to get sum of weights once */
555
556 GMT_Report (API, GMT_MSG_INFORMATION, "Entering scan mode: poles\n");
557 for (d = POLESPOTTER_AH; d <= POLESPOTTER_FZ; d++) {
558 if (In[d] == NULL) continue; /* Don't have this data set */
559 weight = (d == POLESPOTTER_AH) ? Ctrl->A.weight : Ctrl->F.weight;
560 for (tbl = 0; tbl < In[d]->n_tables; tbl++) {
561 for (seg = 0; seg < In[d]->table[tbl]->n_segments; seg++) { /* For each segment in the table */
562 S = In[d]->table[tbl]->segment[seg]; /* Set shortcut to current segment */
563 /* Convert the entire segment to geocentric latitude as we go through */
564 S->data[GMT_Y][0] = gmt_lat_swap (GMT, S->data[GMT_Y][0], GMT_LATSWAP_G2O);
565 gmt_geo_to_cart (GMT, S->data[GMT_Y][0], S->data[GMT_X][0], P1, true); /* get x/y/z of first point P1 */
566 for (row = 1; row < S->n_rows; row++) {
567 S->data[GMT_Y][row] = gmt_lat_swap (GMT, S->data[GMT_Y][row], GMT_LATSWAP_G2O);
568 gmt_geo_to_cart (GMT, S->data[GMT_Y][row], S->data[GMT_X][row], P2, true); /* get x/y/z of 2nd point P2 */
569 L = d_acos (gmt_dot3v (GMT, P1, P2)) * RADIAN2KM; /* Length of this segment */
570 sum_L += L; /* Add up total weight sum */
571 gmt_M_memcpy (P1, P2, 3, double); /* Let old P2 be next P1 */
572 }
573 }
574 }
575 }
576 /* Now we know sum_L which we will divide our grid by at the end */
577
578 plon = gmt_grd_coord (GMT, Grid->header, GMT_X);
579 plat = gmt_grd_coord (GMT, Grid->header, GMT_Y);
580 for (grow = 0; grow < Grid->header->n_rows; grow++) { /* Try all possible pole latitudes in selected region */
581 plat[grow] = gmt_lat_swap (GMT, plat[grow], GMT_LATSWAP_G2O); /* Convert latitude to geodetic */
582 for (gcol = 0; gcol < Grid->header->n_columns; gcol++) { /* Try all possible pole longitudes in selected region */
583 node = gmt_M_ijp (Grid->header, grow, gcol); /* Current grid node */
584 gmt_geo_to_cart (GMT, plat[grow], plon[gcol], X, true); /* Get x/y/z of current pole X */
585 /* Now visit all our segments */
586 for (d = POLESPOTTER_AH; d <= POLESPOTTER_FZ; d++) {
587 if (In[d] == NULL) continue; /* Don't have this data set */
588 weight = (d == POLESPOTTER_AH) ? Ctrl->A.weight : Ctrl->F.weight;
589 for (tbl = 0; tbl < In[d]->n_tables; tbl++) {
590 for (seg = 0; seg < In[d]->table[tbl]->n_segments; seg++) { /* For each segment in the table */
591 S = In[d]->table[tbl]->segment[seg]; /* Set shortcut to current segment */
592 if (gmt_parse_segment_item (GMT, S->header, "-D", header)) /* Found -D<val> */
593 seg_weight = 1.0 / atof (header);
594 else
595 seg_weight = weight; /* Already 1/sigma, actually */
596 /* Reminder, latitudes in segments are now geocentric latitudes */
597 gmt_geo_to_cart (GMT, S->data[GMT_Y][0], S->data[GMT_X][0], P1, true); /* Get x/y/z of first point P1 */
598 for (row = 1; row < S->n_rows; row++) {
599 gmt_geo_to_cart (GMT, S->data[GMT_Y][row], S->data[GMT_X][row], P2, true); /* get x/y/z of 2nd point P2 */
600 L = d_acos (gmt_dot3v (GMT, P1, P2)) * RADIAN2KM; /* Length of this segment */
601 del_angle = polespotter_get_angle_between_trends (GMT, P1, P2, d, X);
602 chi2 = L * pow (del_angle * seg_weight, 2.0); /* The weighted chi2 increment from this line */
603 Grid->data[node] += (gmt_grdfloat)chi2; /* Add to total chi2 misfit for this pole */
604 gmt_M_memcpy (P1, P2, 3, double); /* Let old P2 be next P1 */
605 }
606 }
607 }
608 }
609 }
610 }
611 gmt_M_free (GMT, plon);
612 gmt_M_free (GMT, plat);
613
614 for (node = 0; node < Grid->header->size; node++) Grid->data[node] /= (gmt_grdfloat)sum_L; /* Correct for weight sum */
615 }
616 if (Ctrl->G.active) { /* Write the spotting grid */
617 double max = -DBL_MAX;
618 if (Ctrl->N.active) { /* Normalize grid */
619 for (node = 0; node < Grid->header->size; node++) if (Grid->data[node] > max) max = Grid->data[node]; /* Find max value */
620 max = 1.0 / max; /* Do division here */
621 for (node = 0; node < Grid->header->size; node++) Grid->data[node] *= (gmt_grdfloat)max; /* Normalize */
622 }
623 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
624 Return (API->error);
625 }
626 }
627
628 Return (GMT_NOERROR);
629 }
630