Skip to content

Add a new option to grdfft and gravfft to let change the MGAL_AT_45 value #8699

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
joa-quim merged 2 commits into from
Mar 10, 2025
+65 −20
Merged

Add a new option to grdfft and gravfft to let change the MGAL_AT_45 value #8699

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
77d4d33
Add a new option to grdfft and gravfft to let change the MGAL_AT_45
joa-quim Mar 10, 2025
bc858e4
Need to skip -M parsing on the second loop
joa-quim Mar 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions doc/rst/source/grdfft.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ Synopsis
[ |-E|\ [**r**\|\ **x**\|\ **y**][**+n**][**+w**\ [**k**]] ]
[ |-F|\ [**r**\|\ **x**\|\ **y**]\ *params* ]
[ |-I|\ [*scale*\|\ **g**] ]
[ |-M|\ *mgal_at_45* ]
[ |-N|\ *params* ]
[ |-Q|\ ]
[ |-S|\ *scale*\|\ **d** ]
Expand Down Expand Up @@ -164,6 +165,13 @@ Optional Arguments
is gravity anomalies in mGal and output should be geoid heights in
meters. Repeatable. [Default is no scale].

.. _-M:

**-M**\ *mgal_at_45*
Specify the value of the gravity in mili Gals at 45 degrees latitude (used to convert
gravity anomalies to geoid heights). Default is 980619.9203 mGal (Moritz's 1980 IGF value).
This value needs to be changed accordingly when using data from other planets.

.. _-N:

.. include:: explain_fft.rst_
Expand Down
8 changes: 8 additions & 0 deletions doc/rst/source/supplements/potential/gravfft.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ Synopsis
[ |-E|\ *n_terms* ]
[ |-F|\ [**f**\ [**+s**\|\ **z**]\|\ **b**\|\ **g**\|\ **v**\|\ **n**\|\ **e**] ]
[ |-I|\ **w**\|\ **b**\|\ **c**\|\ **t**\|\ **k** ]
[ |-M|\ *mgal_at_45* ]
[ |-N|\ *params* ]
[ |-Q| ]
[ |-S| ]
Expand Down Expand Up @@ -138,6 +139,13 @@ Optional Arguments
theoretical admittance, and **t** writes a fourth column with "elastic
plate" theoretical admittance.

.. _-M:

**-M**\ *mgal_at_45*
Specify the value of the gravity in mili Gals at 45 degrees latitude (used to convert
gravity anomalies to geoid heights). Default is 980619.9203 mGal (Moritz's 1980 IGF value).
This value needs to be changed accordingly when using data from other planets.

.. _-N:

.. include:: ../../explain_fft.rst_
Expand Down
44 changes: 31 additions & 13 deletions src/grdfft.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,10 @@ struct GRDFFT_CTRL {
bool active;
double value;
} I;
struct GRDFFT_M { /* -M<mgal_at_45> */
bool active;
double mgal_at_45;
} M;
struct GRDFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
bool active;
struct GMT_FFT_INFO *info;
Expand Down Expand Up @@ -115,8 +119,6 @@ enum Grdfft_operators {
GRDFFT_ISOSTASY ,
GRDFFT_NOTHING };

#define MGAL_AT_45 980619.9203 /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */

struct F_INFO {
double lc[3]; /* Low-cut frequency for r, x, and y */
double lp[3]; /* Low-pass frequency for r, x, and y */
Expand Down Expand Up @@ -571,7 +573,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 %s [<ingrid2>] [-A<azimuth>] [-C<zlevel>] [-D[<scale>|g]] "
"[-E[r|x|y][+n][+w[k]]] [-F[r|x|y]<parameters>] [-G<outgrid>|<table>] [-I[<scale>|g]] [-N%s] [-Q] "
"[-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] "
"[-S<scale>|d] [%s] [-fg] [%s] [%s]\n", name, GMT_INGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_ho_OPT, GMT_PAR_OPT);

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

for (opt = options; opt; opt = opt->next) { /* Process -M first because we may need its value for -D and/or -I */
switch (opt->option) {
case 'M': /* Geographic data */
if (gmt_M_compat_check (GMT, 4)) {
if (sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45))
Ctrl->M.active = true;
else {
GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
if (gmt_M_is_cartesian (GMT, GMT_IN))
gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
}
}
else {
n_errors += gmt_M_repeated_module_option(API, Ctrl->M.active);
sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45);
}
break;
}
}
if (!Ctrl->M.active) Ctrl->M.mgal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */

