Skip to content

Commit b3580a7

Browse files
joa-quimjoa-quim
authored
Update to #8699. There was another #def variable holding Earth's grav. (#8703)
Integrate it too in the -M option.
1 parent e1f3593 commit b3580a7

File tree

2 files changed

+17
-14
lines changed

2 files changed

+17
-14
lines changed

‎src/grdfft.c‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -227,7 +227,6 @@ GMT_LOCAL unsigned int grdfft_do_azimuthal_derivative (struct GMT_GRID *Grid, do
227227
/* Now obsolete but left for backwards compatibility. Users are encouraged to use gravfft */
228228
#define POISSONS_RATIO 0.25
229229
#define YOUNGS_MODULUS 1.0e11 /* Pascal = Nt/m**2 */
230-
#define NORMAL_GRAVITY 9.806199203 /* m/s**2 */
231230

232231
GMT_LOCAL unsigned int grdfft_do_isostasy (struct GMT_GRID *Grid, struct GRDFFT_CTRL *Ctrl, double *par, struct GMT_FFT_WAVENUMBER *K) {
233232
/* Do the isostatic response function convolution in the Freq domain.
@@ -255,7 +254,7 @@ GMT_LOCAL unsigned int grdfft_do_isostasy (struct GMT_GRID *Grid, struct GRDFFT_
255254
}
256255

257256
rigidity_d = (YOUNGS_MODULUS * pow (te, 3.0)) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
258-
d_over_restoring_force = rigidity_d / ((rm - ri) * NORMAL_GRAVITY);
257+
d_over_restoring_force = rigidity_d / ((rm - ri) * Ctrl->M.mgal_at_45 * 1e-5);
259258

260259
for (k = 0; k < Grid->header->size; k += 2) {
261260
mk = gmt_fft_get_wave (k, K);

‎src/potential/gravfft.c‎

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -93,7 +93,7 @@ struct GRAVFFT_CTRL {
9393
} I;
9494
struct GRAVFFT_M { /* -M<mgal_at_45> */
9595
bool active;
96-
double mgal_at_45;
96+
double mgal_at_45, normal_gravity;
9797
} M;
9898
struct GRAVFFT_N { /* -N[f|q|s<n_columns>/<n_rows>][+e|m|n][+t<width>][+w[<suffix>]][+z[p]] */
9999
bool active;
@@ -137,7 +137,7 @@ struct GRAVFFT_CTRL {
137137
#endif
138138

139139
#define YOUNGS_MODULUS 7.0e10 /* Pascal = Nt/m**2 */
140-
#define NORMAL_GRAVITY 9.806199203 /* Moritz's 1980 IGF value for gravity at 45 degrees latitude (m/s) */
140+
//#define NORMAL_GRAVITY 9.806199203 /* Moritz's 1980 IGF value for gravity at 45 degrees latitude (m/s) */
141141
#define POISSONS_RATIO 0.25
142142

143143
static bool sphericity = false;
@@ -322,6 +322,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
322322
case 'M': /* Miligals at 45 degrees */
323323
n_errors += gmt_M_repeated_module_option(API, Ctrl->M.active);
324324
sscanf(opt->arg, "%lf", &Ctrl->M.mgal_at_45);
325+
Ctrl->M.normal_gravity = Ctrl->M.mgal_at_45 * 1e-5;
325326
break;
326327
case 'N':
327328
n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
@@ -423,7 +424,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVFFT_CTRL *Ctrl, struct GMT_OP
423424
!(Ctrl->F.mode == GRAVFFT_FAA || Ctrl->F.mode == GRAVFFT_GEOID),
424425
"Theoretical admittances are only defined for FAA or GEOID.\n");
425426

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 */
427+
if (!Ctrl->M.active) {
428+
Ctrl->M.mgal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
429+
Ctrl->M.normal_gravity = 9.806199203;
430+
}
427431

428432
return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
429433
}
@@ -875,7 +879,7 @@ GMT_LOCAL void gravfft_do_isostasy (struct GMT_CTRL *GMT, struct GMT_GRID *Grid,
875879
}
876880
GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Using effective load density rho_l = %g and Airy boost factor A = %g\n", rho_load, A);
877881
rigidity_d = (YOUNGS_MODULUS * Ctrl->T.te * Ctrl->T.te * Ctrl->T.te) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
878-
d_over_restoring_force = rigidity_d / ( (Ctrl->T.rhom - rho_load) * NORMAL_GRAVITY);
882+
d_over_restoring_force = rigidity_d / ( (Ctrl->T.rhom - rho_load) * Ctrl->M.normal_gravity);
879883
airy_ratio = -A * (rho_load - Ctrl->T.rhow)/(Ctrl->T.rhom - rho_load);
880884

881885
if (Ctrl->T.te == 0.0) { /* Airy isostasy; scale global variable scale_out and return */
@@ -1120,7 +1124,7 @@ GMT_LOCAL void gravfft_load_from_below_admitt(struct GMT_CTRL *GMT, struct GRAVF
11201124

11211125
delta_k /= TWO_PI; /* Use frequency, not wavenumber */
11221126
D = (YOUNGS_MODULUS * Ctrl->T.te * Ctrl->T.te * Ctrl->T.te) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
1123-
alfa = pow (TWO_PI, 4.) * D / (NORMAL_GRAVITY * Ctrl->T.rho_mc);
1127+
alfa = pow (TWO_PI, 4.) * D / (Ctrl->M.normal_gravity * Ctrl->T.rho_mc);
11241128

11251129
for (k = 0; k < nk; k++) {
11261130
freq = (k + 1) * delta_k;
@@ -1129,7 +1133,7 @@ GMT_LOCAL void gravfft_load_from_below_admitt(struct GMT_CTRL *GMT, struct GRAVF
11291133
if (Ctrl->F.mode == GRAVFFT_FAA)
11301134
t1 *= 1.0e5; /* to have it in mGals */
11311135
else /* Must be the GEOID case */
1132-
t1 /= (NORMAL_GRAVITY * freq * TWO_PI);
1136+
t1 /= (Ctrl->M.normal_gravity * freq * TWO_PI);
11331137
t2 = Ctrl->T.rho_cw * exp(-TWO_PI * freq * Ctrl->misc.z_level) +
11341138
Ctrl->T.rho_mc * exp(-TWO_PI * freq * Ctrl->Z.zm);
11351139
t3 = -(Ctrl->T.rho_mw + Ctrl->T.rho_mc * pow(freq,4.) * alfa) * exp(-TWO_PI * freq * Ctrl->Z.zl);
@@ -1156,7 +1160,7 @@ GMT_LOCAL void gravfft_load_from_top_admitt (struct GMT_CTRL *GMT, struct GRAVFF
11561160

11571161
delta_k /= TWO_PI; /* Use frequency, not wavenumber */
11581162
D = (YOUNGS_MODULUS * Ctrl->T.te * Ctrl->T.te * Ctrl->T.te) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
1159-
alfa = pow (TWO_PI,4.) * D / (NORMAL_GRAVITY * Ctrl->T.rho_mc);
1163+
alfa = pow (TWO_PI,4.) * D / (Ctrl->M.normal_gravity * Ctrl->T.rho_mc);
11601164

11611165
for (k = 0; k < nk; k++) {
11621166
freq = (k + 1) * delta_k;
@@ -1165,7 +1169,7 @@ GMT_LOCAL void gravfft_load_from_top_admitt (struct GMT_CTRL *GMT, struct GRAVFF
11651169
if (Ctrl->F.mode == GRAVFFT_FAA)
11661170
t1 *= 1.0e5; /* to have it in mGals */
11671171
else /* Must be the GEOID case */
1168-
t1 /= (NORMAL_GRAVITY * freq * TWO_PI);
1172+
t1 /= (Ctrl->M.normal_gravity * freq * TWO_PI);
11691173
t2 = exp(-TWO_PI * freq * Ctrl->misc.z_level) - exp(-TWO_PI * freq * Ctrl->Z.zm) / (1 + alfa*pow(freq,4.));
11701174
z_from_top[k] = t1 * Ctrl->T.rho_cw * t2;
11711175
}
@@ -1188,7 +1192,7 @@ GMT_LOCAL void gravfft_load_from_top_grid (struct GMT_CTRL *GMT, struct GMT_GRID
11881192
p = n - 1.0;
11891193

11901194
D = (YOUNGS_MODULUS * Ctrl->T.te * Ctrl->T.te * Ctrl->T.te) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
1191-
alfa = pow(TWO_PI,4.) * D / (NORMAL_GRAVITY * Ctrl->T.rho_mc);
1195+
alfa = pow(TWO_PI,4.) * D / (Ctrl->M.normal_gravity * Ctrl->T.rho_mc);
11921196
raised[0] = 0.0f; raised[1] = 0.0f;
11931197

11941198
for (k = 0; k < Grid->header->size; k+= 2) {
@@ -1213,7 +1217,7 @@ GMT_LOCAL void gravfft_load_from_top_grid (struct GMT_CTRL *GMT, struct GMT_GRID
12131217
datac[k+1] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);
12141218
break;
12151219
case GRAVFFT_GEOID:
1216-
if (mk > 0.0) t1 /= (NORMAL_GRAVITY * mk * TWO_PI);
1220+
if (mk > 0.0) t1 /= (Ctrl->M.normal_gravity * mk * TWO_PI);
12171221
datac[k] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k]);
12181222
datac[k+1] += (gmt_grdfloat) ((Ctrl->T.rho_cw * t1 * t2) * t / f * raised[k+1]);
12191223
break;
@@ -1259,7 +1263,7 @@ GMT_LOCAL void gravfft_load_from_below_grid (struct GMT_CTRL *GMT, struct GMT_GR
12591263
p = n - 1.0;
12601264

12611265
D = (YOUNGS_MODULUS * Ctrl->T.te * Ctrl->T.te * Ctrl->T.te) / (12.0 * (1.0 - POISSONS_RATIO * POISSONS_RATIO));
1262-
alfa = pow(TWO_PI,4.) * D / (NORMAL_GRAVITY * Ctrl->T.rho_mc);
1266+
alfa = pow(TWO_PI,4.) * D / (Ctrl->M.normal_gravity * Ctrl->T.rho_mc);
12631267
raised[0] = 0.0f; raised[1] = 0.0f;
12641268

12651269
for (k = 0; k < Grid->header->size; k+= 2) {
@@ -1275,7 +1279,7 @@ GMT_LOCAL void gravfft_load_from_below_grid (struct GMT_CTRL *GMT, struct GMT_GR
12751279
if (Ctrl->F.mode == GRAVFFT_FAA)
12761280
t1 *= 1.0e5; /* to have it in mGals */
12771281
else /* Must be the GEOID case */
1278-
t1 /= (NORMAL_GRAVITY * mk * TWO_PI);
1282+
t1 /= (Ctrl->M.normal_gravity * mk * TWO_PI);
12791283
t2 = Ctrl->T.rho_cw * exp(-TWO_PI * mk * Ctrl->misc.z_level) +
12801284
Ctrl->T.rho_mc * exp(-TWO_PI * mk * Ctrl->Z.zm);
12811285
t3 = -(Ctrl->T.rho_mw + Ctrl->T.rho_mc * pow(mk,4.) * alfa) * exp(-TWO_PI * mk * Ctrl->Z.zl);

0 commit comments

Comments
 (0)