diff --git a/doc/rst/source/grdfft.rst b/doc/rst/source/grdfft.rst index ec9eaf836c5..3674b5bacff 100644 --- a/doc/rst/source/grdfft.rst +++ b/doc/rst/source/grdfft.rst @@ -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** ] @@ -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_ diff --git a/doc/rst/source/supplements/potential/gravfft.rst b/doc/rst/source/supplements/potential/gravfft.rst index beed5a24b17..00030a4b5be 100644 --- a/doc/rst/source/supplements/potential/gravfft.rst +++ b/doc/rst/source/supplements/potential/gravfft.rst @@ -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| ] @@ -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_ diff --git a/src/grdfft.c b/src/grdfft.c index f4138f29952..03ae0512249 100644 --- a/src/grdfft.c +++ b/src/grdfft.c @@ -85,6 +85,10 @@ struct GRDFFT_CTRL { bool active; double value; } I; + struct GRDFFT_M { /* -M */ + bool active; + double mgal_at_45; + } M; struct GRDFFT_N { /* -N[f|q|s/][+e|m|n][+t][+w[]][+z[p]] */ bool active; struct GMT_FFT_INFO *info; @@ -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 */ @@ -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 [] [-A] [-C] [-D[|g]] " - "[-E[r|x|y][+n][+w[k]]] [-F[r|x|y]] [-G|] [-I[|g]] [-N%s] [-Q] " + "[-E[r|x|y][+n][+w[k]]] [-F[r|x|y]] [-G|
] [-I[|g]] [-M] [-N%s] [-Q] " "[-S|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); @@ -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[|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"); + 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|d"); @@ -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); @@ -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; @@ -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; @@ -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); diff --git a/src/potential/gravfft.c b/src/potential/gravfft.c index 5945fb12c81..b3ed4f513b5 100644 --- a/src/potential/gravfft.c +++ b/src/potential/gravfft.c @@ -91,6 +91,10 @@ struct GRAVFFT_CTRL { bool active; double value; } I; + struct GRAVFFT_M { /* -M */ + bool active; + double mgal_at_45; + } M; struct GRAVFFT_N { /* -N[f|q|s/][+e|m|n][+t][+w[]][+z[p]] */ bool active; struct GMT_FFT_INFO *info; @@ -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 */ @@ -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); } @@ -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 [] -G%s [-C] " - "[-D] [-E] [-F[f[+s|z]|b|g|v|n|e]] [-I] [-N%s] [-Q] [-S] " + "[-D] [-E] [-F[f[+s|z]|b|g|v|n|e]] [-I] [-M] [-N%s] [-Q] [-S] " "[-T[/][+m]] [%s] [-W[k]] [-Z[/]] [-fg] [%s]\n", name, GMT_INGRID, GMT_OUTGRID, GMT_FFT_OPT, GMT_V_OPT, GMT_PAR_OPT); @@ -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"); + 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) " @@ -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; @@ -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; @@ -921,7 +932,7 @@ 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]); @@ -929,7 +940,7 @@ GMT_LOCAL void gravfft_do_parker (struct GMT_CTRL *GMT, struct GMT_GRID *Grid, s 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]); @@ -1214,7 +1225,7 @@ 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]); @@ -1222,7 +1233,7 @@ GMT_LOCAL void gravfft_load_from_top_grid (struct GMT_CTRL *GMT, struct GMT_GRID 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]);