Add option to the IGRF calculations to limit the number of used harmonics.

doc/rst/source/changes.rst

@@ -6,6 +6,14 @@
 Changelog
 =========
 
+New Features in GMT 6.8.0
+=========================
+
+* :doc:`/supplements/mgd77/mgd77magref`: Add the **-F**...\ /**0+l**\ *low/high* modifier that limits the IGRF
+  evaluation to a band of spherical harmonic degrees, in the same spirit as **-S** does for the CM4 core and
+  lithospheric fields.  It is only accepted for the IGRF model (**-F**...\ /**0**) and is not available in
+  GMT 6.7.1 or earlier.
+
 New Features in GMT 6.7.1
 =========================
 

doc/rst/source/supplements/mgd77/mgd77magref.rst

@@ -17,7 +17,7 @@ Synopsis
 [ |-C|\ *cm4file* ]
 [ |-D|\ *Dstfile* ]
 [ |-E|\ *f107file* ]
-[ |-F|\ *flags* ]
+[ |-F|\ *flags*\ [**+l**\ *low/high*] ]
 [ |-G| ]
 [ |-S|\ **c**\|\ **l**\ *low/high* ]
 [ |SYN_OPT-V| ]
@@ -116,7 +116,7 @@ Optional Arguments
 
 .. _-F:
 
-**-F**\ *flags*
+**-F**\ *flags*\ [**+l**\ *low/high*]
     Selects output items; *flags* is a string made up of one or more of
     these characters:
 
@@ -167,6 +167,20 @@ Optional Arguments
     order they appear in *flags* [Default is **-Frthxyzdi**/**1**]. **Note**:
     You can further select a subset of fields via **-o**.
 
+    Append **+l**\ *low/high* to limit the IGRF contribution to the spherical
+    harmonic degrees in the *low/high* band [all degrees]. This modifier is only
+    available for the IGRF model, that is, when **-F**...\ /**0** is used.
+    **Note**: The IGRF only has degrees up to 10 before 1995 and up to 13 from
+    1995 onwards, so a *high* beyond that limit is truncated.
+
+    **Note**: The **+l** modifier requires a version later than GMT 6.7.0. Earlier versions do
+    not know it and will read its digits as CM4 field source codes instead: most
+    bands then fail with a misleading message about selecting both the CM4 core
+    and the IGRF, but a band such as **+l5/2**, whose digits happen to be valid
+    source codes other than the core, is silently ignored and the complete IGRF
+    field is returned. Check with **gmt --version** before using **+l** in a
+    script that may run against an older installation.
+
 .. _-G:
 
 **-G**

src/mgd77/mgd77.c

@@ -4965,7 +4965,7 @@ double MGD77_carter_correction (struct GMT_CTRL *GMT, double lon, double lat, do
   *
 *--------------------------------------------------------------------*/
 