for (opt = options; opt; opt = opt->next) { /* Process all the options given */
gmt_M_memset (par, 5, double);

Expand Down Expand Up @@ -697,7 +722,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
break;
case 'D': /* d/dz [repeatable] */
Ctrl->D.active = true;
par[0] = (opt->arg[0]) ? ((opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg)) : 1.0;
par[0] = (opt->arg[0]) ? ((opt->arg[0] == 'g' || opt->arg[0] == 'G') ? Ctrl->M.mgal_at_45 : atof (opt->arg)) : 1.0;
n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -D: scale must be nonzero\n");
grdfft_add_operation (GMT, Ctrl, GRDFFT_DIFFERENTIATE, 1, par);
break;
Expand Down Expand Up @@ -750,7 +775,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
break;
case 'I': /* Integrate [repeatable]*/
Ctrl->I.active = true;
par[0] = (opt->arg[0] == 'g' || opt->arg[0] == 'G') ? MGAL_AT_45 : atof (opt->arg);
par[0] = (opt->arg[0] == 'g' || opt->arg[0] == 'G') ? Ctrl->M.mgal_at_45 : atof (opt->arg);
n_errors += gmt_M_check_condition (GMT, par[0] == 0.0, "Option -I: scale must be nonzero\n");
grdfft_add_operation (GMT, Ctrl, GRDFFT_INTEGRATE, 1, par);
break;
Expand All @@ -760,14 +785,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDFFT_CTRL *Ctrl, struct F_INFO
else
n_errors += gmt_default_option_error (GMT, opt);
break;
case 'M': /* Geographic data */
if (gmt_M_compat_check (GMT, 4)) {
GMT_Report (API, GMT_MSG_COMPAT, "Option -M is deprecated; -fg was set instead, use this in the future.\n");
if (gmt_M_is_cartesian (GMT, GMT_IN))
gmt_parse_common_options (GMT, "f", 'f', "g"); /* Set -fg unless already set */
}
else
n_errors += gmt_default_option_error (GMT, opt);
case 'M': /* This was already parsed above */
break;
case 'N': /* Grid dimension setting or inquiery */
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
Expand Down
25 changes: 18 additions & 7 deletions src/potential/gravfft.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,10 @@ struct GRAVFFT_CTRL {
bool active;
double value;
} I;
struct GRAVFFT_M { /* -M<mgal_at_45> */
bool active;
double mgal_at_45;
} M;
struct GRAVFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
bool active;
struct GMT_FFT_INFO *info;
Expand Down Expand Up @@ -315,6 +319,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
else
n_errors += gmt_default_option_error (GMT, opt);
break;
case 'M': /* Miligals at 45 degrees */
n_errors += gmt_M_repeated_module_option(API, Ctrl->M.active);
sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45);
break;
case 'N':
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
if (popt && gmt_M_compat_check (GMT, 4)) { /* Got both old -L and -N; append */
Expand Down Expand Up @@ -414,6 +422,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
n_errors += gmt_M_check_condition (GMT, (Ctrl->misc.from_top || Ctrl->misc.from_below) &&
!(Ctrl->F.mode == GRAVFFT_FAA || Ctrl->F.mode == GRAVFFT_GEOID),
"Theoretical admittances are only defined for FAA or GEOID.\n");

if (!Ctrl->M.active) Ctrl->M.mgal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */

return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}
Expand All @@ -422,7 +432,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 %s [<ingrid2>] -G%s [-C<n/wavelength/mean_depth/tbw>] "
"[-D<density>] [-E<n_terms>] [-F[f[+s|z]|b|g|v|n|e]] [-I<cbktw>] [-N%s] [-Q] [-S] "
"[-D<density>] [-E<n_terms>] [-F[f[+s|z]|b|g|v|n|e]] [-I<cbktw>] [-M<mgal_at_45>] [-N%s] [-Q] [-S] "
"[-T<te/rl/rm/rw>[/<ri>][+m]] [%s] [-W<wd>[k]] [-Z<zm>[/<zl>]] [-fg] [%s]\n",
name, GMT_INGRID, GMT_OUTGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_PAR_OPT);

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

