1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 2009-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See README file for copying and redistribution conditions.
5 *--------------------------------------------------------------------*/
6 /*
7 * mgd77magref produces output derived from input locations and time and
8 * the CM4 or IGRF magnetic field models.
9 *
10 * Author: Joaquim Luis and Paul Wessel
11 * Date: 1-MAY-2009
12 * Version: 1.0
13 *
14 */
15
16 #include "gmt_dev.h"
17 #include "mgd77.h"
18
19 #define THIS_MODULE_CLASSIC_NAME "mgd77magref"
20 #define THIS_MODULE_MODERN_NAME "mgd77magref"
21 #define THIS_MODULE_LIB "mgd77"
22 #define THIS_MODULE_PURPOSE "Evaluate the IGRF or CM4 magnetic field models"
23 #define THIS_MODULE_KEYS "<D{,>D}"
24 #define THIS_MODULE_NEEDS ""
25 #define THIS_MODULE_OPTIONS "-Vbdho" GMT_OPT("Hm")
26
27 struct MGD77MAGREF_CTRL { /* All control options for this program (except common args) */
28 /* active is true if the option has been activated */
29 struct MGD77_CM4 *CM4;
30 bool copy_input;
31 bool do_IGRF;
32 bool do_CM4;
33 bool joint_IGRF_CM4;
34 struct MGD77MAGREF_A { /* -A */
35 bool active;
36 bool fixed_alt;
37 bool fixed_time;
38 int years;
39 int t_col;
40 double altitude;
41 double time;
42 } A;
43 struct MGD77MAGREF_C { /* -C */
44 bool active;
45 } C;
46 struct MGD77MAGREF_D { /* -D */
47 bool active;
48 } D;
49 struct MGD77MAGREF_F { /* -F */
50 bool active;
51 } F;
52 struct MGD77MAGREF_G { /* -G */
53 bool active;
54 } G;
55 struct MGD77MAGREF_L { /* -L */
56 bool active;
57 } L;
58 struct MGD77MAGREF_S { /* -S */
59 bool active;
60 } S;
61 };
62
New_Ctrl(struct GMT_CTRL * GMT)63 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
64 struct MGD77MAGREF_CTRL *C = NULL;
65
66 C = gmt_M_memory (GMT, NULL, 1, struct MGD77MAGREF_CTRL);
67 C->CM4 = calloc (1U, sizeof (struct MGD77_CM4));
68
69 /* Initialize values whose defaults are not 0/false/NULL */
70
71 C->do_CM4 = true;
72 return (C);
73 }
74
Free_Ctrl(struct GMT_CTRL * GMT,struct MGD77MAGREF_CTRL * C)75 static void Free_Ctrl (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *C) { /* Deallocate control structure */
76 if (!C) return;
77 gmt_M_str_free (C->CM4->CM4_M.path);
78 gmt_M_str_free (C->CM4->CM4_D.path);
79 gmt_M_str_free (C->CM4->CM4_I.path);
80 gmt_M_str_free (C->CM4);
81 gmt_M_free (GMT, C);
82 }
83
usage(struct GMTAPI_CTRL * API,int level)84 static int usage (struct GMTAPI_CTRL *API, int level) {
85 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
86 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
87 GMT_Usage (API, 0, "usage: %s [<table>] [-A+a<alt>+t<date>+y] [-C<cm4file>] [-D<dstfile>] [-E<f107file>] "
88 "[-Frthxyzdi[/[0|9]1234567]] [-G] [-Lrtxyz[/1234]] [-Sc|l<low>/<high>] [%s "
89 "[%s] [%s] [%s] [%s] [%s] [%s]\n",
90 name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_h_OPT, GMT_o_OPT, GMT_colon_OPT, GMT_PAR_OPT);
91
92 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
93
94 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
95 GMT_Usage (API, 1, "\n<table>");
96 GMT_Usage (API, -2, "File with records that must contain <lon>, <lat>, <alt>, <time>[, other cols]. "
97 "Here, (<lon>, <lat>) is the geocentric position on the ellipsoid [but see -G], "
98 "<alt> is the altitude in km positive above the ellipsoid, and "
99 "<time> is the time of data acquisition, in <date>T<clock> format (but see -A+y). "
100 "We read <stdin> if no input file is given.");
101 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
102 GMT_Usage (API, 1, "\n-A+a<alt>+t<date>+y");
103 GMT_Usage (API, -2, "Adjust how the input records are interpreted. Append modifiers:");
104 GMT_Usage (API, 3, "+a Append <alt> to indicate a constant altitude [Default is 3rd column].");
105 GMT_Usage (API, 3, "+t Append <time> to indicate a constant time [Default is 4th column].");
106 GMT_Usage (API, 3, "+y Indicate times are given in decimal years [Default is ISO <date>T<clock> format].");
107 GMT_Usage (API, 1, "\n-C<cm4file>");
108 GMT_Usage (API, -2, "Select an alternate file with coefficients for the CM4 model [%s/umdl.CM4].",
109 API->GMT->session.SHAREDIR);
110 GMT_Usage (API, 1, "\n-D<dstfile>");
111 GMT_Usage (API, -2, "Select an alternate file with hourly means of the Dst index for CM4 [%s/Dst_all.wdc], "
112 "OR a single Dst index to apply for all records.",
113 API->GMT->session.SHAREDIR);
114 GMT_Usage (API, 1, "\n-E<f107file>");
115 GMT_Usage (API, -2, "Select an alternate file with monthly means of absolute F10.7 solar radio flux for CM4 [%s/F107_mon.plt], "
116 "OR a single solar radio flux to apply for all records.",
117 API->GMT->session.SHAREDIR);
118 GMT_Usage (API, 1, "\n-Frthxyzdi[/[0|9]1234567]");
119 GMT_Usage (API, -2, "Dataflags is a string made up of one or more of these codes:");
120 GMT_Usage (API, 3, "r: Output all input columns before adding the items below (all in nTesla).");
121 GMT_Usage (API, 3, "t: List total field.");
122 GMT_Usage (API, 3, "h: List horizontal field.");
123 GMT_Usage (API, 3, "x: List X component.");
124 GMT_Usage (API, 3, "y: List Y component.");
125 GMT_Usage (API, 3, "z: List Z component.");
126 GMT_Usage (API, 3, "d: List declination.");
127 GMT_Usage (API, 3, "i: List inclination.");
128 GMT_Usage (API, -2, "Optionally, append one or more numbers to indicate the requested field contribution(s):");
129 GMT_Usage (API, 3, "0: Core field from IGRF only (no CM4 evaluation).");
130 GMT_Usage (API, 3, "1: Core field.");
131 GMT_Usage (API, 3, "2: Lithospheric field.");
132 GMT_Usage (API, 3, "3: Primary Magnetospheric field.");
133 GMT_Usage (API, 3, "4: Induced Magnetospheric field.");
134 GMT_Usage (API, 3, "5: Primary ionospheric field.");
135 GMT_Usage (API, 3, "6: Induced ionospheric field.");
136 GMT_Usage (API, 3, "7: Toroidal field.");
137 GMT_Usage (API, 3, "9: Core field from IGRF and other contributions from CM4. DO NOT USE BOTH 1 AND 9.");
138 GMT_Usage (API, -2, "Note: Append several numbers to add up the different contributions. For example, "
139 "-Ft/12 computes the total field due to CM4 Core and Lithospheric sources. "
140 "Two special cases are allowed which mix Core field from IGRF and other sources from CM4: "
141 "-Ft/934 computes Core field due to IGRF plus terms 3 and 4 from CM4. "
142 "-Fxyz/934 the same as above but output the field components. "
143 "The data are written out in the order specified "
144 "[Default is -Frthxyzdi/1].");
145 GMT_Usage (API, 1, "\n-G Specify that coordinates are geocentric [geodetic].");
146 GMT_Usage (API, 1, "\n-Lrtxyz[/1234]");
147 GMT_Usage (API, -2, "Compute J field vectors from certain external sources. "
148 "Append a string made up of one or more of these codes:");
149 GMT_Usage (API, 3, "r: Output all input columns before adding the items below (all in Ampere/m).");
150 GMT_Usage (API, 3, "t: List magnitude field.");
151 GMT_Usage (API, 3, "x: List X component.");
152 GMT_Usage (API, 3, "y: List Y component.");
153 GMT_Usage (API, 3, "z: List Z or current function Psi.");
154 GMT_Usage (API, -2, "Optionally, append a number to indicate the requested J contribution(s):");
155 GMT_Usage (API, -2, "1: Induced Magnetospheric field.");
156 GMT_Usage (API, -2, "2: Primary ionospheric field.");
157 GMT_Usage (API, -2, "3: Induced ionospheric field.");
158 GMT_Usage (API, -2, "4: Poloidal field.");
159 GMT_Usage (API, 1, "\n-Sc|l<low>/<high>");
160 GMT_Usage (API, -2, "Limit the CM4 contributions from core and lithosphere to certain harmonic degree bands. "
161 "Append c(ore) or l(ithosphere) and the <low> and <high> degrees to use [-Sc1/13 -Sl14/65].");
162 GMT_Option (API, "V,bi0");
163 if (gmt_M_showusage (API)) {
164 GMT_Usage (API, -2, "Default is 4 input columns (unless -A is used). Note for binary input, absolute time must "
165 "be in the UNIX time-system (unless -A+y is used).");
166 }
167 GMT_Option (API, "bo,d,h,o,:,.");
168
169 return (GMT_MODULE_USAGE);
170 }
171
parse(struct GMT_CTRL * GMT,struct MGD77MAGREF_CTRL * Ctrl,struct GMT_OPTION * options)172 static int parse (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *Ctrl, struct GMT_OPTION *options) {
173 /* This parses the options provided to mgd77magref and sets parameters in CTRL.
174 * Any GMT common options will override values set previously by other commands.
175 * It also replaces any file names specified as input or output with the data ID
176 * returned when registering these sources/destinations with the API.
177 */
178
179 unsigned int n_errors = 0, pos, n_out, lfval = 0, pos_slash = 0, nval = 0, nfval = 0, lval = 0;
180 int j;
181 char p[GMT_BUFSIZ] = {""}, tfixed[GMT_LEN64] = {""};
182 bool do_CM4core = false;
183 struct GMT_OPTION *opt = NULL;
184 struct GMTAPI_CTRL *API = GMT->parent;
185
186 for (opt = options; opt; opt = opt->next) {
187 switch (opt->option) {
188
189 case '<': /* Skip input files */
190 break;
191
192 /* Processes program-specific parameters */
193
194 case 'A':
195 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
196 Ctrl->A.active = true;
197 pos = 0;
198 while ((gmt_strtok (opt->arg, "+", &pos, p))) {
199 switch (p[0]) {
200 case 'a':
201 Ctrl->A.fixed_alt = true;
202 Ctrl->A.altitude = atof (&p[1]);
203 break;
204 case 't':
205 Ctrl->A.fixed_time = true;
206 strncpy (tfixed, &p[1], GMT_LEN64-1);
207 gmt_set_column_type (GMT, GMT_OUT, 3, GMT_IS_FLOAT);
208 break;
209 case 'y':
210 Ctrl->A.years = true;
211 gmt_set_column_type (GMT, GMT_IO, 2, GMT_IS_FLOAT);
212 gmt_set_column_type (GMT, GMT_IO, 3, GMT_IS_FLOAT);
213 break;
214 default:
215 break;
216 }
217 }
218 if (Ctrl->A.fixed_time) {
219 if (Ctrl->A.years)
220 Ctrl->A.time = atof (tfixed);
221 else
222 gmt_scanf_arg (GMT, tfixed, GMT_IS_ABSTIME, false, &Ctrl->A.time);
223 }
224 break;
225 case 'C': /* Alternate CM4 coefficient file */
226 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
227 Ctrl->C.active = true;
228 gmt_M_str_free (Ctrl->CM4->CM4_M.path);
229 Ctrl->CM4->CM4_M.path = strdup (opt->arg);
230 break;
231 case 'D':
232 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
233 Ctrl->D.active = true;
234 j = 0;
235 if (opt->arg[j] == '-') j++;
236 if ((opt->arg[j] > 47) && (opt->arg[j] < 58)) { /* arg is numeric -> Dst Index */
237 Ctrl->CM4->CM4_D.dst[0] = atof (&opt->arg[j]);
238 Ctrl->CM4->CM4_D.index = false;
239 }
240 else {
241 gmt_M_str_free (Ctrl->CM4->CM4_D.path);
242 Ctrl->CM4->CM4_D.path = strdup (opt->arg);
243 Ctrl->CM4->CM4_D.load = true;
244 }
245 break;
246 case 'E':
247 if ((opt->arg[0] > 47) && (opt->arg[0] < 58)) { /* arg is numeric -> Dst Index */
248 Ctrl->CM4->CM4_I.F107 = atof (opt->arg);
249 Ctrl->CM4->CM4_I.index = false;
250 }
251 else {
252 gmt_M_str_free (Ctrl->CM4->CM4_I.path);
253 Ctrl->CM4->CM4_I.path = strdup (opt->arg);
254 Ctrl->CM4->CM4_I.load = true;
255 }
256 break;
257 case 'F':
258 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
259 Ctrl->CM4->CM4_F.active = Ctrl->F.active = true;
260
261 pos_slash = 0;
262 for (j = 0; opt->arg[j]; j++) {
263 if (opt->arg[j] == '/') {
264 pos_slash = j;
265 break;
266 }
267 switch (opt->arg[j]) {
268 case 'r': /* Echo input record */
269 Ctrl->copy_input = true;
270 break;
271 case 't': /* Total field is requested */
272 Ctrl->CM4->CM4_F.field_components[nval++] = 0;
273 break;
274 case 'h': /* Horizontal field is requested */
275 Ctrl->CM4->CM4_F.field_components[nval++] = 1;
276 break;
277 case 'x': /* X component is requested */
278 Ctrl->CM4->CM4_F.field_components[nval++] = 2;
279 break;
280 case 'y': /* Y component is requested */
281 Ctrl->CM4->CM4_F.field_components[nval++] = 3;
282 break;
283 case 'z': /* Z component is requested */
284 Ctrl->CM4->CM4_F.field_components[nval++] = 4;
285 break;
286 case 'd': /* Declination is requested */
287 Ctrl->CM4->CM4_F.field_components[nval++] = 5;
288 break;
289 case 'i': /* Inclination is requested */
290 Ctrl->CM4->CM4_F.field_components[nval++] = 6;
291 break;
292 }
293 }
294 Ctrl->CM4->CM4_F.n_field_components = (int)nval;
295
296 if (pos_slash) {
297 for (j = pos_slash; opt->arg[j]; j++) {
298 switch (opt->arg[j]) {
299 case '0': /* IGRF */
300 Ctrl->do_IGRF = true;
301 Ctrl->do_CM4 = false;
302 Ctrl->joint_IGRF_CM4 = false;
303 break;
304 case '1': /* Main field 1 */
305 Ctrl->CM4->CM4_F.field_sources[nfval++] = 0;
306 do_CM4core = true;
307 break;
308 case '2': /* Main field 2 */
309 Ctrl->CM4->CM4_F.field_sources[nfval++] = 1;
310 break;
311 case '3': /* Primary Magnetospheric field */
312 Ctrl->CM4->CM4_F.field_sources[nfval++] = 2;
313 break;
314 case '4': /* Induced Magnetospheric field */
315 Ctrl->CM4->CM4_F.field_sources[nfval++] = 3;
316 break;
317 case '5': /* Primary ionospheric field */
318 Ctrl->CM4->CM4_F.field_sources[nfval++] = 4;
319 break;
320 case '6': /* Induced ionospheric field */
321 Ctrl->CM4->CM4_F.field_sources[nfval++] = 5;
322 break;
323 case '7': /* Toroidal field */
324 Ctrl->CM4->CM4_F.field_sources[nfval++] = 6;
325 break;
326 case '9': /* Main field is computed with the IGRF */
327 Ctrl->do_IGRF = false;/* No contradiction, the way will be through joint_IGRF_CM4 */
328 Ctrl->do_CM4 = true; /* Well maybe, if some other source is selected also */
329 Ctrl->joint_IGRF_CM4 = true;
330 break;
331 }
332 }
333 Ctrl->CM4->CM4_F.n_field_sources = (int)nfval;
334 }
335 break;
336 case 'G':
337 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
338 Ctrl->G.active = true;
339 Ctrl->CM4->CM4_G.geodetic = false;
340 break;
341 case 'L':
342 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
343 Ctrl->CM4->CM4_L.curr = Ctrl->L.active = true;
344
345 pos_slash = 0;
346 for (j = 0; opt->arg[j]; j++) {
347 if (opt->arg[j] == '/') {
348 pos_slash = j;
349 break;
350 }
351 switch (opt->arg[j]) {
352 case 'r': /* Echo input record */
353 Ctrl->copy_input = true;
354 break;
355 case 't': /* Total field is requested */
356 Ctrl->CM4->CM4_L.curr_components[lval++] = 0;
357 break;
358 case 'x': /* X component is requested */
359 Ctrl->CM4->CM4_L.curr_components[lval++] = 1;
360 break;
361 case 'y': /* Y component is requested */
362 Ctrl->CM4->CM4_L.curr_components[lval++] = 2;
363 break;
364 case 'z': /* Z component is requested */
365 Ctrl->CM4->CM4_L.curr_components[lval++] = 3;
366 break;
367 }
368 }
369 Ctrl->CM4->CM4_L.n_curr_components = (int)lval;
370
371 if (pos_slash) {
372 for (j = pos_slash; opt->arg[j]; j++) {
373 switch (opt->arg[j]) {
374 case '1': /* Induced Magnetospheric field */
375 Ctrl->CM4->CM4_L.curr_sources[lfval++] = 0;
376 break;
377 case '2': /* Primary ionospheric field */
378 Ctrl->CM4->CM4_L.curr_sources[lfval++] = 1;
379 break;
380 case '3': /* Induced ionospheric field */
381 Ctrl->CM4->CM4_L.curr_sources[lfval++] = 2;
382 break;
383 case '4': /* Poloidal field */
384 Ctrl->CM4->CM4_L.curr_sources[lfval++] = 3;
385 break;
386 }
387 }
388 Ctrl->CM4->CM4_L.n_curr_sources = (int)lfval;
389 }
390 break;
391 case 'S':
392 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
393 Ctrl->S.active = true;
394 if (opt->arg[0] == 'c') {
395 j = sscanf (&opt->arg[1], "%d/%d", &Ctrl->CM4->CM4_S.nlmf[0], &Ctrl->CM4->CM4_S.nhmf[0]);
396 if (j != 2) {
397 GMT_Report (API, GMT_MSG_ERROR, "-Sc option usage is -Sc<low/high>\n");
398 n_errors++;
399 }
400 }
401 if (opt->arg[0] == 'l') {
402 j = sscanf (&opt->arg[1], "%d/%d", &Ctrl->CM4->CM4_S.nlmf[1], &Ctrl->CM4->CM4_S.nhmf[1]);
403 if (j != 2) {
404 GMT_Report (API, GMT_MSG_ERROR, "-Sl option usage is -Sl<low/high>\n");
405 n_errors++;
406 }
407 }
408 break;
409 default: /* Report bad options */
410 n_errors += gmt_default_error (GMT, opt->option);
411 break;
412 }
413 }
414
415 n_out = 4 - (Ctrl->A.fixed_alt + Ctrl->A.fixed_time); /* Minimum input columns (could be more) */
416 if (GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0)
417 GMT->common.b.ncol[GMT_IN] = n_out;
418 n_errors += gmt_M_check_condition (GMT, GMT->common.b.active[GMT_IN] && GMT->common.b.ncol[GMT_IN] == 0,
419 "Binary input data (-bi) must have at least %d columns\n", n_out);
420 n_errors += gmt_M_check_condition (GMT, Ctrl->CM4->CM4_F.active && Ctrl->CM4->CM4_L.curr,
421 "You cannot select both -F and -L options\n");
422 n_errors += gmt_M_check_condition (GMT, (do_CM4core && Ctrl->do_IGRF) || (do_CM4core && Ctrl->joint_IGRF_CM4),
423 "You cannot select both CM4 core (1) and IGRF as they are both core fields.\n");
424
425 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
426 }
427
428 #define bailout(code) {gmt_M_free_options (mode); return (code);}
429 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
430
GMT_mgd77magref(void * V_API,int mode,void * args)431 EXTERN_MSC int GMT_mgd77magref (void *V_API, int mode, void *args) {
432 unsigned int j, nval = 0, nfval = 0, error = 0;
433 unsigned int lval = 0, lfval = 0, n_field_components, tbl;
434 unsigned int n_out = 0, n_in, t_col = 3;
435 bool cm4_igrf_T = false;
436
437 size_t i, s, need = 0, n_alloc = 0;
438
439 double the_altitude, the_time, *time_array = NULL, *alt_array = NULL, *time_years = NULL, IGRF[7], out[GMT_MAX_COLUMNS];
440 double *igrf_xyz = NULL; /* Temporary storage for the joint_IGRF_CM4 case */
441
442 struct MGD77_CONTROL M;
443 struct MGD77MAGREF_CTRL *Ctrl = NULL;
444 struct GMT_DATASET *Din = NULL;
445 struct GMT_DATATABLE *T = NULL;
446 struct GMT_RECORD *Out = NULL;
447 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
448 struct GMT_OPTION *options = NULL;
449 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
450
451 /*----------------------- Standard module initialization and parsing ----------------------*/
452
453 if (API == NULL) return (GMT_NOT_A_SESSION);
454 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
455 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
456
457 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
458
459 /* Parse the command-line arguments */
460
461 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 */
462 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
463 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
464 MGD77_Init (GMT, &M); /* Initialize MGD77 Machinery */
465 MGD77_CM4_init (GMT, &M, Ctrl->CM4); /* Presets path using strdup */
466 if ((error = parse (GMT, Ctrl, options)) != 0) {
467 MGD77_end (GMT, &M);
468 Return (error);
469 }
470
471 /*---------------------------- This is the mgd77magref main code ----------------------------*/
472
473 n_in = 4 - (Ctrl->A.fixed_alt + Ctrl->A.fixed_time);
474
475 if (Ctrl->A.fixed_alt) t_col = 2; /* Since we are missing the altitude column */
476
477 Ctrl->CM4->CM4_D.dst = calloc (1U, sizeof(double)); /* We need at least a size of one in case a value is given in input */
478 if (!Ctrl->A.fixed_time) /* Otherwise we don't print the time */
479 gmt_set_column_type (GMT, GMT_IO, t_col, GMT_IS_ABSTIME);
480
481 /* Shorthand for these */
482 nval = Ctrl->CM4->CM4_F.n_field_components;
483 nfval = Ctrl->CM4->CM4_F.n_field_sources;
484 lval = Ctrl->CM4->CM4_L.n_curr_components;
485 lfval = Ctrl->CM4->CM4_L.n_curr_sources;
486
487 if (nfval && Ctrl->do_IGRF) GMT_Message (API, GMT_TIME_NONE, "Warning. Source fields other than IGRF will be ignored. It's in the manual\n");
488
489 /* ------------- Test, and take measures, if mix mode IGRF/CM4 is to be used -------------------------- */
490 if (Ctrl->joint_IGRF_CM4) {
491 if (nfval == 0) { /* It means we had a -F.../9 which is exactly iqual to -F.../0 */
492 Ctrl->do_IGRF = true;
493 Ctrl->do_CM4 = false;
494 Ctrl->joint_IGRF_CM4 = false;
495 }
496 else {
497 cm4_igrf_T = false;
498 if ( (nval == 1) && (Ctrl->CM4->CM4_F.field_components[0] == 0) ) {
499 for (i = 0; i < 3; i++) Ctrl->CM4->CM4_F.field_components[i] = (int)i+2; /* Force the x,y,z request */
500 cm4_igrf_T = true;
501 }
502 else if (!((nval == 3) && (Ctrl->CM4->CM4_F.field_components[0] == 2) && (Ctrl->CM4->CM4_F.field_components[1] == 3) &&
503 (Ctrl->CM4->CM4_F.field_components[2] == 4)) ) {
504 GMT_Report (API, GMT_MSG_ERROR, "In mix CM4/IGRF mode -F option can only be -Ft[r]/... or -Fxyz[r]/...\n");
505 error++;
506 }
507
508 nval = 3;
509 Ctrl->CM4->CM4_F.n_field_components = (int)nval;
510 }
511 }
512 /* ----------------------------------------------------------------------------------------------------- */
513
514 if (error) Return (GMT_RUNTIME_ERROR);
515
516 if (!Ctrl->CM4->CM4_F.active && !Ctrl->CM4->CM4_L.curr) Ctrl->CM4->CM4_F.active = true;
517
518 /* Sort the order in which the parameters appear */
519 if (Ctrl->CM4->CM4_F.active) {
520 if (nval == 0) { /* No components selected, default used */
521 Ctrl->copy_input = true;
522 Ctrl->CM4->CM4_F.n_field_components = 7;
523 for (i = 0; i < 7; i++) Ctrl->CM4->CM4_F.field_components[i] = (int)i;
524 }
525 if (nfval == 0 && !Ctrl->joint_IGRF_CM4) { /* No sources selected, retain only the main field */
526 Ctrl->CM4->CM4_F.field_sources[0] = 0;
527 Ctrl->CM4->CM4_F.n_field_sources = 1;
528 }
529 n_field_components = Ctrl->CM4->CM4_F.n_field_components;
530 }
531 else {
532 if (lval == 0) { /* Nothing selected, default used */
533 Ctrl->copy_input = true;
534 Ctrl->CM4->CM4_L.n_curr_components = 1;
535 Ctrl->CM4->CM4_L.curr_components[0] = 0;
536 }
537 if (lfval == 0) { /* Nothing selected, retain only the induced magnetospheric field */
538 Ctrl->CM4->CM4_L.curr_sources[0] = 0;
539 Ctrl->CM4->CM4_L.n_curr_sources = 1;
540 }
541 n_field_components = Ctrl->CM4->CM4_L.n_curr_components;
542 }
543
544 if (Ctrl->A.fixed_alt) { /* A single altitude should apply to all records; set array to point to this single altitude */
545 alt_array = &Ctrl->A.altitude;
546 Ctrl->CM4->CM4_DATA.n_altitudes = 1;
547 }
548 if (Ctrl->A.fixed_time) { /* A single time should apply to all records; set array to point to this single time */
549 if (!Ctrl->A.years) {
550 if (M.adjust_time) Ctrl->A.time = MGD77_time2utime (GMT, &M, Ctrl->A.time); /* Convert to Unix time if need be */
551 Ctrl->A.time = MGD77_time_to_fyear (GMT, &M, Ctrl->A.time); /* Get decimal year */
552 }
553 time_array = &Ctrl->A.time;
554 Ctrl->CM4->CM4_DATA.n_times = 1;
555 }
556 else /* Make sure input time columns are encoded/decoded properly since here we know t_col is set. */
557 gmt_set_column_type (GMT, GMT_IO, t_col, (Ctrl->A.years) ? GMT_IS_FLOAT : GMT_IS_ABSTIME);
558
559 gmt_set_column_type (GMT, GMT_IO, t_col+1, GMT_IS_FLOAT); /* Override any previous t_col = 3 settings */
560 if (!Ctrl->copy_input) { /* No time on output */
561 gmt_set_column_type (GMT, GMT_OUT, 2, GMT_IS_FLOAT);
562 gmt_set_column_type (GMT, GMT_OUT, 3, GMT_IS_FLOAT);
563 }
564
565 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_PLP, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */
566 Return (API->error);
567 }
568 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_PLP, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
569 Return (API->error);
570 }
571
572 if ((Din = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) {
573 Return (API->error);
574 }
575 if (Din->n_columns < n_in) {
576 GMT_Report (API, GMT_MSG_ERROR, "Input data have %d column(s) but at least %d are needed\n", (int)Din->n_columns, n_in);
577 Return (GMT_DIM_TOO_SMALL);
578 }
579 n_out = n_field_components + ((Ctrl->copy_input) ? (unsigned int)Din->n_columns : 0);
580 if (cm4_igrf_T) n_out -= 2; /* Decrease by 2 because the x,y,z were imposed internally only. i.e not for output */
581 if ((error = GMT_Set_Columns (API, GMT_OUT, n_out, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
582 Return (error);
583 }
584
585 if (GMT->common.b.active[GMT_OUT] && GMT->common.b.ncol[GMT_OUT] > 0 && n_out > GMT->common.b.ncol[GMT_OUT]) {
586 GMT_Report (API, GMT_MSG_ERROR, "Binary output must have at least %d columns (your -bo option only set %d)\n", n_out, GMT->common.b.ncol[GMT_OUT]);
587 Return (GMT_RUNTIME_ERROR);
588 }
589
590 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
591 Return (API->error);
592 }
593 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
594 Return (API->error);
595 }
596 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
597
598 for (tbl = 0; tbl < Din->n_tables; tbl++) { /* Loop over all input tables */
599 T = Din->table[tbl]; /* Current table */
600
601 if (T->n_columns < n_in) {
602 GMT_Report (API, GMT_MSG_ERROR, "Table %d has %d columns, but from the used options we expect %d\n",
603 tbl + 1, T->n_columns, n_in);
604 continue;
605 }
606
607 for (s = 0; s < T->n_segments; s++) { /* Process each file segment separately */
608 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, T->segment[s]->header);
609 need = T->segment[s]->n_rows; /* Size of output array needed in MGD77_cm4field */
610 if (need > n_alloc) { /* Need to reallocate */
611 n_alloc = need;
612 Ctrl->CM4->CM4_DATA.out_field = gmt_M_memory (GMT, Ctrl->CM4->CM4_DATA.out_field, n_alloc * n_field_components, double);
613 if (!(Ctrl->A.years || Ctrl->A.fixed_time))
614 time_years = gmt_M_memory (GMT, time_years, n_alloc, double);
615
616 if (Ctrl->joint_IGRF_CM4)
617 igrf_xyz = gmt_M_memory (GMT, igrf_xyz, n_alloc * 3, double);
618 }
619
620 if (!Ctrl->A.fixed_alt) { /* Assign the alt_array to the provided altitude array */
621 alt_array = T->segment[s]->data[GMT_Z];
622 Ctrl->CM4->CM4_DATA.n_altitudes = (int)T->segment[s]->n_rows;
623 }
624
625 if (!Ctrl->A.fixed_time) { /* Assign the time_array to the provided time array */
626 Ctrl->CM4->CM4_DATA.n_times = (int)T->segment[s]->n_rows;
627 if (Ctrl->A.years)
628 time_array = T->segment[s]->data[t_col];
629 else { /* Must convert internal GMT time to decimal years first */
630 for (i = 0; i < T->segment[s]->n_rows; i++)
631 time_years[i] = MGD77_time_to_fyear (GMT, &M, T->segment[s]->data[t_col][i]);
632 time_array = time_years;
633 }
634 }
635
636 if (!(Ctrl->do_IGRF || Ctrl->joint_IGRF_CM4 ) && !s && time_array[0] > 2002.7) { /* Only atmospheric terms may be reliable */
637 GMT_Report (API, GMT_MSG_WARNING, "Time is outside the CM4 strict validity domain [1960.0-2002.7].\n");
638 GMT_Report (API, GMT_MSG_WARNING, "\tThe secular variation estimation will be unreliable. In this"
639 "\n\tcase you really should use the IGRF to estimate the core contribution\n");
640 }
641
642 Ctrl->CM4->CM4_DATA.n_pts = (int)T->segment[s]->n_rows;
643 if (Ctrl->do_IGRF || Ctrl->joint_IGRF_CM4) {
644 int type;
645 type = (Ctrl->CM4->CM4_G.geodetic) ? 1 : 2;
646 for (i = 0; i < T->segment[s]->n_rows; i++) {
647 the_altitude = (Ctrl->A.fixed_alt) ? alt_array[0] : alt_array[i];
648 the_time = (Ctrl->A.fixed_time) ? time_array[0] : time_array[i];
649 if (type == 2) the_altitude += 6371.2;
650 MGD77_igrf10syn (GMT, 0, the_time, type, the_altitude, T->segment[s]->data[GMT_X][i],
651 T->segment[s]->data[GMT_Y][i], IGRF);
652 if (!Ctrl->joint_IGRF_CM4) { /* IGRF only */
653 int jj;
654 for (jj = 0; jj < Ctrl->CM4->CM4_F.n_field_components; jj++)
655 Ctrl->CM4->CM4_DATA.out_field[i*n_field_components+jj] = IGRF[Ctrl->CM4->CM4_F.field_components[jj]];
656 }
657 else { /* Store the IGRF x,y,z components for later use */
658 for (j = 0; j < 3; j++)
659 igrf_xyz[i*3+j] = IGRF[Ctrl->CM4->CM4_F.field_components[j]];
660 }
661 }
662 }
663
664 if (Ctrl->do_CM4) { /* DO CM4 only. Eval CM4 at all points */
665 int err;
666 if ((err = MGD77_cm4field (GMT, Ctrl->CM4, T->segment[s]->data[GMT_X],
667 T->segment[s]->data[GMT_Y], alt_array, time_array)) != 0) {
668 GMT_Report (API, GMT_MSG_ERROR, "this segment has a record generating an error.\n"
669 "Unfortunately, this means all other eventually good\n"
670 "records are also ignored. Fix the bad record and rerun the command.\n");
671 continue;
672 }
673 }
674
675 if ((Ctrl->do_CM4 || Ctrl->do_IGRF) && !Ctrl->joint_IGRF_CM4) { /* DID CM4 or (exclusive) IGRF only. */
676 for (i = 0; i < T->segment[s]->n_rows; i++) { /* Output the requested columns */
677 n_out = 0;
678 if (Ctrl->copy_input)
679 for (j = 0; j < T->segment[s]->n_columns; j++)
680 out[n_out++] = T->segment[s]->data[j][i];
681 for (j = 0; j < n_field_components; j++)
682 out[n_out++] = Ctrl->CM4->CM4_DATA.out_field[i*n_field_components+j];
683
684 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
685 }
686 }
687 else { /* DID CM4 and IGRF */
688 double x, y, z;
689 for (i = 0; i < T->segment[s]->n_rows; i++) { /* Output the requested columns */
690 n_out = 0;
691 if (Ctrl->copy_input) for (j = 0; j < T->segment[s]->n_columns; j++) out[n_out++] = T->segment[s]->data[j][i];
692 if (cm4_igrf_T) {
693 x = Ctrl->CM4->CM4_DATA.out_field[i*3 ] + igrf_xyz[i*3 ];
694 y = Ctrl->CM4->CM4_DATA.out_field[i*3+1] + igrf_xyz[i*3+1];
695 z = Ctrl->CM4->CM4_DATA.out_field[i*3+2] + igrf_xyz[i*3+2];
696 out[n_out++] = sqrt(x*x + y*y + z*z);
697 }
698 else {
699 for (j = 0; j < 3; j++)
700 out[n_out++] = Ctrl->CM4->CM4_DATA.out_field[i*3+j] + igrf_xyz[i*3+j];
701 }
702
703 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
704 }
705 }
706
707 }
708 }
709 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data input */
710 Return (API->error);
711 }
712
713 gmt_M_free (GMT, Out);
714 gmt_M_str_free (Ctrl->CM4->CM4_D.dst);
715 gmt_M_free (GMT, Ctrl->CM4->CM4_DATA.out_field);
716 if (!(Ctrl->A.years || Ctrl->A.fixed_time)) gmt_M_free (GMT, time_years);
717 if (Ctrl->joint_IGRF_CM4) gmt_M_free (GMT, igrf_xyz);
718
719 MGD77_end (GMT, &M);
720
721 Return (GMT_NOERROR);
722 }
723