-int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, double alt, double elong, double lat, double *out) {
+int MGD77_igrf10syn_band (struct GMT_CTRL *GMT, int isv, double date, int itype, double alt, double elong, double lat, int nlow, int nhigh, double *out) {
  /*     This is a synthesis routine for the 10th generation IGRF as agreed
   *     in December 2004 by IAGA Working Group V-MOD. It is valid 1900.0 to
   *     2010.0 inclusive. Values for dates from 1945.0 to 2000.0 inclusive are
@@ -4982,6 +4982,9 @@ int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, doub
   *           = distance from centre of Earth in km if itype = 2 (>3485 km)
   *     lat   = latitude (90-90)
   *     elong = east-longitude (0-360) -- it works also in [-180;+180]
+  *     nlow  = lowest spherical harmonic degree to include in the synthesis [1]
+  *     nhigh = highest spherical harmonic degree to include in the synthesis
+  *             Pass nlow = 1 and nhigh = 0 (or 13) to get the complete field.
   *   OUTPUT
   *     out[0] F  = total intensity (nT) if isv = 0, rubbish if isv = 1
   *     out[1] H  = horizontal intensity (nT)
@@ -5534,6 +5537,7 @@ int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, doub
 	 };
 
 	int i, j, k, l, m, n, ll, lm, kmx, nmx, nc;
+	bool in_band;
 	double cd, cl[13], tc, ct, sd, fn = 0.0, gn = 0.0, fm, sl[13];
 	double rr, st, one, gmm, rho, two, three, ratio;
 	double p[105], q[105], r, t, a2, b2;
@@ -5579,6 +5583,18 @@ int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, doub
 		nc = nmx * (nmx + 2);
 		kmx = (nmx + 1) * (nmx + 2) / 2;
 	}
+	/* Check the requested harmonic degree band against what this particular model actually has */
+	if (nlow < 1) nlow = 1;
+	if (nhigh < 1) nhigh = nmx;	/* Not set, so use all the degrees available */
+	if (nhigh > nmx) {
+		GMT_Report(GMT->parent, GMT_MSG_WARNING, "The IGRF model for %g only has harmonic degrees up to %d, so resetting the upper degree from %d to %d\n", date, nmx, nhigh, nmx);
+		nhigh = nmx;
+	}
+	if (nlow > nhigh) {
+		GMT_Report(GMT->parent, GMT_MSG_ERROR, "The lower harmonic degree (%d) exceeds the upper degree (%d)\n", nlow, nhigh);
+		return MGD77_BAD_IGRFDATE;
+	}
+
 	r = alt;
 	sincosd (90.0 - lat, &st, &ct);
 	sincosd (elong, &(sl[0]), &(cl[0]));
@@ -5644,20 +5660,25 @@ int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, doub
 
 		lm = ll + l;
 		one = (tc * gh[lm-1] + t * gh[lm+nc-1]) * rr;
+		in_band = (n >= nlow && n <= nhigh);	/* Only degrees inside the band contribute, but p, q and l must be advanced regardless */
 		if (m == 0) {
-			X += one * q[k-1];
-			Z -= (fn + 1.) * one * p[k-1];
+			if (in_band) {
+				X += one * q[k-1];
+				Z -= (fn + 1.) * one * p[k-1];
+			}
 			l++;
 		}
 		else {
 			two = (tc * gh[lm] + t * gh[lm+nc]) * rr;
 			three = one * cl[m-1] + two * sl[m - 1];
-			X += three * q[k-1];
-			Z -= (fn + 1.) * three * p[k-1];
-			if (st != 0.)
-				Y += (one * sl[m-1] - two * cl[m-1]) * fm * p[k-1] / st;
-			else
-				Y += (one * sl[m-1] - two * cl[m-1]) * q[k-1] * ct;
+			if (in_band) {
+				X += three * q[k-1];
+				Z -= (fn + 1.) * three * p[k-1];
+				if (st != 0.)
+					Y += (one * sl[m-1] - two * cl[m-1]) * fm * p[k-1] / st;
+				else
+					Y += (one * sl[m-1] - two * cl[m-1]) * q[k-1] * ct;
+			}
 			l += 2;
 		}
 		m++;
@@ -5678,6 +5699,11 @@ int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, doub
 	return (MGD77_NO_ERROR);
 }
 
+int MGD77_igrf10syn(struct GMT_CTRL *GMT, int isv, double date, int itype, double alt, double elong, double lat, double *out) {
+	/* Evaluate the complete IGRF field, i.e., using all the harmonic degrees available */
+	return (MGD77_igrf10syn_band (GMT, isv, date, itype, alt, elong, lat, 1, 0, out));
+}
+
 void MGD77_IGF_text (struct GMTAPI_CTRL *API, int indent, int version) {
 	switch (version) {
 		case 1:	/* Heiskanen 1924 model */

src/mgd77/mgd77.h

@@ -565,6 +565,7 @@ EXTERN_MSC double MGD77_carter_correction (struct GMT_CTRL *GMT, double lon, dou
 /* User functions for direct use of IGRF corrections, theoretical gravity */
 
 EXTERN_MSC int MGD77_igrf10syn (struct GMT_CTRL *GMT, int isv, double date, int itype, double alt, double lon, double lat, double *out);
+EXTERN_MSC int MGD77_igrf10syn_band(struct GMT_CTRL *GMT, int isv, double date, int itype, double alt, double lon, double lat, int nlow, int nhigh, double *out);
 EXTERN_MSC double MGD77_Theoretical_Gravity (struct GMT_CTRL *GMT, double lon, double lat, int version);
 EXTERN_MSC void MGD77_IGF_text (struct GMTAPI_CTRL *API, int indent, int version);
 EXTERN_MSC double MGD77_Recalc_Mag_Anomaly_IGRF (struct GMT_CTRL *GMT, struct MGD77_CONTROL *F, double time, double lon, double lat, double obs, bool calc_date);

src/mgd77/mgd77magref.c

@@ -51,6 +51,8 @@ struct MGD77MAGREF_CTRL {	/* All control options for this program (except common
 	} D;
 	struct MGD77MAGREF_F {	/* -F */
 		bool active;
+		bool band;	/* True if +l<low>/<high> was given */
+		int nlow, nhigh;	/* The IGRF harmonic degree band to use */
 	} F;
 	struct MGD77MAGREF_G {	/* -G */
 		bool active;
@@ -72,6 +74,7 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) {	/* Allocate and initialize a new
 	/* Initialize values whose defaults are not 0/false/NULL */
 
 	C->do_CM4 = true;
+	C->F.nlow = 1;	C->F.nhigh = 0;	/* Zero means use all the degrees the IGRF model has */
 	return (C);
 }
 
@@ -88,7 +91,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
 	GMT_Usage (API, 0, "usage: %s [<table>] [-A+a<alt>+t<date>+y] [-C<cm4file>] [-D<dstfile>] [-E<f107file>] "
-		"[-Frthxyzdi[/[0|9]1234567]] [-G] [-Lrtxyz[/1234]] [-Sc|l<low>/<high>] [%s] "
+		"[-Frthxyzdi[/[0|9]1234567][+l<low>/<high>]] [-G] [-Lrtxyz[/1234]] [-Sc|l<low>/<high>] [%s] "
 		"[%s] [%s] [%s] [%s] [%s] [%s]\n",
 		name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_h_OPT, GMT_o_OPT, GMT_colon_OPT, GMT_PAR_OPT);
 
@@ -118,7 +121,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, -2, "Select an alternate file with monthly means of absolute F10.7 solar radio flux for CM4 [%s/F107_mon.plt], "
 		"OR a single solar radio flux to apply for all records.",
 		API->GMT->session.SHAREDIR);
-	GMT_Usage (API, 1, "\n-Frthxyzdi[/[0|9]1234567]");
+	GMT_Usage (API, 1, "\n-Frthxyzdi[/[0|9]1234567][+l<low>/<high>]");
 	GMT_Usage (API, -2, "Dataflags is a string made up of one or more of these codes:");
 	GMT_Usage (API, 3, "r: Output all input columns before adding the items below (all in nTesla).");
 	GMT_Usage (API, 3, "t: List total field.");
@@ -145,6 +148,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 		"-Fxyz/934 the same as above but output the field components. "
 		"The data are written out in the order specified "
 		"[Default is -Frthxyzdi/1].");
+	GMT_Usage(API, 3, "+l Limit the IGRF to the harmonic degrees in the <low>/<high> band [all degrees]. "
+		"Only available for the IGRF model, i.e., with -F.../0. Note: The IGRF only has degrees "
+		"up to 10 before 1995 and up to 13 from 1995 onwards.");
 	GMT_Usage (API, 1, "\n-G Specify that coordinates are geocentric [geodetic].");
 	GMT_Usage (API, 1, "\n-Lrtxyz[/1234]");
 	GMT_Usage (API, -2, "Compute J field vectors from certain external sources. "
@@ -181,7 +187,7 @@ static int parse (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *Ctrl, struct GM
 
 	unsigned int n_errors = 0, pos, n_out, lfval = 0, pos_slash = 0, nval = 0, nfval = 0, lval = 0;
 	int j;
-	char p[GMT_BUFSIZ] = {""}, tfixed[GMT_LEN64] = {""};
+	char p[GMT_BUFSIZ] = {""}, tfixed[GMT_LEN64] = {""}, *c = NULL;
 	bool do_CM4core = false;
 	struct GMT_OPTION *opt = NULL;
 	struct GMTAPI_CTRL *API = GMT->parent;
@@ -258,6 +264,17 @@ static int parse (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *Ctrl, struct GM
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
 				Ctrl->CM4->CM4_F.active = true;
 
+				if ((c = strstr(opt->arg, "+l"))) {	/* Want to limit the IGRF to a harmonic degree band */
+					if (sscanf(&c[2], "%d/%d", &Ctrl->F.nlow, &Ctrl->F.nhigh) != 2) {
+						GMT_Report(API, GMT_MSG_ERROR, "Option -F: The +l modifier usage is +l<low>/<high>\n");
+						Ctrl->F.nlow = 1;	Ctrl->F.nhigh = 0;
+						n_errors++;
+					}
+					else
+						Ctrl->F.band = true;
+					c[0] = '\0';	/* Chop off the modifier while we parse the rest */
+				}
+
 				pos_slash = 0;
 				for (j = 0; opt->arg[j]; j++) {
 					if (opt->arg[j] == '/') {
@@ -332,6 +349,7 @@ static int parse (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *Ctrl, struct GM
 					}
 					Ctrl->CM4->CM4_F.n_field_sources = (int)nfval;
 				}
+				if (c) c[0] = '+';	/* Restore the modifier */
 				break;
 			case 'G':
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
@@ -419,6 +437,10 @@ static int parse (struct GMT_CTRL *GMT, struct MGD77MAGREF_CTRL *Ctrl, struct GM
 			"You cannot select both -F and -L options\n");
 	n_errors += gmt_M_check_condition (GMT, (do_CM4core && Ctrl->do_IGRF) || (do_CM4core && Ctrl->joint_IGRF_CM4),
 			"You cannot select both CM4 core (1) and IGRF as they are both core fields.\n");
+	n_errors += gmt_M_check_condition(GMT, Ctrl->F.band && !Ctrl->do_IGRF,
+			"Option -F: The +l modifier is only available for the IGRF model, i.e., with -F.../0\n");
+	n_errors += gmt_M_check_condition(GMT, Ctrl->F.band && (Ctrl->F.nlow < 1 || Ctrl->F.nhigh < Ctrl->F.nlow),
+			"Option -F: The +l modifier requires 1 <= <low> <= <high>\n");
 
 	return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
 }
@@ -645,8 +667,8 @@ EXTERN_MSC int GMT_mgd77magref (void *V_API, int mode, void *args) {
 					the_altitude = (Ctrl->A.fixed_alt) ? alt_array[0] : alt_array[i];
 					the_time = (Ctrl->A.fixed_time) ? time_array[0] : time_array[i];
 					if (type == 2) the_altitude += 6371.2;
-					MGD77_igrf10syn (GMT, 0, the_time, type, the_altitude, T->segment[s]->data[GMT_X][i],
-							T->segment[s]->data[GMT_Y][i], IGRF);
+					MGD77_igrf10syn_band(GMT, 0, the_time, type, the_altitude, T->segment[s]->data[GMT_X][i],
+							T->segment[s]->data[GMT_Y][i], Ctrl->F.nlow, Ctrl->F.nhigh, IGRF);
 					if (!Ctrl->joint_IGRF_CM4) {		/* IGRF only */
 						int jj;
 						for (jj = 0; jj < Ctrl->CM4->CM4_F.n_field_components; jj++)

test/mgd77/igrf_band.sh

@@ -0,0 +1,61 @@
+#!/usr/bin/env bash
+#
+# Tests the IGRF harmonic degree band limiting, i.e. gmt mgd77magref -F.../0+l<low>/<high>
+#
+# The X,Y,Z components are linear in the spherical harmonic coefficients, so any split of
+# the degree range must add back up to the unrestricted field.  That is checked here for
+# two different splits, plus a few pinned values and the option error checking.
+
+log=igrf_band.log
+
+cat << EOF > pts.txt
+-20	35	0	2020.0
+110	-45	10	2020.5
+0	90	0	2005.0
+30	-10	400	1980.0
+EOF
+
+# -Vq silences the warning issued when the requested degree 13 is truncated to 10 for the
+# pre-1995 models, which is itself tested further down.  We ask for all the digits we can
+# get since the residuals below are computed from these printed values.
+fmt=--FORMAT_FLOAT_OUT=%.15e
+gmt mgd77magref pts.txt -A+y -Fxyz/0       -Vq $fmt > full.dat
+gmt mgd77magref pts.txt -A+y -Fxyz/0+l1/5  -Vq $fmt > lo.dat
+gmt mgd77magref pts.txt -A+y -Fxyz/0+l6/13 -Vq $fmt > hi.dat
+gmt mgd77magref pts.txt -A+y -Fxyz/0+l1/1  -Vq $fmt > dip.dat
+gmt mgd77magref pts.txt -A+y -Fxyz/0+l2/13 -Vq $fmt > nodip.dat
+
+echo "Degrees 1-5 plus 6-13 must equal all degrees:" > $log
+paste lo.dat hi.dat full.dat | $AWK '{for (j = 1; j <= 3; j++) {d = $j + $(j+3) - $(j+6); if (d < 0) d = -d; if (d > max) max = d}}
+	END {printf "%s\n", (max < 1e-6) ? "residual < 1e-6 nT" : "FAIL: residual is " max " nT"}' >> $log
+echo "Degree 1 plus degrees 2-13 must equal all degrees:" >> $log
+paste dip.dat nodip.dat full.dat | $AWK '{for (j = 1; j <= 3; j++) {d = $j + $(j+3) - $(j+6); if (d < 0) d = -d; if (d > max) max = d}}
+	END {printf "%s\n", (max < 1e-6) ? "residual < 1e-6 nT" : "FAIL: residual is " max " nT"}' >> $log
+
+# Asking for all the degrees the model has must reproduce the unrestricted field exactly
+echo "Explicitly asking for degrees 1-13 must be the same as not restricting at all:" >> $log
+gmt mgd77magref pts.txt -A+y -Fxyz/0+l1/13 -Vq $fmt | diff - full.dat --strip-trailing-cr > /dev/null && echo "identical" >> $log
+
+# Pinned values (F, X, Y, Z) at lon = -20, lat = 35, sea level, 2020.0
+echo "Total field and components for the full field, degree 1 only, and degrees 6-13:" >> $log
+for band in "" "+l1/1" "+l6/13"; do
+	echo -20 35 0 2020.0 | gmt mgd77magref -A+y -Ftxyz/0$band --FORMAT_FLOAT_OUT=%.4f >> $log
+done
+
+# The modifier only applies to the IGRF, so it must be refused for any CM4 combination
+echo "Number of errors when +l is used with CM4 (-Ft/1) or with the mixed mode (-Ft/934):" >> $log
+gmt mgd77magref pts.txt -A+y -Ft/1+l1/5   2>&1 | grep -c "only available for the IGRF model" >> $log
+gmt mgd77magref pts.txt -A+y -Ft/934+l1/5 2>&1 | grep -c "only available for the IGRF model" >> $log
+echo "Number of errors for a malformed band (+l5), a reversed band (+l5/2) and a zero degree (+l0/3):" >> $log
+gmt mgd77magref pts.txt -A+y -Ft/0+l5   2>&1 | grep -c "usage is +l<low>/<high>" >> $log
+gmt mgd77magref pts.txt -A+y -Ft/0+l5/2 2>&1 | grep -c "requires 1 <= <low> <= <high>" >> $log
+gmt mgd77magref pts.txt -A+y -Ft/0+l0/3 2>&1 | grep -c "requires 1 <= <low> <= <high>" >> $log
+
+# The IGRF only has degrees up to 10 before 1995, so a higher upper degree must be truncated
+echo "Number of truncation warnings for a 1980 date with +l1/13:" >> $log
+echo 30 -10 400 1980.0 | gmt mgd77magref -A+y -Ft/0+l1/13 2>&1 | grep -c "only has harmonic degrees up to 10" >> $log
+echo "Truncating degree 13 to 10 for 1980 must give the complete 1980 field:" >> $log
+echo 30 -10 400 1980.0 | gmt mgd77magref -A+y -Ftxyz/0+l1/13 -Vq --FORMAT_FLOAT_OUT=%.4f > trunc.dat
+echo 30 -10 400 1980.0 | gmt mgd77magref -A+y -Ftxyz/0 --FORMAT_FLOAT_OUT=%.4f | diff - trunc.dat --strip-trailing-cr > /dev/null && echo "identical" >> $log
+
+diff $log "${src:-.}"/igrf_band.txt --strip-trailing-cr | tee fail > /dev/null

test/mgd77/igrf_band.txt

@@ -0,0 +1,21 @@
+Degrees 1-5 plus 6-13 must equal all degrees:
+residual < 1e-6 nT
+Degree 1 plus degrees 2-13 must equal all degrees:
+residual < 1e-6 nT
+Explicitly asking for degrees 1-13 must be the same as not restricting at all:
+identical
+Total field and components for the full field, degree 1 only, and degrees 6-13:
+41836.4995	28009.0175	-2809.0860	30949.9058
+44677.9387	22573.0698	-3876.4283	38360.7616
+222.4216	128.1612	170.1325	-64.0390
+Number of errors when +l is used with CM4 (-Ft/1) or with the mixed mode (-Ft/934):
+1
+1
+Number of errors for a malformed band (+l5), a reversed band (+l5/2) and a zero degree (+l0/3):
+1
+1
+1
+Number of truncation warnings for a 1980 date with +l1/13:
+1
+Truncating degree 13 to 10 for 1980 must give the complete 1980 field:
+identical