diff --git a/doc/rst/source/changes.rst b/doc/rst/source/changes.rst index d1796311473..b515a97f1b2 100644 --- a/doc/rst/source/changes.rst +++ b/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 ========================= diff --git a/doc/rst/source/supplements/mgd77/mgd77magref.rst b/doc/rst/source/supplements/mgd77/mgd77magref.rst index 015f6a38bad..b76d9797a40 100644 --- a/doc/rst/source/supplements/mgd77/mgd77magref.rst +++ b/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** diff --git a/src/mgd77/mgd77.c b/src/mgd77/mgd77.c index c9d296fb167..6a7e8ad4ba5 100644 --- a/src/mgd77/mgd77.c +++ b/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 */ diff --git a/src/mgd77/mgd77.h b/src/mgd77/mgd77.h index e6970f5ffaf..199528500fd 100644 --- a/src/mgd77/mgd77.h +++ b/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); diff --git a/src/mgd77/mgd77magref.c b/src/mgd77/mgd77magref.c index d3eff87d3ae..e53fdd3896c 100644 --- a/src/mgd77/mgd77magref.c +++ b/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/ 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 [] [-A+a+t+y] [-C] [-D] [-E] " - "[-Frthxyzdi[/[0|9]1234567]] [-G] [-Lrtxyz[/1234]] [-Sc|l/] [%s] " + "[-Frthxyzdi[/[0|9]1234567][+l/]] [-G] [-Lrtxyz[/1234]] [-Sc|l/] [%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/]"); 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 / 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/\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 <= <= \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++) diff --git a/test/mgd77/igrf_band.sh b/test/mgd77/igrf_band.sh new file mode 100644 index 00000000000..f5923139251 --- /dev/null +++ b/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/ +# +# 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/" >> $log +gmt mgd77magref pts.txt -A+y -Ft/0+l5/2 2>&1 | grep -c "requires 1 <= <= " >> $log +gmt mgd77magref pts.txt -A+y -Ft/0+l0/3 2>&1 | grep -c "requires 1 <= <= " >> $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 diff --git a/test/mgd77/igrf_band.txt b/test/mgd77/igrf_band.txt new file mode 100644 index 00000000000..705dbd99ad6 --- /dev/null +++ b/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