#define MGAL_AT_45 980619.9203 /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
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) {
uint64_t i, k;
double f, p, t, mk, kx, ky, v, c;
Expand Down Expand Up @@ -909,7 +920,7 @@ GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, s
datac[k+1] += (gmt_grdfloat) (v * raised[k+1]);
break;
case GRAVFFT_GEOID:
if (mk > 0.0) v /= (MGAL_AT_45 * mk);
if (mk > 0.0) v /= (Ctrl->M.mgal_at_45 * mk);
datac[k] += (gmt_grdfloat) (v * raised[k]);
datac[k+1] += (gmt_grdfloat) (v * raised[k+1]);
break;
Expand All @@ -921,15 +932,15 @@ GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, s
case GRAVFFT_DEFL_EAST:
if (mk > 0.0) { /* Scale tan (xslope) ~ slope to microradians */
kx = gmt_fft_any_wave (k, GMT_FFT_K_IS_KX, K);
v *= 1.e6 * (-kx / (MGAL_AT_45 * mk));
v *= 1.e6 * (-kx / (Ctrl->M.mgal_at_45 * mk));
}
datac[k] += (gmt_grdfloat) (-v * raised[k+1]);
datac[k+1] += (gmt_grdfloat) ( v * raised[k]);
break;
case GRAVFFT_DEFL_NORTH:
if (mk > 0.0) { /* Scale tan (yslope) ~ slope to microradians */
ky = gmt_fft_any_wave (k, GMT_FFT_K_IS_KY, K);
v *= 1.e6 * (-ky / (MGAL_AT_45 * mk));
v *= 1.e6 * (-ky / (Ctrl->M.mgal_at_45 * mk));
}
datac[k] += (gmt_grdfloat) ( v * raised[k+1]);
datac[k+1] += (gmt_grdfloat) (-v * raised[k]);
Expand Down Expand Up @@ -1214,15 +1225,15 @@ GMT_LOCAL void gravfft_load_from_top_grid (struct GMT_CTRL *GMT, struct GMT_GRID
case GRAVFFT_DEFL_EAST:
if (mk > 0.0) { /* Scale tan (xslope) ~ slope to microradians */
double kx = gmt_fft_any_wave (k, GMT_FFT_K_IS_KX, K);
t1 *= 1.e6 * (-kx / (MGAL_AT_45 * mk));
t1 *= 1.e6 * (-kx / (Ctrl->M.mgal_at_45 * mk));
}
datac[k] += (gmt_grdfloat) (-(Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k]);
datac[k+1] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);
break;
case GRAVFFT_DEFL_NORTH:
if (mk > 0.0) { /* Scale tan (yslope) ~ slope to microradians */
double ky = gmt_fft_any_wave (k, GMT_FFT_K_IS_KY, K);
t1 *= 1.e6 * (-ky / (MGAL_AT_45 * mk));
t1 *= 1.e6 * (-ky / (Ctrl->M.mgal_at_45 * mk));
}
datac[k] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k]);
datac[k+1] += (gmt_grdfloat) (-(Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);
Expand Down
Loading