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