Skip to content

Commit 77d4d33

Browse files
joa-quimjoa-quim
committed
Add a new option to grdfft and gravfft to let change the MGAL_AT_45
See this Forum post. https://forum.generic-mapping-tools.org/t/gravitational-admittance-estimate-for-mars/5776/4
1 parent b7fe1ad commit 77d4d33

File tree

4 files changed

+64
-21
lines changed

4 files changed

+64
-21
lines changed

‎doc/rst/source/grdfft.rst

+8
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ Synopsis
2020
[ |-E|\ [**r**\|\ **x**\|\ **y**][**+n**][**+w**\ [**k**]] ]
2121
[ |-F|\ [**r**\|\ **x**\|\ **y**]\ *params* ]
2222
[ |-I|\ [*scale*\|\ **g**] ]
23+
[ |-M|\ *mgal_at_45* ]
2324
[ |-N|\ *params* ]
2425
[ |-Q|\ ]
2526
[ |-S|\ *scale*\|\ **d** ]
@@ -164,6 +165,13 @@ Optional Arguments
164165
is gravity anomalies in mGal and output should be geoid heights in
165166
meters. Repeatable. [Default is no scale].
166167

168+
.. _-M:
169+
170+
**-M**\ *mgal_at_45*
171+
Specify the value of the gravity in mili Gals at 45 degrees latitude (used to convert
172+
gravity anomalies to geoid heights). Default is 980619.9203 mGal (Moritz's 1980 IGF value).
173+
This value needs to be changed accordingly when using data from other planets.
174+
167175
.. _-N:
168176

169177
.. include:: explain_fft.rst_

‎doc/rst/source/supplements/potential/gravfft.rst

+8
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ Synopsis
1919
[ |-E|\ *n_terms* ]
2020
[ |-F|\ [**f**\ [**+s**\|\ **z**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**] ]
2121
[ |-I|\ **w**\|\ **b**\|\ **c**\|\ **t**\|\ **k** ]
22+
[ |-M|\ *mgal_at_45* ]
2223
[ |-N|\ *params* ]
2324
[ |-Q| ]
2425
[ |-S| ]
@@ -138,6 +139,13 @@ Optional Arguments
138139
theoretical admittance, and **t** writes a fourth column with "elastic
139140
plate" theoretical admittance.
140141

142+
.. _-M:
143+
144+
**-M**\ *mgal_at_45*
145+
Specify the value of the gravity in mili Gals at 45 degrees latitude (used to convert
146+
gravity anomalies to geoid heights). Default is 980619.9203 mGal (Moritz's 1980 IGF value).
147+
This value needs to be changed accordingly when using data from other planets.
148+
141149
.. _-N:
142150

143151
.. include:: ../../explain_fft.rst_

‎src/grdfft.c

+30-14
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,10 @@ struct GRDFFT_CTRL {
8585
bool active;
8686
double value;
8787
} I;
88+
struct GRDFFT_M { /* -M<mgal_at_45> */
89+
bool active;
90+
double mgal_at_45;
91+
} M;
8892
struct GRDFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
8993
bool active;
9094
struct GMT_FFT_INFO *info;
@@ -115,8 +119,6 @@ enum Grdfft_operators {
115119
GRDFFT_ISOSTASY ,
116120
GRDFFT_NOTHING };
117121

118-
#define MGAL_AT_45 980619.9203 /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
119-
120122
struct F_INFO {
121123
double lc[3]; /* Low-cut frequency for r, x, and y */
122124
double lp[3]; /* Low-pass frequency for r, x, and y */
@@ -571,7 +573,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
571573
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
572574
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
573575
GMT_Usage (API, 0, "usage: %s %s [<ingrid2>] [-A<azimuth>] [-C<zlevel>] [-D[<scale>|g]] "
574-
"[-E[r|x|y][+n][+w[k]]] [-F[r|x|y]<parameters>] [-G<outgrid>|<table>] [-I[<scale>|g]] [-N%s] [-Q] "
576+
"[-E[r|x|y][+n][+w[k]]] [-F[r|x|y]<parameters>] [-G<outgrid>|<table>] [-I[<scale>|g]] [-M<mgal_at_45>] [-N%s] [-Q] "
575577
"[-S<scale>|d] [%s] [-fg] [%s] [%s]\n", name, GMT_INGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_ho_OPT, GMT_PAR_OPT);
576578

577579
if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
@@ -609,6 +611,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
609611
"Note: Optional for -E (spectrum written to standard output); required otherwise.");
610612
GMT_Usage (API, 1, "\n-I[<scale>|g]");
611613
GMT_Usage (API, -2, "Integrate, i.e., divide by kr [ and optionally by scale]. Use -Ig to get m from mGal] (repeatable).");
614+
GMT_Usage (API, 1, "\n-M<mgal_at_45>");
615+
GMT_Usage (API, -2, "Value for gravity in mGal at 45 degrees latitude. Default is 980619.9203 mGal (Moritz's 1980 IGF value).");
612616
GMT_FFT_Option (API, 'N', GMT_FFT_DIM, "Choose or inquire about suitable grid dimensions for FFT, and set modifiers.");
613617
GMT_Usage (API, 1, "\n-Q Perform no operations, just do forward FFF and write output if selected in -N.");
614618
GMT_Usage (API, 1, "\n-S<scale>|d");
@@ -665,6 +669,27 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
665669
f_info->hp[j] = f_info->hc[j] = DBL_MAX; /* Set huge positive, above valid frequency range */
666670
}
667671

672+
for (opt = options; opt; opt = opt->next) { /* Process -M first because we may need its value for -D and/or -I */
673+
switch (opt->option) {
674+
case 'M': /* Geographic data */
675+
if (gmt_M_compat_check (GMT, 4)) {
676+
if (sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45))
677+
Ctrl->M.active = true;
678+
else {
679+
GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
680+
if (gmt_M_is_cartesian (GMT, GMT_IN))
681+
gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
682+
}
683+
}
684+
else {
685+
n_errors += gmt_M_repeated_module_option(API, Ctrl->M.active);
686+
sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45);
687+
}
688+
break;
689+
}
690+
}
691+
if (!Ctrl->M.active) Ctrl->M.mgal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
692+
668693
for (opt = options; opt; opt = opt->next) { /* Process all the options given */
669694
gmt_M_memset (par, 5, double);
670695

@@ -697,7 +722,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
697722
break;
698723
case 'D': /* d/dz [repeatable] */
699724
Ctrl->D.active = true;
700-
par[0] = (opt->arg[0]) ? ((opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg)) : 1.0;
725+
par[0] = (opt->arg[0]) ? ((opt->arg[0] == 'g' || opt->arg[0] == 'G') ? Ctrl->M.mgal_at_45 : atof (opt->arg)) : 1.0;
701726
n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -D: scale must be nonzero\n");
702727
grdfft_add_operation (GMT, Ctrl, GRDFFT_DIFFERENTIATE, 1, par);
703728
break;
@@ -750,7 +775,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
750775
break;
751776
case 'I': /* Integrate [repeatable]*/
752777
Ctrl->I.active = true;
753-
par[0] = (opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg);
778+
par[0] = (opt->arg[0] == 'g' || opt->arg[0] == 'G') ? Ctrl->M.mgal_at_45 : atof (opt->arg);
754779
n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -I: scale must be nonzero\n");
755780
grdfft_add_operation (GMT, Ctrl, GRDFFT_INTEGRATE, 1, par);
756781
break;
@@ -760,15 +785,6 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
760785
else
761786
n_errors += gmt_default_option_error (GMT, opt);
762787
break;
763-
case 'M': /* Geographic data */
764-
if (gmt_M_compat_check (GMT, 4)) {
765-
GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
766-
if (gmt_M_is_cartesian (GMT, GMT_IN))
767-
gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
768-
}
769-
else
770-
n_errors += gmt_default_option_error (GMT, opt);
771-
break;
772788
case 'N': /* Grid dimension setting or inquiery */
773789
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
774790
if (ptr && gmt_M_compat_check (GMT, 4)) { /* Got both old -L and -N; append */

‎src/potential/gravfft.c

+18-7
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,10 @@ struct GRAVFFT_CTRL {
9191
bool active;
9292
double value;
9393
} I;
94+
struct GRAVFFT_M { /* -M<mgal_at_45> */
95+
bool active;
96+
double mgal_at_45;
97+
} M;
9498
struct GRAVFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
9599
bool active;
96100
struct GMT_FFT_INFO *info;
@@ -315,6 +319,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
315319
else
316320
n_errors += gmt_default_option_error (GMT, opt);
317321
break;
322+
case 'M': /* Miligals at 45 degrees */
323+
n_errors += gmt_M_repeated_module_option(API, Ctrl->M.active);
324+
sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45);
325+
break;
318326
case 'N':
319327
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
320328
if (popt && gmt_M_compat_check (GMT, 4)) { /* Got both old -L and -N; append */
@@ -414,6 +422,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
414422
n_errors += gmt_M_check_condition (GMT, (Ctrl->misc.from_top || Ctrl->misc.from_below) &&
415423
!(Ctrl->F.mode == GRAVFFT_FAA || Ctrl->F.mode == GRAVFFT_GEOID),
416424
"Theoretical admittances are only defined for FAA or GEOID.\n");
425+
426+
if (!Ctrl->M.active) Ctrl->M.mgal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
417427

418428
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
419429
}
@@ -422,7 +432,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
422432
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
423433
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
424434
GMT_Usage (API, 0, "usage: %s %s [<ingrid2>] -G%s [-C<n/wavelength/mean_depth/tbw>] "
425-
"[-D<density>] [-E<n_terms>] [-F[f[+s|z]|b|g|v|n|e]] [-I<cbktw>] [-N%s] [-Q] [-S] "
435+
"[-D<density>] [-E<n_terms>] [-F[f[+s|z]|b|g|v|n|e]] [-I<cbktw>] [-M<mgal_at_45>] [-N%s] [-Q] [-S] "
426436
"[-T<te/rl/rm/rw>[/<ri>][+m]] [%s] [-W<wd>[k]] [-Z<zm>[/<zl>]] [-fg] [%s]\n",
427437
name, GMT_INGRID, GMT_OUTGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_PAR_OPT);
428438

@@ -468,6 +478,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
468478
GMT_Usage (API, 3, "k: Use km or wavelength unit [m].");
469479
GMT_Usage (API, 3, "t: Write a forth column with \"elastic plate\" theoretical admittance.");
470480
GMT_Usage (API, 3, "w: Write wavelength instead of wavenumber.");
481+
GMT_Usage (API, 1, "\n-M<mgal_at_45>");
482+
GMT_Usage (API, -2, "Value for gravity in mGal at 45 degrees latitude. Default is 980619.9203 mGal (Moritz's 1980 IGF value).");
471483
GMT_FFT_Option (API, 'N', GMT_FFT_DIM, "Choose or inquire about suitable grid dimensions for FFT, and set modifiers.");
472484
GMT_Usage (API, -2, "Warning: both -D -T...+m and -Q will implicitly set -N's +h.");
473485
GMT_Usage (API, 1, "\n-Q Writes out a grid with the flexural topography (with z positive up) "
@@ -880,7 +892,6 @@ GMT_LOCAL void gravfft_do_isostasy (struct GMT_CTRL *GMT, struct GMT_GRID *Grid,
880892
}
881893
}
882894

883-
#define MGAL_AT_45 980619.9203 /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
884895
GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, struct GRAVFFT_CTRL *Ctrl, struct GMT_FFT_WAVENUMBER *K, gmt_grdfloat *raised, uint64_t n, double rho) {
885896
uint64_t i, k;
886897
double f, p, t, mk, kx, ky, v, c;
@@ -909,7 +920,7 @@ GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, s
909920
datac[k+1] += (gmt_grdfloat) (v * raised[k+1]);
910921
break;
911922
case GRAVFFT_GEOID:
912-
if (mk > 0.0) v /= (MGAL_AT_45 * mk);
923+
if (mk > 0.0) v /= (Ctrl->M.mgal_at_45 * mk);
913924
datac[k] += (gmt_grdfloat) (v * raised[k]);
914925
datac[k+1] += (gmt_grdfloat) (v * raised[k+1]);
915926
break;
@@ -921,15 +932,15 @@ GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, s
921932
case GRAVFFT_DEFL_EAST:
922933
if (mk > 0.0) { /* Scale tan (xslope) ~ slope to microradians */
923934
kx = gmt_fft_any_wave (k, GMT_FFT_K_IS_KX, K);
924-
v *= 1.e6 * (-kx / (MGAL_AT_45 * mk));
935+
v *= 1.e6 * (-kx / (Ctrl->M.mgal_at_45 * mk));
925936
}
926937
datac[k] += (gmt_grdfloat) (-v * raised[k+1]);
927938
datac[k+1] += (gmt_grdfloat) ( v * raised[k]);
928939
break;
929940
case GRAVFFT_DEFL_NORTH:
930941
if (mk > 0.0) { /* Scale tan (yslope) ~ slope to microradians */
931942
ky = gmt_fft_any_wave (k, GMT_FFT_K_IS_KY, K);
932-
v *= 1.e6 * (-ky / (MGAL_AT_45 * mk));
943+
v *= 1.e6 * (-ky / (Ctrl->M.mgal_at_45 * mk));
933944
}
934945
datac[k] += (gmt_grdfloat) ( v * raised[k+1]);
935946
datac[k+1] += (gmt_grdfloat) (-v * raised[k]);
@@ -1214,15 +1225,15 @@ GMT_LOCAL void gravfft_load_from_top_grid (struct GMT_CTRL *GMT, struct GMT_GRID
12141225
case GRAVFFT_DEFL_EAST:
12151226
if (mk > 0.0) { /* Scale tan (xslope) ~ slope to microradians */
12161227
double kx = gmt_fft_any_wave (k, GMT_FFT_K_IS_KX, K);
1217-
t1 *= 1.e6 * (-kx / (MGAL_AT_45 * mk));
1228+
t1 *= 1.e6 * (-kx / (Ctrl->M.mgal_at_45 * mk));
12181229
}
12191230
datac[k] += (gmt_grdfloat) (-(Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k]);
12201231
datac[k+1] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);
12211232
break;
12221233
case GRAVFFT_DEFL_NORTH:
12231234
if (mk > 0.0) { /* Scale tan (yslope) ~ slope to microradians */
12241235
double ky = gmt_fft_any_wave (k, GMT_FFT_K_IS_KY, K);
1225-
t1 *= 1.e6 * (-ky / (MGAL_AT_45 * mk));
1236+
t1 *= 1.e6 * (-ky / (Ctrl->M.mgal_at_45 * mk));
12261237
}
12271238
datac[k] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k]);
12281239
datac[k+1] += (gmt_grdfloat) (-(Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);

0 commit comments

Comments
 (0)