From 01187361011874dbf97b54fe5ab84ac975084fe0 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sat, 12 Mar 2022 17:29:44 +0000 Subject: [PATCH 01/27] WIP Add module for geopotential calculations over vertical prisms Named gravprisms for now but may change as features are added. --- src/potential/CMakeLists.txt | 2 +- src/potential/gravprisms.c | 656 +++++++++++++++++++++++++++++++++++ 2 files changed, 657 insertions(+), 1 deletion(-) create mode 100644 src/potential/gravprisms.c diff --git a/src/potential/CMakeLists.txt b/src/potential/CMakeLists.txt index bdc557d037a..adcb270a00b 100644 --- a/src/potential/CMakeLists.txt +++ b/src/potential/CMakeLists.txt @@ -27,6 +27,6 @@ set (SUPPL_NAME potential) set (SUPPL_HEADERS okbfuns.h talwani.h) set (SUPPL_PROGS_SRCS grdredpol.c gmtgravmag3d.c gmtflexure.c gravfft.c - grdflexure.c grdgravmag3d.c grdseamount.c talwani2d.c talwani3d.c) + gravprisms.c grdflexure.c grdgravmag3d.c grdseamount.c talwani2d.c talwani3d.c) set (SUPPL_LIB_SRCS ${SUPPL_PROGS_SRCS} okbfuns.c) set (SUPPL_EXAMPLE_FILES README.potential) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c new file mode 100644 index 00000000000..f767d4e6ce4 --- /dev/null +++ b/src/potential/gravprisms.c @@ -0,0 +1,656 @@ +/*-------------------------------------------------------------------- + * + * Copyright (c) 1991-2022 by the GMT Team (https://www.generic-mapping-tools.org/team.html) + * See LICENSE.TXT file for copying and redistribution conditions. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; version 3 or any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * Contact info: www.generic-mapping-tools.org + *--------------------------------------------------------------------*/ +/* + * Authord: Paul Wessel and Seung-Sep Kim + * Date: 12-MAR-2022 + * + * + * Calculates gravity due to any number of individual vertical prisms + * that may have their own unique dimensions (or all are constant) and + * individual density contrasts (or all the same). + * Based on methods by + * + * Accelerated with OpenMP; see -x. + */ + +#include "gmt_dev.h" +#include "talwani.h" + +#define THIS_MODULE_CLASSIC_NAME "gravprisms" +#define THIS_MODULE_MODERN_NAME "gravprisms" +#define THIS_MODULE_LIB "potential" +#define THIS_MODULE_PURPOSE "Compute gravity anomalies over 3-D vertical prisms" +#define THIS_MODULE_KEYS " fixed density to override individual prisms */ + bool active; + double rho; + } D; + struct GRAVPRISMS_E { /* -E/ fixed prism dimensions [read from file] */ + bool active; + double dx, dy, dz; + } E; + struct GRAVPRISMS_F { /* -F[f|n[]|v] */ + bool active, lset; + unsigned int mode; + double lat; + } F; + struct GRAVPRISMS_G { /* Output file */ + bool active; + char *file; + } G; + struct GRAVPRISMS_I { /* -Idx[/dy] */ + bool active; + double inc[2]; + } I; + struct GRAVPRISMS_M { /* -Mh|z */ + bool active[2]; /* True if km, else m */ + } M; + struct GRAVPRISMS_N { /* Desired output points */ + bool active; + char *file; + } N; + struct GRAVPRISMS_Z { /* Observation level file or constant */ + bool active; + double level; + unsigned int mode; + char *file; + } Z; +}; + +static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ + struct GRAVPRISMS_CTRL *C = NULL; + + C = gmt_M_memory (GMT, NULL, 1, struct GRAVPRISMS_CTRL); + + /* Initialize values whose defaults are not 0/false/NULL */ + + C->F.lat = 45.0; /* So we compute normal gravity at 45 */ + + return (C); +} + +static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *C) { /* Deallocate control structure */ + if (!C) return; + gmt_M_str_free (C->N.file); + gmt_M_str_free (C->G.file); + gmt_M_str_free (C->Z.file); + gmt_M_free (GMT, C); +} + +static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT_OPTION *options) { + unsigned int k, n_errors = 0; + struct GMT_OPTION *opt = NULL; + struct GMTAPI_CTRL *API = GMT->parent; + + for (opt = options; opt; opt = opt->next) { /* Process all the options given */ + switch (opt->option) { + + case '<': /* Input file(s) */ + if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;; + break; + + case 'A': /* Specify z-axis is positive up [Default is down] */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active); + Ctrl->A.active = true; + break; + case 'D': + n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active); + Ctrl->D.active = true; + Ctrl->D.rho = atof (opt->arg); + break; + case 'E': + n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); + Ctrl->E.active = true; + sscanf (opt->arg, "%lg/%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy, &Ctrl->E.dz); + break; + case 'F': + n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active); + Ctrl->F.active = true; + switch (opt->arg[0]) { + case 'v': Ctrl->F.mode = GRAVPRISMS_VGG; break; + case 'n': Ctrl->F.mode = GRAVPRISMS_GEOID; + if (opt->arg[1]) Ctrl->F.lat = atof (&opt->arg[1]), Ctrl->F.lset = true; + break; + case 'g': Ctrl->F.mode = GRAVPRISMS_FAA; break; + default: + GMT_Report (API, GMT_MSG_WARNING, "Option -F: Unrecognized field %c\n", opt->arg[0]); + n_errors++; + break; + } + break; + case 'G': + n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active); + Ctrl->G.active = true; + Ctrl->G.file = strdup (opt->arg); + break; + case 'I': + n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active); + Ctrl->I.active = true; + n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg); + break; + case 'M': /* Length units */ + k = 0; + while (opt->arg[k]) { + switch (opt->arg[k]) { + case 'h': + n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active[GRAVPRISMS_HOR]); + Ctrl->M.active[GRAVPRISMS_HOR] = true; + break; + case 'z': + n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active[GRAVPRISMS_VER]); + Ctrl->M.active[GRAVPRISMS_VER] = true; + break; + default: + n_errors++; + GMT_Report (API, GMT_MSG_WARNING, "Option -M: Unrecognized modifier %c\n", opt->arg[k]); + break; + } + k++; + } + break; + case 'N': + n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active); + Ctrl->N.active = true; + Ctrl->N.file = strdup (opt->arg); + break; + case 'Z': + n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active); + Ctrl->Z.active = true; + if (!gmt_access (GMT, opt->arg, F_OK)) { /* file exists */ + Ctrl->Z.file = strdup (opt->arg); + Ctrl->Z.mode = 1; + } + else { + Ctrl->Z.mode = 0; + Ctrl->Z.level = atof (opt->arg); + } + break; + default: + n_errors += gmt_default_option_error (GMT, opt); + break; + } + } + if (GMT->common.R.active[RSET]) { + n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[ISET], + "Option -R: Must specify both -R and -I (and optionally -r)\n"); + } + n_errors += gmt_M_check_condition (GMT, (GMT->common.R.active[RSET] && GMT->common.R.active[ISET]) && Ctrl->Z.mode == 1, + "Option -Z: Cannot also specify -R -I\n"); + n_errors += gmt_M_check_condition (GMT, !Ctrl->N.active && !Ctrl->G.active, + "Option -G: Must specify output gridfile name.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && !Ctrl->G.file, + "Option -G: Must specify output gridfile name.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->N.file, + "Option -N: Must specify output gridfile name.\n"); + + return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); +} + +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] [-D] [-Ff|n[]|v] [-G] [%s] " + "[-M[hz]] [-N] [%s] [%s] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", + name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, + GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT); + + if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); + + GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n"); + GMT_Usage (API, 1, "\n"); + GMT_Usage (API, -2, "One or more multiple-segment ASCII data files. If no files are given, standard " + "input is read. Contains (x,y,z) center coordinates of prisms with optional [ dx dy dz] [rho] if not set via -D and ."); + GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); + GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); + GMT_Usage (API, 1, "\n-D"); + GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column."); + GMT_Usage (API, 1, "\n-Ff|n[]|v]"); + GMT_Usage (API, -2, "Specify desired geopotential field component:"); + GMT_Usage (API, 3, "f: Free-air anomalies (mGal) [Default]."); + GMT_Usage (API, 3, "n: Geoid anomalies (meter). Optionally append latitude for evaluation of normal gravity " + "[Default latitude is mid-grid for grid output or mid-latitude if -N is used]."); + GMT_Usage (API, 3, "v: Vertical Gravity Gradient anomalies (Eotvos = 0.1 mGal/km)."); + GMT_Usage (API, 1, "\n-G"); + GMT_Usage (API, -2, "Set name of output file. Grid file name is requires unless -N is used."); + GMT_Option (API, "I"); + GMT_Usage (API, 1, "\n-M[hz]"); + GMT_Usage (API, -2, "Change units used, via one or two directives:"); + GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); + GMT_Usage (API, 3, "z: All z-distances are given in km [meters]."); + GMT_Usage (API, 1, "\n-N"); + GMT_Usage (API, -2, "File with output locations where a calculation is requested. No grid " + "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); + GMT_Option (API, "R,V"); + GMT_Usage (API, 1, "\n-Z"); + GMT_Usage (API, -2, "Set observation level for output locations [0]. " + "Append either a constant or the name of grid file with variable levels. " + "If given a grid then it also defines the output grid."); + GMT_Usage (API, -2, "Note: Cannot use both -Z and -R -I [-r]."); + GMT_Option (API, "bo,d,e"); + GMT_Usage (API, 1, "\n-fg Map units (lon, lat in degree, else in m [but see -Mh])."); + GMT_Option (API, "h,i,o,r,x,."); + return (GMT_MODULE_USAGE); +} + +#define GRAVITATIONAL_CONST 6.674e-11 +#define GRAVITATIONAL_CONST_MGAL 6.674e-6 + +GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { + /* Geoid anomaly from a single prism [Nagy et al, 2000] */ + double n, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; + double R111, R112, R121, R122, R211, R212, R221, R222; + double n111, n112, n121, n122, n211, n212, n221, n222; + + /* Square distances */ + dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; + dy1_sq = dy1 * dy1; dy2_sq = dy2 * dy2; + dz1_sq = dz1 * dz1; dz2_sq = dz2 * dz2; + /* Get radii */ + R111 = sqrt (dx1_sq + dy1_sq + dz1_sq); + R112 = sqrt (dx2_sq + dy1_sq + dz1_sq); + R121 = sqrt (dx1_sq + dy2_sq + dz1_sq); + R122 = sqrt (dx2_sq + dy2_sq + dz1_sq); + R211 = sqrt (dx1_sq + dy1_sq + dz2_sq); + R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); + R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); + R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Evaluate at dz1 */ + n111 = -(0.5 * (dx1_sq * atan ((dy1 * dz1) / (dz1 * R111)) + dy1_sq * atan ((dx1 * dz1) / (dz1 * R111)) + dz1_sq * atan ((dx1 * dy1) / (dz1 * R111))) - dx1 * dz1 * log(R111 + dy1) - dy1 * dz1 * log (R111 + dx1) - dx1 * dy1 * log (R111 + dz1)); + n112 = +(0.5 * (dx2_sq * atan ((dy1 * dz1) / (dz1 * R112)) + dy1_sq * atan ((dx2 * dz1) / (dz1 * R112)) + dz1_sq * atan ((dx2 * dy1) / (dz1 * R112))) - dx2 * dz1 * log(R112 + dy1) - dy1 * dz1 * log (R112 + dx2) - dx2 * dy1 * log (R112 + dz1)); + n121 = +(0.5 * (dx1_sq * atan ((dy2 * dz1) / (dz1 * R121)) + dy2_sq * atan ((dx1 * dz1) / (dz1 * R121)) + dz1_sq * atan ((dx1 * dy2) / (dz1 * R121))) - dx1 * dz1 * log(R121 + dy2) - dy2 * dz1 * log (R121 + dx1) - dx1 * dy2 * log (R121 + dz1)); + n122 = -(0.5 * (dx2_sq * atan ((dy2 * dz1) / (dz1 * R122)) + dy2_sq * atan ((dx2 * dz1) / (dz1 * R122)) + dz1_sq * atan ((dx2 * dy2) / (dz1 * R122))) - dx2 * dz1 * log(R122 + dy2) - dy2 * dz1 * log (R122 + dx2) - dx2 * dy2 * log (R122 + dz1)); + /* Evaluate at dz2 */ + n211 = +(0.5 * (dx1_sq * atan ((dy1 * dz2) / (dz2 * R211)) + dy1_sq * atan ((dx1 * dz2) / (dz2 * R211)) + dz2_sq * atan ((dx1 * dy1) / (dz2 * R211))) - dx1 * dz2 * log (R211 + dy1) - dy1 * dz2 * log (R211 + dx1) - dx1 * dy1 * log (R111 + dz2)); + n212 = -(0.5 * (dx2_sq * atan ((dy1 * dz2) / (dz2 * R212)) + dy1_sq * atan ((dx2 * dz2) / (dz2 * R212)) + dz2_sq * atan ((dx2 * dy1) / (dz2 * R212))) - dx2 * dz2 * log (R212 + dy1) - dy1 * dz2 * log (R212 + dx2) - dx2 * dy1 * log (R112 + dz2)); + n221 = -(0.5 * (dx1_sq * atan ((dy2 * dz2) / (dz2 * R221)) + dy2_sq * atan ((dx1 * dz2) / (dz2 * R221)) + dz2_sq * atan ((dx1 * dy2) / (dz2 * R221))) - dx1 * dz2 * log (R221 + dy2) - dy2 * dz2 * log (R221 + dx1) - dx1 * dy2 * log (R121 + dz2)); + n222 = +(0.5 * (dx2_sq * atan ((dy2 * dz2) / (dz2 * R222)) + dy2_sq * atan ((dx2 * dz2) / (dz2 * R222)) + dz2_sq * atan ((dx2 * dy2) / (dz2 * R222))) - dx2 * dz2 * log (R222 + dy2) - dy2 * dz2 * log (R222 + dx2) - dx2 * dy2 * log (R122 + dz2)); + + n = rho * GRAVITATIONAL_CONST_MGAL * (n111 + n112 + n121 + n122 + n211 + n212 + n221 + n222); + + return (n); +} + +GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { + /* Gravity anomaly from a single prism */ + double g, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; + double R111, R112, R121, R122, R211, R212, R221, R222; + double g111, g112, g121, g122, g211, g212, g221, g222; + + /* Square distances */ + dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; + dy1_sq = dy1 * dy1; dy2_sq = dy2 * dy2; + dz1_sq = dz1 * dz1; dz2_sq = dz2 * dz2; + /* Get radii */ + R111 = sqrt (dx1_sq + dy1_sq + dz1_sq); + R112 = sqrt (dx2_sq + dy1_sq + dz1_sq); + R121 = sqrt (dx1_sq + dy2_sq + dz1_sq); + R122 = sqrt (dx2_sq + dy2_sq + dz1_sq); + R211 = sqrt (dx1_sq + dy1_sq + dz2_sq); + R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); + R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); + R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Evaluate at dz1 */ + g111 = -(dz1 * atan ((dx1 * dy1) / (dz1 * R111)) - dx1 * log(R111 + dy1) - dy1 * log (R111 + dx1)); + g112 = +(dz1 * atan ((dx2 * dy1) / (dz1 * R112)) - dx2 * log(R112 + dy1) - dy1 * log (R112 + dx2)); + g121 = +(dz1 * atan ((dx1 * dy2) / (dz1 * R121)) - dx1 * log(R121 + dy2) - dy2 * log (R121 + dx1)); + g122 = -(dz1 * atan ((dx2 * dy2) / (dz1 * R122)) - dx2 * log(R122 + dy2) - dy2 * log (R122 + dx2)); + /* Evaluate at dz2 */ + g211 = +(dz2 * atan ((dx1 * dy1) / (dz2 * R211)) - dx1 * log (R211 + dy1) - dy1 * log (R211 + dx1)); + g212 = -(dz2 * atan ((dx2 * dy1) / (dz2 * R212)) - dx2 * log (R212 + dy1) - dy1 * log (R212 + dx2)); + g221 = -(dz2 * atan ((dx1 * dy2) / (dz2 * R221)) - dx1 * log (R221 + dy2) - dy2 * log (R221 + dx1)); + g222 = +(dz2 * atan ((dx2 * dy2) / (dz2 * R222)) - dx2 * log (R222 + dy2) - dy2 * log (R222 + dx2)); + + g = rho * GRAVITATIONAL_CONST_MGAL * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); + + return (g); +} + +GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { + /* Vertical gravity gradient from a single prism [Kim & Wessel, 2016] */ + double v, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; + double R111, R112, R121, R122, R211, R212, R221, R222; + double v111, v112, v121, v122, v211, v212, v221, v222; + + /* Square distances */ + dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; + dy1_sq = dy1 * dy1; dy2_sq = dy2 * dy2; + dz1_sq = dz1 * dz1; dz2_sq = dz2 * dz2; + /* Get radii */ + R111 = sqrt (dx1_sq + dy1_sq + dz1_sq); + R112 = sqrt (dx2_sq + dy1_sq + dz1_sq); + R121 = sqrt (dx1_sq + dy2_sq + dz1_sq); + R122 = sqrt (dx2_sq + dy2_sq + dz1_sq); + R211 = sqrt (dx1_sq + dy1_sq + dz2_sq); + R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); + R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); + R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Evaluate at dz1 */ + v111 = -atan ((dx1 * dy1) / (dz1 * R111)); + v112 = +atan ((dx2 * dy1) / (dz1 * R112)); + v121 = +atan ((dx1 * dy2) / (dz1 * R121)); + v122 = -atan ((dx2 * dy2) / (dz1 * R122)); + /* Evaluate at dz2 */ + v211 = +atan ((dx1 * dy1) / (dz2 * R211)); + v212 = -atan ((dx2 * dy1) / (dz2 * R212)); + v221 = -atan ((dx1 * dy2) / (dz2 * R221)); + v222 = +atan ((dx2 * dy2) / (dz2 * R222)); + + v = rho * GRAVITATIONAL_CONST_MGAL * (v111 + v112 + v121 + v122 + v211 + v212 + v221 + v222); + + return (v); +} + +GMT_LOCAL double gravprisms_get_one_g_output (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { + /* (x, y, z) is the observation point */ + double g = 0.0; + gmt_M_unused (unused); + for (uint64_t k = 0; k < n_prisms; k++) + g += gravprism (P[0][k]-x, P[1][k]-x, P[2][k]-y, P[3][k]-y, P[4][k]-z, P[5][k]-z, P[6][k]); + return (g); +} + +GMT_LOCAL double gravprisms_get_one_n_output (double x, double y, double z, uint64_t n_prisms, double **P, double constant) { + /* (x, y, z) is the observation point */ + double n = 0.0; + for (uint64_t k = 0; k < n_prisms; k++) + n += geoidprism (P[0][k]-x, P[1][k]-x, P[2][k]-y, P[3][k]-y, P[4][k]-z, P[5][k]-z, P[6][k]); + return (n * constant * 0.01); /* To get geoid in meter */ +} + +GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { + /* (x, y, z) is the observation point */ + double v = 0.0; + gmt_M_unused (unused); + for (uint64_t k = 0; k < n_prisms; k++) + v += vggprism (P[0][k]-x, P[1][k]-x, P[2][k]-y, P[3][k]-y, P[4][k]-z, P[5][k]-z, P[6][k]); + return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ +} + +#define bailout(code) {gmt_M_free_options (mode); return (code);} +#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} + +EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { + int error = 0; + unsigned int n_expected = 7; + uint64_t tbl, seg, row, col, k, node, n_prisms; + + bool flat_earth = false; + + char *uname[2] = {"meter", "km"}, *kind[3] = {"FAA", "VGG", "GEOID"}, remark[GMT_LEN64] = {""}; + double z_level, rho, dx, dy, dz, lat, G0; + double *prism[7]; + double (*eval) (double, double, double, uint64_t, double **, double); + + struct GRAVPRISMS_CTRL *Ctrl = NULL; + struct GMT_GRID *G = NULL; + struct GMT_DATASET *D = NULL, *P = NULL; + struct GMT_DATASEGMENT *S = NULL; + struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; + struct GMT_OPTION *options = NULL; + struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ + + /*----------------------- Standard module initialization and parsing ----------------------*/ + + if (API == NULL) return (GMT_NOT_A_SESSION); + if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ + options = GMT_Create_Options (API, mode, args); + if (API->error) return (API->error); /* Set or get option list */ + + if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */ + + /* Parse the command-line arguments */ + + if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */ + if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error); + Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */ + if ((error = parse (GMT, Ctrl, options)) != 0) Return (error); + + /*---------------------------- This is the gravprisms main code ----------------------------*/ + + gmt_enable_threads (GMT); /* Set number of active threads, if supported */ + /* Specify input expected columns to be at least 3 */ + if (Ctrl->D.active) n_expected--; /* Not reading density */ + if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz */ + if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { + Return (error); + } + /* Register likely model files unless the caller has already done so */ + if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */ + Return (API->error); + } + if ((P = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) { + Return (API->error); + } + + for (k = 0; k < 7; k++) + prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); + + /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free */ + dx = Ctrl->E.dx; dy = Ctrl->E.dy; dz = Ctrl->E.dz; rho = Ctrl->D.rho; + col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if present */ + for (tbl = k = 0; tbl < P->n_tables; tbl++) { + for (seg = 0; seg < P->table[tbl]->n_segments; seg++) { + S = P->table[tbl]->segment[seg]; + for (row = 0; row < S->n_rows; row++, k++) { + if (!Ctrl->E.active) { + dx = S->data[3][row]; + dy = S->data[4][row]; + dz = S->data[5][row]; + } + if (!Ctrl->D.active) rho = S->data[col][row]; + prism[0][k] = S->data[GMT_X][row] - dx; /* Set the x-coordinates of the x-edges of prism */ + prism[1][k] = S->data[GMT_X][row] + dx; + prism[2][k] = S->data[GMT_Y][row] - dy; /* Set the y-coordinates of the y-edges of prism */ + prism[3][k] = S->data[GMT_Y][row] + dy; + prism[4][k] = S->data[GMT_Z][row] - dz; /* Set the z-coordinates of the z-edges of prism */ + prism[5][k] = S->data[GMT_Z][row] + dz; + prism[6][k] = rho; + } + } + } + n_prisms = k; + if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) { + Return (GMT_MEMORY_ERROR); + } + + if (Ctrl->Z.mode == 1) { /* Got grid with observation levels which also sets output locations; it could also set -fg so do this first */ + if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Z.file, NULL)) == NULL) + Return (API->error); + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); + } + else if (GMT->common.R.active[RSET]) { /* Gave -R -I [-r] and possibly -fg indirectly via geographic coordinates in -R */ + if ((G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, + GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) + Return (API->error); + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); + } + else { /* Got a dataset with output locations via -N */ + gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from the -N file */ + if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) + Return (API->error); + if (D->n_columns < 2) { + GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->N.file, (int)D->n_columns); + Return (GMT_DIM_TOO_SMALL); + } + gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */ + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (D->min[GMT_Y] + D->max[GMT_Y]); + } + + flat_earth = gmt_M_is_geographic (GMT, GMT_IN); /* If true then input is in degrees and we must convert to km later on */ + + if (flat_earth && Ctrl->M.active[GRAVPRISMS_HOR]) { + GMT_Report (API, GMT_MSG_ERROR, "Option -M: Cannot specify both geographic coordinates (degrees) AND -Mh\n"); + Return (GMT_RUNTIME_ERROR); + } + + if (Ctrl->A.active) Ctrl->Z.level = -Ctrl->Z.level; + + /* Read polygon information from multiple segment file */ + GMT_Report (API, GMT_MSG_INFORMATION, "All x/y-values are assumed to be given in %s\n", uname[Ctrl->M.active[GRAVPRISMS_HOR]]); + GMT_Report (API, GMT_MSG_INFORMATION, "All z-values are assumed to be given in %s\n", uname[Ctrl->M.active[GRAVPRISMS_VER]]); + + /* Now we can write (if -V) to the screen the user's polygon model characteristics. */ + + GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms: %" PRIu64 "\n", n_prisms); + GMT_Report (API, GMT_MSG_INFORMATION, "Start calculating %s\n", kind[Ctrl->F.mode]); + + switch (Ctrl->F.mode) { /* Set pointer to chosen evaluation function */ + case GRAVPRISMS_VGG: + eval = &gravprisms_get_one_v_output; + break; + case GRAVPRISMS_GEOID: + eval = &gravprisms_get_one_n_output; + G0 = (Ctrl->F.lset) ? g_normal (Ctrl->F.lat) : g_normal (lat); + G0 = 1.0 / G0; /* Since we will be dividing */ + break; + case GRAVPRISMS_FAA: + eval = &gravprisms_get_one_g_output; + break; + default: + /* Just for Coverity */ + break; + } + + if (Ctrl->N.active) { /* Single loop over specified output locations */ + unsigned int wmode = GMT_ADD_DEFAULT; + double scl = (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) ? (1.0 / METERS_IN_A_KM) : 1.0; /* Perhaps convert to km */ + double out[4]; + struct GMT_RECORD *Rec = gmt_new_record (GMT, out, NULL); + /* Must register Ctrl->G.file first since we are going to writing rec-by-rec */ + if (Ctrl->G.active) { + int out_ID; + if ((out_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_OUT, NULL, Ctrl->G.file)) == GMT_NOTSET) { + Return (API->error); + } + wmode = GMT_ADD_EXISTING; + } + if ((error = GMT_Set_Columns (API, GMT_OUT, 4, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { + Return (error); + } + if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, wmode, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */ + Return (API->error); + } + if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */ + Return (API->error); + } + if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */ + Return (API->error); + } + if (D->n_segments > 1) gmt_set_segmentheader (GMT, GMT_OUT, true); + for (tbl = 0; tbl < D->n_tables; tbl++) { + for (seg = 0; seg < D->table[tbl]->n_segments; seg++) { + int64_t row; + S = D->table[tbl]->segment[seg]; /* Current segment */ + GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, S->header); + gmt_prep_tmp_arrays (GMT, GMT_OUT, S->n_rows, 1); /* Init or reallocate tmp vector */ +#ifdef _OPENMP + /* Spread calculation over selected cores */ +#pragma omp parallel for private(row,z_level) shared(GMT,Ctrl,S,eval,scl,n_prisms,prism,flat_earth, G0) +#endif + /* Separate the calculation from the output in two separate row-loops since cannot do rec-by-rec output + * with OpenMP due to race condiations that would mess up the output order */ + for (row = 0; row < (int64_t)S->n_rows; row++) { /* Calculate attraction at all output locations for this segment */ + z_level = (S->n_columns == 3 && !Ctrl->Z.active) ? S->data[GMT_Z][row] : Ctrl->Z.level; /* Default observation z level unless provided in input file */ + GMT->hidden.mem_coord[GMT_X][row] = eval (S->data[GMT_X][row] * scl, S->data[GMT_Y][row] * scl, z_level, n_prisms, prism, G0); + } + /* This loop is not under OpenMP */ + out[GMT_Z] = Ctrl->Z.level; /* Default observation z level unless provided in input file */ + for (row = 0; row < (int64_t)S->n_rows; row++) { /* Loop over output locations */ + out[GMT_X] = S->data[GMT_X][row]; + out[GMT_Y] = S->data[GMT_Y][row]; + if (S->n_columns == 3 && !Ctrl->Z.active) out[GMT_Z] = S->data[GMT_Z][row]; + out[3] = GMT->hidden.mem_coord[GMT_X][row]; + GMT_Put_Record (API, GMT_WRITE_DATA, Rec); /* Write this to output */ + } + } + } + gmt_M_free (GMT, Rec); + if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */ + Return (API->error); + } + } + else { /* Dealing with a grid */ + openmp_int row, col, n_columns = (openmp_int)G->header->n_columns, n_rows = (openmp_int)G->header->n_rows; /* To shut up compiler warnings */ + double y_obs, *x_obs = gmt_M_memory (GMT, NULL, G->header->n_columns, double); + for (col = 0; col < n_columns; col++) { + x_obs[col] = gmt_M_grd_col_to_x (GMT, col, G->header); + if (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) x_obs[col] /= METERS_IN_A_KM; /* Convert to km */ + } +#ifdef _OPENMP + /* Spread calculation over selected cores */ +#pragma omp parallel for private(row,col,node,y_obs,z_level) shared(API,GMT,Ctrl,G,eval,x_obs,n_prisms,prism,flat_earth, G0) +#endif + for (row = 0; row < n_rows; row++) { /* Do row-by-row and report on progress if -V */ + y_obs = gmt_M_grd_row_to_y (GMT, row, G->header); + if (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) y_obs /= METERS_IN_A_KM; /* Convert to km */ +#ifndef _OPENMP + GMT_Report (API, GMT_MSG_INFORMATION, "Finished row %5d\n", row); +#endif + for (col = 0; col < (openmp_int)G->header->n_columns; col++) { + /* Loop over cols; always save the next level before we update the array at that col */ + node = gmt_M_ijp (G->header, row, col); + z_level = (Ctrl->A.active) ? -G->data[node] : G->data[node]; /* Get observation z level and possibly flip direction */ + G->data[node] = (gmt_grdfloat) eval (x_obs[col], y_obs, z_level, n_prisms, prism, G0); + } + } + gmt_M_free (GMT, x_obs); + GMT_Report (API, GMT_MSG_INFORMATION, "Create %s\n", Ctrl->G.file); + sprintf (remark, "Calculated 3-D %s", kind[Ctrl->F.mode]); + if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, remark, G)) { + Return (API->error); + } + if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G)) { + Return (API->error); + } + if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, G) != GMT_NOERROR) { + Return (API->error); + } + } + + for (k = 0; k < 7; k++) + gmt_M_free (GMT, prism[k]); + + Return (GMT_NOERROR); +} From c86942745664f90e54ace82154366a97ccef4703 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sat, 12 Mar 2022 18:18:14 +0000 Subject: [PATCH 02/27] Update gravprisms.c --- src/potential/gravprisms.c | 226 ++++++++++++++++++++++++++++++------- 1 file changed, 183 insertions(+), 43 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index f767d4e6ce4..07b7d764454 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -58,7 +58,6 @@ enum gravprisms_fields { GRAVPRISMS_VER=1 }; - struct GRAVPRISMS_CTRL { struct GRAVPRISMS_A { /* -A Set positive up */ bool active; @@ -67,6 +66,11 @@ struct GRAVPRISMS_CTRL { bool active; double rho; } D; + struct GRAVPRISMS_C { /* -C//] creates prisms between surfaces */ + bool active; + char *base, *top, *height; + double dz; + } C; struct GRAVPRISMS_E { /* -E/ fixed prism dimensions [read from file] */ bool active; double dx, dy, dz; @@ -80,7 +84,14 @@ struct GRAVPRISMS_CTRL { bool active; char *file; } G; - struct GRAVPRISMS_I { /* -Idx[/dy] */ + struct GRAVPRISMS_H { /* -H//[+d][+p] */ + bool active; + double H, rho_l, rho_h; + double p, densify; + double del_rho; /* Computed as rho_h - rho_l */ + double p1; /* Will be p + 1 to avoid addition in loops */ + } H; +struct GRAVPRISMS_I { /* -Idx[/dy] */ bool active; double inc[2]; } I; @@ -114,13 +125,18 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *C) { /* Deallocate control structure */ if (!C) return; gmt_M_str_free (C->N.file); + gmt_M_str_free (C->C.base); + gmt_M_str_free (C->C.top); + gmt_M_str_free (C->C.height); gmt_M_str_free (C->G.file); gmt_M_str_free (C->Z.file); gmt_M_free (GMT, C); } static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT_OPTION *options) { + int nc = 0; unsigned int k, n_errors = 0; + char *c = NULL, A[GMT_LEN256] = {""}, B[GMT_LEN256] = {""}, C[GMT_LEN256] = {""}; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; @@ -135,6 +151,14 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active); Ctrl->A.active = true; break; + case 'C': + n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active); + Ctrl->C.active = true; + nc = sscanf (opt->arg, "%[^,],%[^,],%[^,]%lg", A, B, C, &Ctrl->C.dz); + if (A[0] != '-') Ctrl->C.base = strdup (A); + if (B[0] != '-') Ctrl->C.top = strdup (B); + if (C[0] != '-') Ctrl->C.height = strdup (C); + break; case 'D': n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active); Ctrl->D.active = true; @@ -165,6 +189,37 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->G.active = true; Ctrl->G.file = strdup (opt->arg); break; + case 'H': /* Reference seamount density parameters */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->H.active); + Ctrl->H.active = true; + if ((c = gmt_first_modifier (GMT, opt->arg, "dp"))) { + unsigned int pos = 0; + char txt[GMT_LEN256] = {""}; + while (gmt_getmodopt (GMT, 'H', c, "dp", &pos, txt, &n_errors) && n_errors == 0) { + switch (txt[0]) { + case 'd': /* Get densify rate over reference water depth H */ + Ctrl->H.densify = atof (&txt[1]); + if (Ctrl->H.densify < 10.0) Ctrl->H.densify *= 1000; /* Gave units of g/cm^3 */ + break; + case 'p': /* Get power coefficient */ + Ctrl->H.p = atof (&txt[1]); + break; + default: + n_errors++; + break; + } + } + c[0] = '\0'; /* Chop off all modifiers so range can be determined */ + } + if (sscanf (opt->arg, "%lg/%lg/%lg", &Ctrl->H.H, &Ctrl->H.rho_l, &Ctrl->H.rho_h) != 3) { + GMT_Report (API, GMT_MSG_ERROR, "Option -H: Unable to parse the three values\n"); + n_errors++; + } + if (Ctrl->H.rho_l < 10.0) Ctrl->H.rho_l *= 1000; /* Gave units of g/cm^3 */ + if (Ctrl->H.rho_h < 10.0) Ctrl->H.rho_h *= 1000; /* Gave units of g/cm^3 */ + Ctrl->H.del_rho = Ctrl->H.rho_h - Ctrl->H.rho_l; + Ctrl->H.p1 = Ctrl->H.p + 1; + break; case 'I': n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active); Ctrl->I.active = true; @@ -224,6 +279,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -G: Must specify output gridfile name.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && !Ctrl->N.file, "Option -N: Must specify output gridfile name.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->H.active, + "Option -H: Cannot be used with -D.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && nc < 4, + "Option -C: Requires when -H is used\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -231,7 +290,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-D] [-Ff|n[]|v] [-G] [%s] " + GMT_Usage (API, 0, "usage: %s [-A] [-C,,,] [-D] [-Ff|n[]|v] [-G] [-H//[+d][+p]] [%s] " "[-M[hz]] [-N] [%s] [%s] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT); @@ -244,6 +303,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) { "input is read. Contains (x,y,z) center coordinates of prisms with optional [ dx dy dz] [rho] if not set via -D and ."); GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); + GMT_Usage (API, 1, "\n-C,,,"); + GMT_Usage (API, -2, "Create prisms from the level to level for full seamount , where arguments are grids or constants. " + "If -H is used then is required for the discretization of rho(r,z). Give = - if is the full height."); GMT_Usage (API, 1, "\n-D"); GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column."); GMT_Usage (API, 1, "\n-Ff|n[]|v]"); @@ -254,6 +316,11 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "v: Vertical Gravity Gradient anomalies (Eotvos = 0.1 mGal/km)."); GMT_Usage (API, 1, "\n-G"); GMT_Usage (API, -2, "Set name of output file. Grid file name is requires unless -N is used."); + GMT_Usage (API, 1, "\n-H//[+d][+p]"); + GMT_Usage (API, -2, "Set parameters for the reference seamount height (in m), flank and deep core densities (kg/m^3 or g/cm^3). " + "Control the density function by these modifiers:"); + GMT_Usage (API, 3, "+d Density increase (kg/m^3 or g/cm^3) due to water pressure over full depth implied by [0]."); + GMT_Usage (API, 3, "+p Exponential coefficient (> 0) for density change with burial depth [1 (linear)]."); GMT_Option (API, "I"); GMT_Usage (API, 1, "\n-M[hz]"); GMT_Usage (API, -2, "Change units used, via one or two directives:"); @@ -408,13 +475,23 @@ GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ } +GMT_LOCAL double gravprisms_rho (struct GRAVPRISMS_CTRL *Ctrl, double z, double h) { + /* Passing in the current seamount's normalized height, h(r) and the current normalized radial point z(r). + * We return the density at this point inside the seamount with reference height 1 */ + double u = (h - z) / Ctrl->H.H; /* Normalized depth into seamount. */ + double q = (Ctrl->H.H - h) / Ctrl->H.H; /* Water depth to flank point */ + + double rho = Ctrl->H.rho_l + Ctrl->H.densify * q + Ctrl->H.del_rho * pow (u, Ctrl->H.p); + return (rho); +} + #define bailout(code) {gmt_M_free_options (mode); return (code);} #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);} EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { int error = 0; unsigned int n_expected = 7; - uint64_t tbl, seg, row, col, k, node, n_prisms; + uint64_t tbl, seg, row, col, k, node, n_prisms = 0; bool flat_earth = false; @@ -450,49 +527,113 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { /*---------------------------- This is the gravprisms main code ----------------------------*/ gmt_enable_threads (GMT); /* Set number of active threads, if supported */ - /* Specify input expected columns to be at least 3 */ - if (Ctrl->D.active) n_expected--; /* Not reading density */ - if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz */ - if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { - Return (error); - } - /* Register likely model files unless the caller has already done so */ - if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */ - Return (API->error); - } - if ((P = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) { - Return (API->error); + + if (Ctrl->C.active) { /* Need to create prisms first */ + struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *G = NULL; + double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid; + size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; + + if (access (Ctrl->C.base, F_OK)) /* No file, just a constant */ + base = atof (Ctrl->C.base); + else if ((B = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.base, NULL)) == NULL) + Return (API->error); + if (access (Ctrl->C.top, F_OK)) /* No file, just a constant */ + top = atof (Ctrl->C.base); + else if ((T = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.top, NULL)) == NULL) + Return (API->error); + if (Ctrl->C.height == NULL) /* No file, so height given by top */ + H = T; + else if ((H = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.height, NULL)) == NULL) + Return (API->error); + G = (B) ? B : T; + dx = 0.5 * G->header->inc[GMT_X]; dy = 0.5 * G->header->inc[GMT_Y]; + for (k = 0; k < 7; k++) + prism[k] = gmt_M_memory (GMT, NULL, n_alloc, double); + + gmt_M_grd_loop (GMT, G, row, col, node) { + z2 = (T) ? top : T->data[node]; + z1 = (B) ? base : B->data[node]; + if (z2 <= z1) continue; /* No layer here */ + z_prev = z1; + do { + z_next = (floor (z_prev / Ctrl->C.dz) + 1.0) * Ctrl->C.dz; /* Presumably next regular z-spacing */ + if (z_next <= z_prev) z_next += Ctrl->C.dz; /* Can happen if z1 is a multiple of dz */ + else if (z_next > z2) z_next = z2; /* At the top, clip to limit */ + z_mid = 0.5 * (z_prev + z_next); + rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); + if (n_prisms == n_alloc) { /* Need more memory */ + n_alloc <<= 1; /* Double it */ + for (k = 0; k < 7; k++) + prism[k] = gmt_M_memory (GMT, prism[k], n_alloc, double); + } + prism[0][n_prisms] = G->x[col] - dx; prism[1][n_prisms] = G->x[col] + dx; + prism[2][n_prisms] = G->y[row] - dy; prism[3][n_prisms] = G->y[row] + dy; + prism[4][n_prisms] = z_prev; prism[5][n_prisms] = z_next; + prism[6][n_prisms] = rho; + n_prisms++; + z_prev = z_next; + } while (z_prev < z2); /* Until we run out of this stack */ + } + /* Finalize allocation */ + for (k = 0; k < 7; k++) + prism[k] = gmt_M_memory (GMT, prism[k], n_prisms, double); + GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms constructed: %" PRIu64 "\n", n_prisms); + if (B && GMT_Destroy_Data (API, &B) != GMT_NOERROR) { + Return (GMT_MEMORY_ERROR); + } + if (T && GMT_Destroy_Data (API, &T) != GMT_NOERROR) { + Return (GMT_MEMORY_ERROR); + } + if (H && GMT_Destroy_Data (API, &H) != GMT_NOERROR) { + Return (GMT_MEMORY_ERROR); + } } + else { /* Read from input */ + /* Specify input expected columns to be at least 3 */ + if (Ctrl->D.active) n_expected--; /* Not reading density */ + if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz */ + if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { + Return (error); + } + /* Register likely model files unless the caller has already done so */ + if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */ + Return (API->error); + } + if ((P = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) { + Return (API->error); + } - for (k = 0; k < 7; k++) - prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); - - /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free */ - dx = Ctrl->E.dx; dy = Ctrl->E.dy; dz = Ctrl->E.dz; rho = Ctrl->D.rho; - col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if present */ - for (tbl = k = 0; tbl < P->n_tables; tbl++) { - for (seg = 0; seg < P->table[tbl]->n_segments; seg++) { - S = P->table[tbl]->segment[seg]; - for (row = 0; row < S->n_rows; row++, k++) { - if (!Ctrl->E.active) { - dx = S->data[3][row]; - dy = S->data[4][row]; - dz = S->data[5][row]; + for (k = 0; k < 7; k++) + prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); + + /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free */ + dx = Ctrl->E.dx; dy = Ctrl->E.dy; dz = Ctrl->E.dz; rho = Ctrl->D.rho; + col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if present */ + for (tbl = k = 0; tbl < P->n_tables; tbl++) { + for (seg = 0; seg < P->table[tbl]->n_segments; seg++) { + S = P->table[tbl]->segment[seg]; + for (row = 0; row < S->n_rows; row++, k++) { + if (!Ctrl->E.active) { + dx = S->data[3][row]; + dy = S->data[4][row]; + dz = S->data[5][row]; + } + if (!Ctrl->D.active) rho = S->data[col][row]; + prism[0][k] = S->data[GMT_X][row] - dx; /* Set the x-coordinates of the x-edges of prism */ + prism[1][k] = S->data[GMT_X][row] + dx; + prism[2][k] = S->data[GMT_Y][row] - dy; /* Set the y-coordinates of the y-edges of prism */ + prism[3][k] = S->data[GMT_Y][row] + dy; + prism[4][k] = S->data[GMT_Z][row] - dz; /* Set the z-coordinates of the z-edges of prism */ + prism[5][k] = S->data[GMT_Z][row] + dz; + prism[6][k] = rho; } - if (!Ctrl->D.active) rho = S->data[col][row]; - prism[0][k] = S->data[GMT_X][row] - dx; /* Set the x-coordinates of the x-edges of prism */ - prism[1][k] = S->data[GMT_X][row] + dx; - prism[2][k] = S->data[GMT_Y][row] - dy; /* Set the y-coordinates of the y-edges of prism */ - prism[3][k] = S->data[GMT_Y][row] + dy; - prism[4][k] = S->data[GMT_Z][row] - dz; /* Set the z-coordinates of the z-edges of prism */ - prism[5][k] = S->data[GMT_Z][row] + dz; - prism[6][k] = rho; } } - } - n_prisms = k; - if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) { - Return (GMT_MEMORY_ERROR); + n_prisms = k; + if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) { + Return (GMT_MEMORY_ERROR); + } + GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms read: %" PRIu64 "\n", n_prisms); } if (Ctrl->Z.mode == 1) { /* Got grid with observation levels which also sets output locations; it could also set -fg so do this first */ @@ -533,7 +674,6 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { /* Now we can write (if -V) to the screen the user's polygon model characteristics. */ - GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms: %" PRIu64 "\n", n_prisms); GMT_Report (API, GMT_MSG_INFORMATION, "Start calculating %s\n", kind[Ctrl->F.mode]); switch (Ctrl->F.mode) { /* Set pointer to chosen evaluation function */ From 786fc9498e86a3fc44de029d4205c966955a6220 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sat, 12 Mar 2022 22:07:07 +0000 Subject: [PATCH 03/27] Update gravprisms.c --- src/potential/gravprisms.c | 142 +++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 60 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 07b7d764454..508f379134f 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -151,25 +151,29 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active); Ctrl->A.active = true; break; - case 'C': + case 'C': /* Create prisms from layer between two surfaces */ n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active); Ctrl->C.active = true; - nc = sscanf (opt->arg, "%[^,],%[^,],%[^,]%lg", A, B, C, &Ctrl->C.dz); + nc = sscanf (opt->arg, "%[^,],%[^,],%[^,],%lg", A, B, C, &Ctrl->C.dz); if (A[0] != '-') Ctrl->C.base = strdup (A); if (B[0] != '-') Ctrl->C.top = strdup (B); if (C[0] != '-') Ctrl->C.height = strdup (C); + if (nc < 3) { + GMT_Report (API, GMT_MSG_WARNING, "Option -C: Not enough arguments supplied (only %d detected but 3 or 4 are needed\n", nc); + n_errors++; + } break; case 'D': n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active); Ctrl->D.active = true; Ctrl->D.rho = atof (opt->arg); break; - case 'E': + case 'E': /* Set fixed prism parameters instead of reading from prismfile */ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); Ctrl->E.active = true; sscanf (opt->arg, "%lg/%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy, &Ctrl->E.dz); break; - case 'F': + case 'F': /* Seldct geopotential type */ n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active); Ctrl->F.active = true; switch (opt->arg[0]) { @@ -184,12 +188,12 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT break; } break; - case 'G': + case 'G': /* Output file given */ n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active); Ctrl->G.active = true; Ctrl->G.file = strdup (opt->arg); break; - case 'H': /* Reference seamount density parameters */ + case 'H': /* Reference seamount density parameters for rho(z) */ n_errors += gmt_M_repeated_module_option (API, Ctrl->H.active); Ctrl->H.active = true; if ((c = gmt_first_modifier (GMT, opt->arg, "dp"))) { @@ -220,12 +224,12 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->H.del_rho = Ctrl->H.rho_h - Ctrl->H.rho_l; Ctrl->H.p1 = Ctrl->H.p + 1; break; - case 'I': + case 'I': /* Grid iincrements */ n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active); Ctrl->I.active = true; n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg); break; - case 'M': /* Length units */ + case 'M': /* Horizontal and vertical length units */ k = 0; while (opt->arg[k]) { switch (opt->arg[k]) { @@ -245,19 +249,19 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT k++; } break; - case 'N': + case 'N': /* Got profile for output locations */ n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active); Ctrl->N.active = true; Ctrl->N.file = strdup (opt->arg); break; - case 'Z': + case 'Z': /* Observation level(s) */ n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active); Ctrl->Z.active = true; - if (!gmt_access (GMT, opt->arg, F_OK)) { /* file exists */ + if (!gmt_access (GMT, opt->arg, F_OK)) { /* File with z-levels exists */ Ctrl->Z.file = strdup (opt->arg); Ctrl->Z.mode = 1; } - else { + else { /* Got a constant z-level */ Ctrl->Z.mode = 0; Ctrl->Z.level = atof (opt->arg); } @@ -267,12 +271,15 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT break; } } + if (GMT->common.R.active[RSET]) { n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[ISET], "Option -R: Must specify both -R and -I (and optionally -r)\n"); } n_errors += gmt_M_check_condition (GMT, (GMT->common.R.active[RSET] && GMT->common.R.active[ISET]) && Ctrl->Z.mode == 1, - "Option -Z: Cannot also specify -R -I\n"); + "Option -Z: Cannot also specify -R -I if a level grid has been supplied\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->Z.mode == 1, + "Option -Z: Cannot also specify -N if a level grid has been supplied\n"); n_errors += gmt_M_check_condition (GMT, !Ctrl->N.active && !Ctrl->G.active, "Option -G: Must specify output gridfile name.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && !Ctrl->G.file, @@ -281,8 +288,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -N: Must specify output gridfile name.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->H.active, "Option -H: Cannot be used with -D.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->D.active && !Ctrl->H.active, + "option -C: Need to select either -D or -H.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && nc < 4, - "Option -C: Requires when -H is used\n"); + "Option -C: Requires when -H is selected\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -305,9 +314,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); GMT_Usage (API, 1, "\n-C,,,"); GMT_Usage (API, -2, "Create prisms from the level to level for full seamount , where arguments are grids or constants. " - "If -H is used then is required for the discretization of rho(r,z). Give = - if is the full height."); + "If -H is used then is required for the discretization of rho(r,z). Give = - if is the full height. No will be read."); GMT_Usage (API, 1, "\n-D"); - GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column."); + GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column or computes it via -H]."); GMT_Usage (API, 1, "\n-Ff|n[]|v]"); GMT_Usage (API, -2, "Specify desired geopotential field component:"); GMT_Usage (API, 3, "f: Free-air anomalies (mGal) [Default]."); @@ -323,16 +332,16 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "+p Exponential coefficient (> 0) for density change with burial depth [1 (linear)]."); GMT_Option (API, "I"); GMT_Usage (API, 1, "\n-M[hz]"); - GMT_Usage (API, -2, "Change units used, via one or two directives:"); + GMT_Usage (API, -2, "Change distance units used, via one or two directives:"); GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); GMT_Usage (API, 3, "z: All z-distances are given in km [meters]."); GMT_Usage (API, 1, "\n-N"); - GMT_Usage (API, -2, "File with output locations where a calculation is requested. No grid " + GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally z may be read as well. No grid " "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); GMT_Option (API, "R,V"); GMT_Usage (API, 1, "\n-Z"); GMT_Usage (API, -2, "Set observation level for output locations [0]. " - "Append either a constant or the name of grid file with variable levels. " + "Append either a constant or the name of a grid file with variable levels. " "If given a grid then it also defines the output grid."); GMT_Usage (API, -2, "Note: Cannot use both -Z and -R -I [-r]."); GMT_Option (API, "bo,d,e"); @@ -490,14 +499,12 @@ GMT_LOCAL double gravprisms_rho (struct GRAVPRISMS_CTRL *Ctrl, double z, double EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { int error = 0; - unsigned int n_expected = 7; uint64_t tbl, seg, row, col, k, node, n_prisms = 0; bool flat_earth = false; char *uname[2] = {"meter", "km"}, *kind[3] = {"FAA", "VGG", "GEOID"}, remark[GMT_LEN64] = {""}; - double z_level, rho, dx, dy, dz, lat, G0; - double *prism[7]; + double z_level, rho, dx, dy, dz, lat, G0, *prism[7]; double (*eval) (double, double, double, uint64_t, double **, double); struct GRAVPRISMS_CTRL *Ctrl = NULL; @@ -528,7 +535,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_enable_threads (GMT); /* Set number of active threads, if supported */ - if (Ctrl->C.active) { /* Need to create prisms first */ + if (Ctrl->C.active) { /* Need to create prisms from two surfaces first */ struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *G = NULL; double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid; size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; @@ -547,21 +554,27 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (API->error); G = (B) ? B : T; dx = 0.5 * G->header->inc[GMT_X]; dy = 0.5 * G->header->inc[GMT_Y]; + rho = Ctrl->D.rho; /* Only helpful if -D is set */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, NULL, n_alloc, double); gmt_M_grd_loop (GMT, G, row, col, node) { - z2 = (T) ? top : T->data[node]; - z1 = (B) ? base : B->data[node]; - if (z2 <= z1) continue; /* No layer here */ - z_prev = z1; - do { - z_next = (floor (z_prev / Ctrl->C.dz) + 1.0) * Ctrl->C.dz; /* Presumably next regular z-spacing */ - if (z_next <= z_prev) z_next += Ctrl->C.dz; /* Can happen if z1 is a multiple of dz */ - else if (z_next > z2) z_next = z2; /* At the top, clip to limit */ - z_mid = 0.5 * (z_prev + z_next); - rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); - if (n_prisms == n_alloc) { /* Need more memory */ + z2 = (T) ? T->data[node] : top; + z1 = (B) ? B->data[node] : base; + if (z2 <= z1) continue; /* No layer thickness detected */ + z_prev = z1; /* We start exactly at z = z1 */ + do { /* Will at least be one prism */ + if (Ctrl->H.active) { /* Variable density means stacked prisms, most of thickness dz, at same point */ + /* First and last prisms in the stack are adjusted to have the height needed to match the surface boundaries */ + z_next = (floor (z_prev / Ctrl->C.dz) + 1.0) * Ctrl->C.dz; /* Presumably next regular z-spacing */ + if (z_next <= z_prev) z_next += Ctrl->C.dz; /* Can happen if z1 is a multiple of dz */ + else if (z_next > z2) z_next = z2; /* At the top, clip to limit */ + z_mid = 0.5 * (z_prev + z_next); /* Middle of prism - used to look up density */ + rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); + } + else /* Constant density rho (set above), just need a single prism per location */ + z_next = z2; + if (n_prisms == n_alloc) { /* Need to allocate more memory for the prisms */ n_alloc <<= 1; /* Double it */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_alloc, double); @@ -571,7 +584,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { prism[4][n_prisms] = z_prev; prism[5][n_prisms] = z_next; prism[6][n_prisms] = rho; n_prisms++; - z_prev = z_next; + z_prev = z_next; /* The the top of this prism be the bottom of the next */ } while (z_prev < z2); /* Until we run out of this stack */ } /* Finalize allocation */ @@ -588,10 +601,11 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (GMT_MEMORY_ERROR); } } - else { /* Read from input */ + else { /* Read prisms from stdin or input file(s) */ + unsigned int n_expected = 7; /* Max number of columns */ /* Specify input expected columns to be at least 3 */ - if (Ctrl->D.active) n_expected--; /* Not reading density */ - if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz */ + if (Ctrl->D.active) n_expected--; /* Not reading density from records */ + if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz from records */ if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { Return (error); } @@ -599,53 +613,61 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default input sources, unless already set */ Return (API->error); } + /* Read the entire prism file(s) */ if ((P = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) { Return (API->error); } - + /* To avoid the need to loop over tables and segments, we extrac the data into a separate array and comptute prism edges instead */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); - /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free */ - dx = Ctrl->E.dx; dy = Ctrl->E.dy; dz = Ctrl->E.dz; rho = Ctrl->D.rho; - col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if present */ + /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free the dataset */ + dx = 0.5 * Ctrl->E.dx; dy = 0.5 * Ctrl->E.dy; dz = 0.5 * Ctrl->E.dz; /* Distances from prism center to respective edges */ + rho = Ctrl->D.rho; + col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if -D is not set */ for (tbl = k = 0; tbl < P->n_tables; tbl++) { for (seg = 0; seg < P->table[tbl]->n_segments; seg++) { S = P->table[tbl]->segment[seg]; for (row = 0; row < S->n_rows; row++, k++) { - if (!Ctrl->E.active) { - dx = S->data[3][row]; - dy = S->data[4][row]; - dz = S->data[5][row]; + if (!Ctrl->E.active) { /* Update the half-dimensions of the sides */ + dx = 0.5 * S->data[3][row]; + dy = 0.5 * S->data[4][row]; + dz = 0.5 * S->data[5][row]; } if (!Ctrl->D.active) rho = S->data[col][row]; - prism[0][k] = S->data[GMT_X][row] - dx; /* Set the x-coordinates of the x-edges of prism */ + prism[0][k] = S->data[GMT_X][row] - dx; /* Store the x-coordinates of the x-edges of prism */ prism[1][k] = S->data[GMT_X][row] + dx; - prism[2][k] = S->data[GMT_Y][row] - dy; /* Set the y-coordinates of the y-edges of prism */ + prism[2][k] = S->data[GMT_Y][row] - dy; /* Store the y-coordinates of the y-edges of prism */ prism[3][k] = S->data[GMT_Y][row] + dy; - prism[4][k] = S->data[GMT_Z][row] - dz; /* Set the z-coordinates of the z-edges of prism */ + prism[4][k] = S->data[GMT_Z][row] - dz; /* Store the z-coordinates of the z-edges of prism */ prism[5][k] = S->data[GMT_Z][row] + dz; - prism[6][k] = rho; + prism[6][k] = rho; /* Store the fixed or variable density of prism */ } } } n_prisms = k; + /* Free the data set */ if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) { Return (GMT_MEMORY_ERROR); } GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms read: %" PRIu64 "\n", n_prisms); } + if (n_prisms == 0) { + GMT_Report (API, GMT_MSG_WARNING, "No prisms found - exiting\n"); + Return (GMT_NOERROR); + } + if (Ctrl->Z.mode == 1) { /* Got grid with observation levels which also sets output locations; it could also set -fg so do this first */ if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Z.file, NULL)) == NULL) Return (API->error); - if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); /* Mid-latitude needed if geoid is selected */ } else if (GMT->common.R.active[RSET]) { /* Gave -R -I [-r] and possibly -fg indirectly via geographic coordinates in -R */ if ((G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error); - if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); /* Mid-latitude needed if geoid is selected */ } else { /* Got a dataset with output locations via -N */ gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from the -N file */ @@ -656,7 +678,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (GMT_DIM_TOO_SMALL); } gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */ - if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (D->min[GMT_Y] + D->max[GMT_Y]); + if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (D->min[GMT_Y] + D->max[GMT_Y]); /* Mid-latitude needed if geoid is selected */ } flat_earth = gmt_M_is_geographic (GMT, GMT_IN); /* If true then input is in degrees and we must convert to km later on */ @@ -676,14 +698,14 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_INFORMATION, "Start calculating %s\n", kind[Ctrl->F.mode]); - switch (Ctrl->F.mode) { /* Set pointer to chosen evaluation function */ + switch (Ctrl->F.mode) { /* Set pointer to chosen geopotential evaluation function */ case GRAVPRISMS_VGG: eval = &gravprisms_get_one_v_output; break; case GRAVPRISMS_GEOID: eval = &gravprisms_get_one_n_output; G0 = (Ctrl->F.lset) ? g_normal (Ctrl->F.lat) : g_normal (lat); - G0 = 1.0 / G0; /* Since we will be dividing */ + G0 = 1.0 / G0; /* To avoid dividing in the loop */ break; case GRAVPRISMS_FAA: eval = &gravprisms_get_one_g_output; @@ -727,7 +749,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_prep_tmp_arrays (GMT, GMT_OUT, S->n_rows, 1); /* Init or reallocate tmp vector */ #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,z_level) shared(GMT,Ctrl,S,eval,scl,n_prisms,prism,flat_earth, G0) +#pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl,n_prisms,prism,flat_earth,G0) #endif /* Separate the calculation from the output in two separate row-loops since cannot do rec-by-rec output * with OpenMP due to race condiations that would mess up the output order */ @@ -760,7 +782,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { } #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,col,node,y_obs,z_level) shared(API,GMT,Ctrl,G,eval,x_obs,n_prisms,prism,flat_earth, G0) +#pragma omp parallel for private(row,y_obs,col,node,z_level) shared(n_rows,GMT,G,flat_earth,Ctrl,n_columns,eval,x_obs,n_prisms,prism,G0) #endif for (row = 0; row < n_rows; row++) { /* Do row-by-row and report on progress if -V */ y_obs = gmt_M_grd_row_to_y (GMT, row, G->header); @@ -768,10 +790,10 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { #ifndef _OPENMP GMT_Report (API, GMT_MSG_INFORMATION, "Finished row %5d\n", row); #endif - for (col = 0; col < (openmp_int)G->header->n_columns; col++) { + for (col = 0; col < n_columns; col++) { /* Loop over cols; always save the next level before we update the array at that col */ node = gmt_M_ijp (G->header, row, col); - z_level = (Ctrl->A.active) ? -G->data[node] : G->data[node]; /* Get observation z level and possibly flip direction */ + z_level = (Ctrl->Z.mode == 1) ? G->data[node] : Ctrl->Z.level; /* Default observation z level unless provided in input grid */ G->data[node] = (gmt_grdfloat) eval (x_obs[col], y_obs, z_level, n_prisms, prism, G0); } } @@ -788,7 +810,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (API->error); } } - + /* Free the prism information */ for (k = 0; k < 7; k++) gmt_M_free (GMT, prism[k]); From 437ff50e983d390ce89d234ba94ca8ecbe54b07b Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 10:20:27 +0000 Subject: [PATCH 04/27] Add documentation --- doc/rst/source/modules-classic.rst | 4 + doc/rst/source/modules.rst | 4 + .../module_supplements_purpose.rst_ | 2 + .../supplements/potential/gravprisms.rst | 239 ++++++++++++++++++ src/potential/gravprisms.c | 41 ++- 5 files changed, 277 insertions(+), 13 deletions(-) create mode 100644 doc/rst/source/supplements/potential/gravprisms.rst diff --git a/doc/rst/source/modules-classic.rst b/doc/rst/source/modules-classic.rst index 5ce60f54349..9faecb80978 100644 --- a/doc/rst/source/modules-classic.rst +++ b/doc/rst/source/modules-classic.rst @@ -139,6 +139,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/potential/gmtflexure supplements/potential/gmtgravmag3d supplements/potential/gravfft + supplements/potential/gravprisms supplements/potential/grdflexure supplements/potential/grdgravmag3d supplements/potential/grdredpol @@ -299,6 +300,7 @@ Supplemental Modules - :doc:`/supplements/potential/gmtflexure` - :doc:`/supplements/potential/gmtgravmag3d` - :doc:`/supplements/potential/gravfft` + - :doc:`/supplements/potential/gravprisms` - :doc:`/supplements/potential/grdflexure` - :doc:`/supplements/potential/grdgravmag3d` - :doc:`/supplements/potential/grdredpol` @@ -646,6 +648,8 @@ potential +--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/gravfft` | |gravfft_purpose| | +--------------------------------------------+--------------------------+ +| :doc:`/supplements/potential/gravprisms` | |gravprisms_purpose| | ++--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/grdflexure` | |grdflexure_purpose| | +--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/grdgravmag3d` | |grdgravmag3d_purpose| | diff --git a/doc/rst/source/modules.rst b/doc/rst/source/modules.rst index 922489407b6..9156ba2eac3 100644 --- a/doc/rst/source/modules.rst +++ b/doc/rst/source/modules.rst @@ -144,6 +144,7 @@ All modules are requested via a call to the :doc:`gmt` program. supplements/potential/gmtflexure supplements/potential/gmtgravmag3d supplements/potential/gravfft + supplements/potential/gravprisms supplements/potential/grdflexure supplements/potential/grdgravmag3d supplements/potential/grdredpol @@ -311,6 +312,7 @@ Supplemental Modules - :doc:`/supplements/potential/gmtflexure` - :doc:`/supplements/potential/gmtgravmag3d` - :doc:`/supplements/potential/gravfft` + - :doc:`/supplements/potential/gravprisms` - :doc:`/supplements/potential/grdflexure` - :doc:`/supplements/potential/grdgravmag3d` - :doc:`/supplements/potential/grdredpol` @@ -675,6 +677,8 @@ potential +--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/gravfft` | |gravfft_purpose| | +--------------------------------------------+--------------------------+ +| :doc:`/supplements/potential/gravprisms` | |gravprisms_purpose| | ++--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/grdflexure` | |grdflexure_purpose| | +--------------------------------------------+--------------------------+ | :doc:`/supplements/potential/grdgravmag3d` | |grdgravmag3d_purpose| | diff --git a/doc/rst/source/supplements/module_supplements_purpose.rst_ b/doc/rst/source/supplements/module_supplements_purpose.rst_ index 4a26acb9ca4..8798f91aacd 100644 --- a/doc/rst/source/supplements/module_supplements_purpose.rst_ +++ b/doc/rst/source/supplements/module_supplements_purpose.rst_ @@ -36,6 +36,8 @@ .. |gravfft_purpose| replace:: Spectral calculations of gravity, isostasy, admittance, and coherence for grids +.. |gravprisms_purpose| replace:: Compute gravity anomalies over 3-D vertical prisms + .. |grdflexure_purpose| replace:: Compute flexural deformation of 3-D surfaces for various rheologies .. |grdgravmag3d_purpose| replace:: Computes the gravity effect of one (or two) grids by the method of Okabe diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst new file mode 100644 index 00000000000..8472ca7c60e --- /dev/null +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -0,0 +1,239 @@ +.. index:: ! gravprisms +.. include:: ../module_supplements_purpose.rst_ + +********** +gravprisms +********** + +|gravprisms_purpose| + +Synopsis +-------- + +.. include:: ../../common_SYN_OPTs.rst_ + +**gmt gravprisms** [ *table* ] +[ |-A| ] +[ |-C|\ [*base*\ /*top*]\ *height*\ [**+z**\ *dz*] ] +[ |-D|\ *density* ] ] +[ |-E|\ *dx*\ /*dy*\ /*dz* ] +[ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] +[ |-G|\ *outfile* ] +[ |-H|\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] ] +[ |SYN_OPT-I| ] +[ |-M|\ [**h**]\ [**v**] ] +[ |-N|\ *trackfile* ] +[ |SYN_OPT-R| ] +[ |-Z|\ *level*\|\ *obsgrid* ] +[ |SYN_OPT-V| ] +[ |SYN_OPT-bo| ] +[ |SYN_OPT-d| ] +[ |SYN_OPT-e| ] +[ |SYN_OPT-f| ] +[ |SYN_OPT-i| ] +[ |SYN_OPT-o| ] +[ |SYN_OPT-r| ] +[ |SYN_OPT-x| ] +[ |SYN_OPT--| ] + +|No-spaces| + +Description +----------- + +**gravprisms** will compute the geopotential field over vertically oriented, rectangular prisms. +We either read the multi-segment *table* from file (or standard input), which may contain up to +7 colums: The first three are the center *x, y, z* coordinates of the prism, while the next +three are the dimensions *dx dy dz* of each prism (see **-E** if all prisms have the same dimesions. +Last column may contain individual prism densities (may be overridden by a fixed density contrast +given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate the space between +two surfaces (one of which may be a constant). If a variable density model (**-H**) is selected +then each vertical prism will be broken into constant-density, stacked sub-prisms using a prescribed +vertical increment *dz*. +We can compute anomalies on an equidistant grid (by specifying a new grid with +**-R** and **-I** or provide an observation grid with desired elevations) or at arbitrary +output points specified via **-N**. Choose between free-air anomalies, vertical +gravity gradient anomalies, or geoid anomalies. Options are available to control +axes units and direction. + + +Required Arguments +------------------ + +*table* + The file describing the prisms with record format *x y z* [ *dx dy dz* ] [ *rho* ], + where the optional items are controlled by options **-E** and **-D**, respectively. + Any density contrast can be given in kg/m^3 of g/cm^3. + +.. _-I: + +.. include:: ../../explain_-I.rst_ + +.. |Add_-R| replace:: |Add_-R_links| +.. include:: ../../explain_-R.rst_ + :start-after: **Syntax** + :end-before: **Description** + +Optional Arguments +------------------ + +.. _-A: + +**-A** + The *z*-axis should be positive upwards [Default is down]. + +.. _-C: + +**-C**\ [*base*\ /*top*]\ *height*\ [**+z**\ *dz*] ] + Create prisms for the layer between the two surfaces *base* and *top*, with *height* describing the full height of the + feature. Either *base* or *top* may be a constant rather than grid. If *top* equals *height* then + give *height* as -. If only *height* is given then we assume we will approximate the entire feature + from *base* = 0 to *height*. If **-H** is used to compute variable density contrasts then we must + split each prism into a stack of sub-prisms with individual densities; thus *dz* needs to be set for the heights + of these sub-prisms (the first and last sub-prisms in the stack may have their heights adjusted to match the + limits of the surfaces). Without **-H** we only create a single uniform-density prism. + +.. _-D: + +**-D**\ *density* + Sets a fixed density contrast that overrides any individual prism settings in the prisms file, in kg/m^3 of g/cm^3. + +.. _-E: + +**-E**\ *dx*\ /*dy*\ /*dz* + If all prisms in *table* have constant dimensions then they can be set here. In that case *table* + only contains the centers of each prism (and optionally *density*; see **-D**). + +.. _-F: + +**-F**\ **f**\|\ **n**\|\ **v** + Specify desired gravitational field component. Choose between **f** (free-air anomaly) [Default], + **n** (geoid; optionally append average latitude for normal gravity reference value [Default is + mid-grid (or mid-profile if **-N**)]) or **v** (vertical gravity gradient). + +.. _-G: + +**-G**\ *outfile* + Specify the name of the output data (for grids, see :ref:`Grid File Formats + `). Required when an equidistant grid is implied for output. + If **-N** is used then output is written to standard output unless **-G** specifies an output file. + +.. _-H: + +**-H**\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] + Set reference seamount parameters for an *ad-hoc* variable radial density function with depth. Give + the low and high seamount densities in kg/m^3 or g/cm^3 and the fixed reference height *H* in meters. + Use modifers **+d** and **+p** to change the water-pressure-driven flank density increase + over the full reference height [0] and the variable density profile exponent *power* [1, i.e., a linear change]. + See :doc:`grdseamount ` for more details. + +.. _-M: + +**-M**\ [**h**]\ [**v**] + Sets distance units used. Append **h** to indicate that both horizontal distances are in km [m], + and append **z** to indicate vertical distances are in km [m]. + +.. _-N: + +**-N**\ *trackfile* + Specifies individual (x, y[, z]) locations where we wish to compute the predicted value. When this option + is used there are no grids and the output data records are written to standard output (see **-bo** for binary output). + If *trackfile* has 3 columns we take the *z* value as our observation level; this level may be overridden via **-Z**. + +.. |Add_-V| replace:: |Add_-V_links| +.. include:: /explain_-V.rst_ + :start-after: **Syntax** + :end-before: **Description** + +.. _-Z: + +**-Z**\ *level*\|\ *obsgrid* + Set observation level, either as a constant or variable by giving the name of a grid with observation + levels. If the latter is used then this grid determines the output grid region as well [0]. + +.. |Add_-bo| replace:: [Default is 2 output columns]. +.. include:: ../../explain_-bo.rst_ + +.. |Add_-d| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-d.rst_ + +.. |Add_-e| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_-e.rst_ + +|SYN_OPT-f| + Geographic grids (i.e., dimensions of longitude, latitude) will be converted to + km via a "Flat Earth" approximation using the current ellipsoidal parameters. + +.. |Add_-h| replace:: Not used with binary data. +.. include:: ../../explain_-h.rst_ + +.. include:: ../../explain_-icols.rst_ + +.. include:: ../../explain_-ocols.rst_ + +.. |Add_nodereg| unicode:: 0x20 .. just an invisible code +.. include:: ../../explain_nodereg.rst_ + +.. include:: ../../explain_core.rst_ + +.. include:: ../../explain_colon.rst_ + +.. include:: ../../explain_help.rst_ + +.. include:: ../../explain_distunits.rst_ + + +Examples +-------- + +To compute the free-air anomalies on a grid over a set of prisms given in prisms.txt, +using 1700 kg/m^3 as the fixed density contrast, with +horizontal distances in km and vertical distances in meters, try + +:: + + gmt gravprisms -R-200/200/-200/200 -I2 -Mh -G3dgrav.nc prisms.txt -D1700 -Ff + +To obtain the vertical gravity gradient anomaly along the track in crossing.txt +for the same model, try + +:: + + gmt gravprisms -Ncrossing.txt -Mh prisms.txt -D1700 -Fv > vgg_crossing.txt + + +Finally, the geoid anomaly along the same track in crossing.txt +for the same model (at 30S) is written to n_crossing.txt by + +:: + + gmt gravprisms -Ncrossing.txt -Mh prisms.txt -D1700 -Fn-30 -Gn_crossing.txt + + +References +---------- + + +Grant, F. S. and West, G. F., 1965, *Interpretation theory in applied geophysics*, +583 pp., McGraw-Hill. + +Kim, S.-S., and P. Wessel, 2016, New analytic solutions for modeling vertical +gravity gradient anomalies, *Geochem. Geophys. Geosyst., 17*, +`https://dx.doi.org/10.1002/2016GC006263 `_. + +Nagy D., Papp G., Benedek J., 2000, The gravitational potential and its derivatives +for the prism, *J. Geod., 74*, 552–560, +`https://dx.doi.org/10.1007/s001900000116 `_. + +See Also +-------- + +:doc:`gmt.conf `, +:doc:`gmt `, +:doc:`grdmath `, +:doc:`gravfft `, +:doc:`gmtgravmag3d `, +:doc:`grdgravmag3d `, +:doc:`grdseamount `, +:doc:`talwani2d `, +:doc:`talwani3d ` diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 508f379134f..203574c1a7f 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -34,7 +34,7 @@ #define THIS_MODULE_MODERN_NAME "gravprisms" #define THIS_MODULE_LIB "potential" #define THIS_MODULE_PURPOSE "Compute gravity anomalies over 3-D vertical prisms" -#define THIS_MODULE_KEYS "//] creates prisms between surfaces */ + struct GRAVPRISMS_C { /* -C[//]] creates prisms between surfaces */ bool active; char *base, *top, *height; double dz; @@ -154,12 +154,26 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'C': /* Create prisms from layer between two surfaces */ n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active); Ctrl->C.active = true; - nc = sscanf (opt->arg, "%[^,],%[^,],%[^,],%lg", A, B, C, &Ctrl->C.dz); - if (A[0] != '-') Ctrl->C.base = strdup (A); - if (B[0] != '-') Ctrl->C.top = strdup (B); - if (C[0] != '-') Ctrl->C.height = strdup (C); - if (nc < 3) { - GMT_Report (API, GMT_MSG_WARNING, "Option -C: Not enough arguments supplied (only %d detected but 3 or 4 are needed\n", nc); + if ((c = strstr (opt->arg, "+z"))) { + Ctrl->C.dz = atof (&c[2]); + c[0] = '\0'; /* Hide odifier */ + } + nc = sscanf (opt->arg, "%[^,],%[^,],%[^,]", A, B, C); + if (nc == 1) { /* -C[+z] */ + if (A[0] != '-') + Ctrl->C.height = strdup (A); + else { + GMT_Report (API, GMT_MSG_ERROR, "Option -C: Must give full seamount height\n"); + n_errors++; + } + } + else if (nc == 3) { /* -C,,[+z] */ + if (A[0] != '-') Ctrl->C.base = strdup (A); + if (B[0] != '-') Ctrl->C.top = strdup (B); + if (C[0] != '-') Ctrl->C.height = strdup (C); + } + else { + GMT_Report (API, GMT_MSG_ERROR, "Option -C: Incorrect number of arguments: %d\n", nc); n_errors++; } break; @@ -290,8 +304,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -H: Cannot be used with -D.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->D.active && !Ctrl->H.active, "option -C: Need to select either -D or -H.\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && nc < 4, - "Option -C: Requires when -H is selected\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && Ctrl->C.dz == 0.0, + "Option -C: Requires +z when -H is selected\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -299,7 +313,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-C,,,] [-D] [-Ff|n[]|v] [-G] [-H//[+d][+p]] [%s] " + GMT_Usage (API, 0, "usage: %s [-A] [-C[,,][+z]] [-D] [-Ff|n[]|v] [-G] [-H//[+d][+p]] [%s] " "[-M[hz]] [-N] [%s] [%s] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT); @@ -312,9 +326,10 @@ static int usage (struct GMTAPI_CTRL *API, int level) { "input is read. Contains (x,y,z) center coordinates of prisms with optional [ dx dy dz] [rho] if not set via -D and ."); GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); - GMT_Usage (API, 1, "\n-C,,,"); + GMT_Usage (API, 1, "\n-C[,,][+z]"); GMT_Usage (API, -2, "Create prisms from the level to level for full seamount , where arguments are grids or constants. " - "If -H is used then is required for the discretization of rho(r,z). Give = - if is the full height. No will be read."); + "If only grid is given then we assume = 0. Other wise, give all three but set = - if is the full height. " + "If -H is used then is required for the discretization of rho(r,z). No will be read."); GMT_Usage (API, 1, "\n-D"); GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column or computes it via -H]."); GMT_Usage (API, 1, "\n-Ff|n[]|v]"); From 3e486d941093bf37d1b5a21b658bc5b80f5b3186 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 11:41:39 +0000 Subject: [PATCH 05/27] Simplify options for externals --- .../module_supplements_purpose.rst_ | 2 +- .../supplements/potential/gravprisms.rst | 59 ++++-- src/potential/gravprisms.c | 192 ++++++++++++------ 3 files changed, 179 insertions(+), 74 deletions(-) diff --git a/doc/rst/source/supplements/module_supplements_purpose.rst_ b/doc/rst/source/supplements/module_supplements_purpose.rst_ index 8798f91aacd..f31d6dfdfd1 100644 --- a/doc/rst/source/supplements/module_supplements_purpose.rst_ +++ b/doc/rst/source/supplements/module_supplements_purpose.rst_ @@ -36,7 +36,7 @@ .. |gravfft_purpose| replace:: Spectral calculations of gravity, isostasy, admittance, and coherence for grids -.. |gravprisms_purpose| replace:: Compute gravity anomalies over 3-D vertical prisms +.. |gravprisms_purpose| replace:: Compute geopotential anomalies over 3-D vertical prisms .. |grdflexure_purpose| replace:: Compute flexural deformation of 3-D surfaces for various rheologies diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 8472ca7c60e..3e3a3237573 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -14,18 +14,22 @@ Synopsis **gmt gravprisms** [ *table* ] [ |-A| ] -[ |-C|\ [*base*\ /*top*]\ *height*\ [**+z**\ *dz*] ] +[ |-C|\ [**+w**\ *file*][**+z**\ *dz*] ] [ |-D|\ *density* ] ] [ |-E|\ *dx*\ /*dy*\ /*dz* ] [ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] [ |-G|\ *outfile* ] [ |-H|\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] ] [ |SYN_OPT-I| ] +[ |-L|\ *base* ] [ |-M|\ [**h**]\ [**v**] ] [ |-N|\ *trackfile* ] [ |SYN_OPT-R| ] -[ |-Z|\ *level*\|\ *obsgrid* ] +[ |-S|\ *shapegrid* ] +[ |-T|\ *top* ] [ |SYN_OPT-V| ] +[ |-W|\ *avedens* ] +[ |-Z|\ *level*\|\ *obsgrid* ] [ |SYN_OPT-bo| ] [ |SYN_OPT-d| ] [ |SYN_OPT-e| ] @@ -46,10 +50,12 @@ We either read the multi-segment *table* from file (or standard input), which ma 7 colums: The first three are the center *x, y, z* coordinates of the prism, while the next three are the dimensions *dx dy dz* of each prism (see **-E** if all prisms have the same dimesions. Last column may contain individual prism densities (may be overridden by a fixed density contrast -given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate the space between -two surfaces (one of which may be a constant). If a variable density model (**-H**) is selected +given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate +the entire feature (**-S**) or just the volume between two surfaces (one of which may be a constant) +that define a layer (set via **-L** and **-T**). If a variable density model (**-H**) is selected then each vertical prism will be broken into constant-density, stacked sub-prisms using a prescribed -vertical increment *dz*. +vertical increment *dz*, otherwise single tall prisms are created with constant (**-D**) or spatially +varying (**-W**) densities. We can compute anomalies on an equidistant grid (by specifying a new grid with **-R** and **-I** or provide an observation grid with desired elevations) or at arbitrary output points specified via **-N**. Choose between free-air anomalies, vertical @@ -84,14 +90,15 @@ Optional Arguments .. _-C: -**-C**\ [*base*\ /*top*]\ *height*\ [**+z**\ *dz*] ] - Create prisms for the layer between the two surfaces *base* and *top*, with *height* describing the full height of the - feature. Either *base* or *top* may be a constant rather than grid. If *top* equals *height* then - give *height* as -. If only *height* is given then we assume we will approximate the entire feature - from *base* = 0 to *height*. If **-H** is used to compute variable density contrasts then we must - split each prism into a stack of sub-prisms with individual densities; thus *dz* needs to be set for the heights - of these sub-prisms (the first and last sub-prisms in the stack may have their heights adjusted to match the - limits of the surfaces). Without **-H** we only create a single uniform-density prism. +**-C**\ [**+w**\ *file*][**+z**\ *dz*] ] + Create prisms for the entire feature given by **-S**\ *height*, or just for the layer between the two surfaces + specified via **-L**\ *base* and **-T**\ *top*. For layers, either *base* or *top* may be a constant rather + than a grid. If only *height* is given then we assume we will approximate the entire feature from *base* = 0 + to *height*. If **-H** is used to compute variable density contrasts then we must split each prism into a stack + of sub-prisms with individual densities; thus *dz* needs to be set for the heights of these sub-prisms (the first + and last sub-prisms in the stack may have their heights adjusted to match the limits of the surfaces). Without + **-H** we only create a single uniform-density prism, but those prisms may have spatially varying densities via + **-W**. Append **+w**\ *file* to save the prisms to a table. .. _-D: @@ -123,10 +130,16 @@ Optional Arguments **-H**\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] Set reference seamount parameters for an *ad-hoc* variable radial density function with depth. Give the low and high seamount densities in kg/m^3 or g/cm^3 and the fixed reference height *H* in meters. - Use modifers **+d** and **+p** to change the water-pressure-driven flank density increase - over the full reference height [0] and the variable density profile exponent *power* [1, i.e., a linear change]. + Use modifers **+d** and **+p** to change the water-pressure-driven flank density increase over the + full reference height [0] and the variable density profile exponent *power* [1, i.e., a linear change]. + Requires **-S** to know the full height of the seamount. See :doc:`grdseamount ` for more details. +.. _-L: + +**-L**\ *base* + Set the base grid surface for a layer we wish to make prisms for, or give a constant level [0]. + .. _-M: **-M**\ [**h**]\ [**v**] @@ -140,11 +153,27 @@ Optional Arguments is used there are no grids and the output data records are written to standard output (see **-bo** for binary output). If *trackfile* has 3 columns we take the *z* value as our observation level; this level may be overridden via **-Z**. +.. _-S: + +**-S**\ *height* + Set the grid with the full seamount height, either to make prisms from or as required by **-H**. + +.. _-T: + +**-T**\ *top* + Set the top grid surface for a layer we wish to make prisms for, or give a constant level. + .. |Add_-V| replace:: |Add_-V_links| .. include:: /explain_-V.rst_ :start-after: **Syntax** :end-before: **Description** +.. _-W: + +**-W**\ *avedens* + Give the name of a grid with spatially varying, vertically-averaged prism densities (see :ref:`Grid File Formats + `). Requires **-C** and must be co-registered with such grids. + .. _-Z: **-Z**\ *level*\|\ *obsgrid* diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 203574c1a7f..6872cf1d9a4 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -15,14 +15,15 @@ * Contact info: www.generic-mapping-tools.org *--------------------------------------------------------------------*/ /* - * Authord: Paul Wessel and Seung-Sep Kim + * Authord: Paul Wessel * Date: 12-MAR-2022 * * - * Calculates gravity due to any number of individual vertical prisms - * that may have their own unique dimensions (or all are constant) and - * individual density contrasts (or all the same). - * Based on methods by + * Calculates geopotential anomalies due to any number of individual + * vertical prisms that may have their own unique dimensions (or all + * are constant) and individual density contrasts (or all the same). + * For vertically varying density we break prisms into a stack of + * short sub-prisms with constant density. * * Accelerated with OpenMP; see -x. */ @@ -33,11 +34,13 @@ #define THIS_MODULE_CLASSIC_NAME "gravprisms" #define THIS_MODULE_MODERN_NAME "gravprisms" #define THIS_MODULE_LIB "potential" -#define THIS_MODULE_PURPOSE "Compute gravity anomalies over 3-D vertical prisms" -#define THIS_MODULE_KEYS "//]] creates prisms between surfaces */ + struct GRAVPRISMS_C { /* -C[+w][+z] creates prisms between surfaces set in -L -S -T */ bool active; - char *base, *top, *height; + char *file; double dz; } C; struct GRAVPRISMS_E { /* -E/ fixed prism dimensions [read from file] */ @@ -91,10 +94,14 @@ struct GRAVPRISMS_CTRL { double del_rho; /* Computed as rho_h - rho_l */ double p1; /* Will be p + 1 to avoid addition in loops */ } H; -struct GRAVPRISMS_I { /* -Idx[/dy] */ + struct GRAVPRISMS_I { /* -Idx[/dy] */ bool active; double inc[2]; } I; + struct GRAVPRISMS_L { /* Low (base) surface(x,y) file */ + bool active; + char *file; + } L; struct GRAVPRISMS_M { /* -Mh|z */ bool active[2]; /* True if km, else m */ } M; @@ -102,6 +109,18 @@ struct GRAVPRISMS_I { /* -Idx[/dy] */ bool active; char *file; } N; + struct GRAVPRISMS_S { /* Full surface(x,y) file */ + bool active; + char *file; + } S; + struct GRAVPRISMS_T { /* Top surface(x,y) file */ + bool active; + char *file; + } T; + struct GRAVPRISMS_W { /* Variable rho(x,y) file */ + bool active; + char *file; + } W; struct GRAVPRISMS_Z { /* Observation level file or constant */ bool active; double level; @@ -125,10 +144,11 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *C) { /* Deallocate control structure */ if (!C) return; gmt_M_str_free (C->N.file); - gmt_M_str_free (C->C.base); - gmt_M_str_free (C->C.top); - gmt_M_str_free (C->C.height); + gmt_M_str_free (C->C.file); gmt_M_str_free (C->G.file); + gmt_M_str_free (C->S.file); + gmt_M_str_free (C->T.file); + gmt_M_str_free (C->W.file); gmt_M_str_free (C->Z.file); gmt_M_free (GMT, C); } @@ -154,28 +174,23 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'C': /* Create prisms from layer between two surfaces */ n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active); Ctrl->C.active = true; - if ((c = strstr (opt->arg, "+z"))) { - Ctrl->C.dz = atof (&c[2]); - c[0] = '\0'; /* Hide odifier */ - } - nc = sscanf (opt->arg, "%[^,],%[^,],%[^,]", A, B, C); - if (nc == 1) { /* -C[+z] */ - if (A[0] != '-') - Ctrl->C.height = strdup (A); - else { - GMT_Report (API, GMT_MSG_ERROR, "Option -C: Must give full seamount height\n"); - n_errors++; + if ((c = gmt_first_modifier (GMT, opt->arg, "wz"))) { + unsigned int pos = 0; + char txt[GMT_LEN256] = {""}; + while (gmt_getmodopt (GMT, 'C', c, "wz", &pos, txt, &n_errors) && n_errors == 0) { + switch (txt[0]) { + case 'w': /* Write created prisms to given file */ + Ctrl->C.file = strdup (&txt[1]); + break; + case 'z': /* Set vertical increment*/ + Ctrl->C.dz = atof (&txt[1]); + break; + default: + n_errors++; + break; + } } } - else if (nc == 3) { /* -C,,[+z] */ - if (A[0] != '-') Ctrl->C.base = strdup (A); - if (B[0] != '-') Ctrl->C.top = strdup (B); - if (C[0] != '-') Ctrl->C.height = strdup (C); - } - else { - GMT_Report (API, GMT_MSG_ERROR, "Option -C: Incorrect number of arguments: %d\n", nc); - n_errors++; - } break; case 'D': n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active); @@ -237,12 +252,18 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT if (Ctrl->H.rho_h < 10.0) Ctrl->H.rho_h *= 1000; /* Gave units of g/cm^3 */ Ctrl->H.del_rho = Ctrl->H.rho_h - Ctrl->H.rho_l; Ctrl->H.p1 = Ctrl->H.p + 1; + if (c) c[0] = '+'; /* Restore modifiers */ break; case 'I': /* Grid iincrements */ n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active); Ctrl->I.active = true; n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg); break; + case 'L': /* Low (base) file (or constant) given */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active); + Ctrl->L.active = true; + Ctrl->L.file = strdup (opt->arg); + break; case 'M': /* Horizontal and vertical length units */ k = 0; while (opt->arg[k]) { @@ -268,6 +289,21 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->N.active = true; Ctrl->N.file = strdup (opt->arg); break; + case 'S': /* Full surface grid file */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active); + Ctrl->S.active = true; + Ctrl->S.file = strdup (opt->arg); + break; + case 'T': /* top of layer file (or constant) given */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active); + Ctrl->T.active = true; + Ctrl->T.file = strdup (opt->arg); + break; + case 'W': /* Grid with vertically-averaged density contrasts */ + n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active); + Ctrl->W.active = true; + Ctrl->W.file = strdup (opt->arg); + break; case 'Z': /* Observation level(s) */ n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active); Ctrl->Z.active = true; @@ -302,10 +338,16 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -N: Must specify output gridfile name.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->H.active, "Option -H: Cannot be used with -D.\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->D.active && !Ctrl->H.active, - "option -C: Need to select either -D or -H.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->W.active && (Ctrl->D.active || Ctrl->H.active), + "Option -W: Cannot be used with -D or -H.\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->D.active && !Ctrl->H.active && !Ctrl->W.active, + "option -C: Need to select either -D, -H or -W.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && Ctrl->C.dz == 0.0, "Option -C: Requires +z when -H is selected\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && !Ctrl->S.active, + "Option -C: Requires -S when -H is set\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->L.active && !Ctrl->T.active, + "Option -L: Requires -T\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -313,8 +355,9 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-C[,,][+z]] [-D] [-Ff|n[]|v] [-G] [-H//[+d][+p]] [%s] " - "[-M[hz]] [-N] [%s] [%s] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", + GMT_Usage (API, 0, "usage: %s [-A] [-C[+w][+z]] [-D] [-Ff|n[]|v] [-G] " + "[-H//[+d][+p]] [-L] [%s] [-M[hz]] [-N] [%s] [-S] " + "[-T] [%s] [-W] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT); @@ -346,6 +389,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "+d Density increase (kg/m^3 or g/cm^3) due to water pressure over full depth implied by [0]."); GMT_Usage (API, 3, "+p Exponential coefficient (> 0) for density change with burial depth [1 (linear)]."); GMT_Option (API, "I"); + GMT_Usage (API, 1, "\n-L"); + GMT_Usage (API, -2, "Set the lower (base) surface grid or constant of a layer to create prisms for; requires -C, -S, and -T [0]"); GMT_Usage (API, 1, "\n-M[hz]"); GMT_Usage (API, -2, "Change distance units used, via one or two directives:"); GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); @@ -353,7 +398,14 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 1, "\n-N"); GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally z may be read as well. No grid " "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); - GMT_Option (API, "R,V"); + GMT_Option (API, "R"); + GMT_Usage (API, 1, "\n-S"); + GMT_Usage (API, -2, "Set the full surface grid height for making prisms of the entire feature (or if -H is set) (or use -L and -T to select a layer); requires -C"); + GMT_Option (API, "V"); + GMT_Usage (API, 1, "\n-T"); + GMT_Usage (API, -2, "Set the top surface grid or constant of a layer to create prisms for; requires -C, -L and -T"); + GMT_Usage (API, 1, "\n-W"); + GMT_Usage (API, -2, "Give grid with vertically-averaged spatially varying densities; requires co-registered grids with -C."); GMT_Usage (API, 1, "\n-Z"); GMT_Usage (API, -2, "Set observation level for output locations [0]. " "Append either a constant or the name of a grid file with variable levels. " @@ -551,29 +603,39 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_enable_threads (GMT); /* Set number of active threads, if supported */ if (Ctrl->C.active) { /* Need to create prisms from two surfaces first */ - struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *G = NULL; + struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *Rho = NULL; double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid; size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; - if (access (Ctrl->C.base, F_OK)) /* No file, just a constant */ - base = atof (Ctrl->C.base); - else if ((B = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.base, NULL)) == NULL) - Return (API->error); - if (access (Ctrl->C.top, F_OK)) /* No file, just a constant */ - top = atof (Ctrl->C.base); - else if ((T = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.top, NULL)) == NULL) - Return (API->error); - if (Ctrl->C.height == NULL) /* No file, so height given by top */ - H = T; - else if ((H = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->C.height, NULL)) == NULL) - Return (API->error); - G = (B) ? B : T; - dx = 0.5 * G->header->inc[GMT_X]; dy = 0.5 * G->header->inc[GMT_Y]; - rho = Ctrl->D.rho; /* Only helpful if -D is set */ + if (Ctrl->L.active) { /* Specifed layer base */ + if (access (Ctrl->L.file, F_OK)) /* No file, just a constant */ + base = atof (Ctrl->L.file); + else if ((B = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->L.file, NULL)) == NULL) + Return (API->error); + } + if (Ctrl->T.active) { /* Specified layer top */ + if (access (Ctrl->T.file, F_OK)) /* No file, just a constant */ + top = atof (Ctrl->T.file); + else if ((T = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->T.file, NULL)) == NULL) + Return (API->error); + } + if (Ctrl->S.active) { /* Full shape needed (via -H or as the entire grid given) */ + if ((H = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->S.file, NULL)) == NULL) + Return (API->error); + if (!Ctrl->T.active) T = H; /* Top and height are the same if just -S is set */ + } + else /* Use H for info for any of the grids */ + H = (B) ? B : T; + if (Ctrl->W.active) { /* Got spatially varying average density contrasts */ + if ((Rho = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->W.file, NULL)) == NULL) + Return (API->error); + } + dx = 0.5 * H->header->inc[GMT_X]; dy = 0.5 * H->header->inc[GMT_Y]; + rho = Ctrl->D.rho; /* Only initialized if -D is set */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, NULL, n_alloc, double); - gmt_M_grd_loop (GMT, G, row, col, node) { + gmt_M_grd_loop (GMT, H, row, col, node) { z2 = (T) ? T->data[node] : top; z1 = (B) ? B->data[node] : base; if (z2 <= z1) continue; /* No layer thickness detected */ @@ -587,15 +649,17 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { z_mid = 0.5 * (z_prev + z_next); /* Middle of prism - used to look up density */ rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); } - else /* Constant density rho (set above), just need a single prism per location */ + else { /* Constant density rho (set above via -D) or by Rho (via -W), just need a single prism per location */ z_next = z2; + if (Ctrl->W.active) rho = Rho->data[node]; /* Update constant density for this prism */ + } if (n_prisms == n_alloc) { /* Need to allocate more memory for the prisms */ n_alloc <<= 1; /* Double it */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_alloc, double); } - prism[0][n_prisms] = G->x[col] - dx; prism[1][n_prisms] = G->x[col] + dx; - prism[2][n_prisms] = G->y[row] - dy; prism[3][n_prisms] = G->y[row] + dy; + prism[0][n_prisms] = H->x[col] - dx; prism[1][n_prisms] = H->x[col] + dx; + prism[2][n_prisms] = H->y[row] - dy; prism[3][n_prisms] = H->y[row] + dy; prism[4][n_prisms] = z_prev; prism[5][n_prisms] = z_next; prism[6][n_prisms] = rho; n_prisms++; @@ -615,6 +679,18 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (H && GMT_Destroy_Data (API, &H) != GMT_NOERROR) { Return (GMT_MEMORY_ERROR); } + if (Ctrl->C.file) { /* Wish to write prisms to a data table */ + uint64_t dim[GMT_DIM_SIZE] = {1, 1, n_prisms, 7}; + if ((P = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) + Return (API->error); + S = P->table[0]->segment[0]; /* The only segment in the table */ + for (k = 0; k < 7; k++) /* Copy over each column */ + gmt_M_memcpy (S->data[k], prism[k], n_prisms, double); + if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->C.file, P) != GMT_NOERROR) { + GMT_Report (API, GMT_MSG_ERROR, "Unable to write prisms to file %s\n", Ctrl->C.file); + Return (GMT_RUNTIME_ERROR); + } + } } else { /* Read prisms from stdin or input file(s) */ unsigned int n_expected = 7; /* Max number of columns */ From 461945cc12eedbd7b0f89413ae5039866c0574a7 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 11:54:16 +0000 Subject: [PATCH 06/27] Polish --- .../supplements/potential/gravprisms.rst | 7 ++-- src/potential/gravprisms.c | 34 +++++++++++++------ 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 3e3a3237573..57f879e9246 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -14,7 +14,7 @@ Synopsis **gmt gravprisms** [ *table* ] [ |-A| ] -[ |-C|\ [**+w**\ *file*][**+z**\ *dz*] ] +[ |-C|\ [**+q**][**+w**\ *file*][**+z**\ *dz*] ] [ |-D|\ *density* ] ] [ |-E|\ *dx*\ /*dy*\ /*dz* ] [ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] @@ -90,7 +90,7 @@ Optional Arguments .. _-C: -**-C**\ [**+w**\ *file*][**+z**\ *dz*] ] +**-C**\ [**+q**][**+w**\ *file*][**+z**\ *dz*] ] Create prisms for the entire feature given by **-S**\ *height*, or just for the layer between the two surfaces specified via **-L**\ *base* and **-T**\ *top*. For layers, either *base* or *top* may be a constant rather than a grid. If only *height* is given then we assume we will approximate the entire feature from *base* = 0 @@ -98,7 +98,8 @@ Optional Arguments of sub-prisms with individual densities; thus *dz* needs to be set for the heights of these sub-prisms (the first and last sub-prisms in the stack may have their heights adjusted to match the limits of the surfaces). Without **-H** we only create a single uniform-density prism, but those prisms may have spatially varying densities via - **-W**. Append **+w**\ *file* to save the prisms to a table. + **-W**. Append **+w**\ *file* to save the prisms to a table, and optionally use **+q** to quit execution once + the file has been saved, i.e., no geopotential calculations will take place. .. _-D: diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 6872cf1d9a4..8812077930a 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -69,8 +69,10 @@ struct GRAVPRISMS_CTRL { bool active; double rho; } D; - struct GRAVPRISMS_C { /* -C[+w][+z] creates prisms between surfaces set in -L -S -T */ + struct GRAVPRISMS_C { /* -C[+q][+w][+z] creates prisms between surfaces set in -L -S -T */ bool active; + bool quit; + bool dump; char *file; double dz; } C; @@ -154,9 +156,8 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *C) { /* Dea } static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT_OPTION *options) { - int nc = 0; unsigned int k, n_errors = 0; - char *c = NULL, A[GMT_LEN256] = {""}, B[GMT_LEN256] = {""}, C[GMT_LEN256] = {""}; + char *c = NULL; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; @@ -174,13 +175,17 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'C': /* Create prisms from layer between two surfaces */ n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active); Ctrl->C.active = true; - if ((c = gmt_first_modifier (GMT, opt->arg, "wz"))) { + if ((c = gmt_first_modifier (GMT, opt->arg, "qwz"))) { unsigned int pos = 0; char txt[GMT_LEN256] = {""}; - while (gmt_getmodopt (GMT, 'C', c, "wz", &pos, txt, &n_errors) && n_errors == 0) { + while (gmt_getmodopt (GMT, 'C', c, "qwz", &pos, txt, &n_errors) && n_errors == 0) { switch (txt[0]) { + case 'q': /* Quit once prism file has been written */ + Ctrl->C.quit = true; + break; case 'w': /* Write created prisms to given file */ - Ctrl->C.file = strdup (&txt[1]); + Ctrl->C.dump = true; + if (txt[1]) Ctrl->C.file = strdup (&txt[1]); break; case 'z': /* Set vertical increment*/ Ctrl->C.dz = atof (&txt[1]); @@ -348,6 +353,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -C: Requires -S when -H is set\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->L.active && !Ctrl->T.active, "Option -L: Requires -T\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.quit && !Ctrl->C.dump, + "Option -C: Modifer +q requires +w\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -355,7 +362,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-C[+w][+z]] [-D] [-Ff|n[]|v] [-G] " + GMT_Usage (API, 0, "usage: %s [-A] [-C[+q][+w][+z]] [-D] [-Ff|n[]|v] [-G] " "[-H//[+d][+p]] [-L] [%s] [-M[hz]] [-N] [%s] [-S] " "[-T] [%s] [-W] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, @@ -369,10 +376,11 @@ static int usage (struct GMTAPI_CTRL *API, int level) { "input is read. Contains (x,y,z) center coordinates of prisms with optional [ dx dy dz] [rho] if not set via -D and ."); GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); - GMT_Usage (API, 1, "\n-C[,,][+z]"); - GMT_Usage (API, -2, "Create prisms from the level to level for full seamount , where arguments are grids or constants. " - "If only grid is given then we assume = 0. Other wise, give all three but set = - if is the full height. " - "If -H is used then is required for the discretization of rho(r,z). No will be read."); + GMT_Usage (API, 1, "\n-C[+q][+w][+z]"); + GMT_Usage (API, -2, "Create prisms from the (-L) level to (-T) level, or for the full seamount (-S). No will be read. Modifiers are:"); + GMT_Usage (API, 3, "+q Quit execution once prism file has been written (see +w). No geopotential calculations are preformed."); + GMT_Usage (API, 3, "+w Write the created prisms to "); + GMT_Usage (API, 3, "+z Set increment for discretization of rho(r,z) when -H is used."); GMT_Usage (API, 1, "\n-D"); GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column or computes it via -H]."); GMT_Usage (API, 1, "\n-Ff|n[]|v]"); @@ -690,6 +698,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_ERROR, "Unable to write prisms to file %s\n", Ctrl->C.file); Return (GMT_RUNTIME_ERROR); } + if (Ctrl->C.quit) goto end_it_all; } } else { /* Read prisms from stdin or input file(s) */ @@ -901,6 +910,9 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (API->error); } } + +end_it_all: + /* Free the prism information */ for (k = 0; k < 7; k++) gmt_M_free (GMT, prism[k]); From 80b34024221d13691c750b4a37328685fe58d790 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 12:02:07 +0000 Subject: [PATCH 07/27] Update gravprisms.c --- src/potential/gravprisms.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 8812077930a..27dcb08d7c4 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -335,7 +335,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -Z: Cannot also specify -R -I if a level grid has been supplied\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->Z.mode == 1, "Option -Z: Cannot also specify -N if a level grid has been supplied\n"); - n_errors += gmt_M_check_condition (GMT, !Ctrl->N.active && !Ctrl->G.active, + n_errors += gmt_M_check_condition (GMT, !Ctrl->N.active && !Ctrl->G.active && !Ctrl->C.quit, "Option -G: Must specify output gridfile name.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && !Ctrl->G.file, "Option -G: Must specify output gridfile name.\n"); @@ -678,6 +678,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_prisms, double); GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms constructed: %" PRIu64 "\n", n_prisms); + if (!Ctrl->T.active && T) T = NULL; /* Undo what we set earlier */ if (B && GMT_Destroy_Data (API, &B) != GMT_NOERROR) { Return (GMT_MEMORY_ERROR); } From d80a8c3bb8815ab885b229c2db11b93ffab95faf Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 12:32:13 +0000 Subject: [PATCH 08/27] Update gravprisms.c --- src/potential/gravprisms.c | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 27dcb08d7c4..b818e7c3845 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -693,8 +693,15 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if ((P = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error); S = P->table[0]->segment[0]; /* The only segment in the table */ - for (k = 0; k < 7; k++) /* Copy over each column */ - gmt_M_memcpy (S->data[k], prism[k], n_prisms, double); + for (k = 0; k < n_prisms; k++) { /* Unscramble edges to coordinates and dimensions */ + S->data[0][k] = 0.5 * (prism[1][k] + prism[0][k]); /* Get x */ + S->data[1][k] = 0.5 * (prism[3][k] + prism[2][k]); /* Get y */ + S->data[2][k] = 0.5 * (prism[5][k] + prism[4][k]); /* Get z */ + S->data[3][k] = prism[1][k] - prism[0][k]; /* Get dx */ + S->data[4][k] = prism[3][k] - prism[2][k]; /* Get dy */ + S->data[5][k] = prism[5][k] - prism[4][k]; /* Get dz */ + } + gmt_M_memcpy (S->data[6], prism[6], n_prisms, double); /* Copy over densities */ if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->C.file, P) != GMT_NOERROR) { GMT_Report (API, GMT_MSG_ERROR, "Unable to write prisms to file %s\n", Ctrl->C.file); Return (GMT_RUNTIME_ERROR); @@ -718,8 +725,8 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if ((P = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, 0, GMT_READ_NORMAL, NULL, NULL, NULL)) == NULL) { Return (API->error); } - /* To avoid the need to loop over tables and segments, we extrac the data into a separate array and comptute prism edges instead */ - for (k = 0; k < 7; k++) + /* To avoid the need to loop over tables and segments, we extrac the data into a separate array and compute prism edges instead */ + for (k = 0; k < n_prisms; k++) prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free the dataset */ From 18440aa38188ac9ba3ee9a4f17654267e8a9940a Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 18:50:43 +0000 Subject: [PATCH 09/27] Clarify -M and -A --- .../supplements/potential/gravprisms.rst | 17 ++- src/potential/gravprisms.c | 138 +++++++++++------- 2 files changed, 96 insertions(+), 59 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 57f879e9246..b51834a364b 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -86,7 +86,8 @@ Optional Arguments .. _-A: **-A** - The *z*-axis should be positive upwards [Default is down]. + The *z*-axis should be positive upwards [Default is down]. All internal *z*-values will have their sign + changed. **Note**: Output *z*-values are not affected. .. _-C: @@ -139,13 +140,15 @@ Optional Arguments .. _-L: **-L**\ *base* - Set the base grid surface for a layer we wish to make prisms for, or give a constant level [0]. + Give name of the base grid surface for a layer we wish to fill with prisms, or give a constant *z*-level [0]. .. _-M: **-M**\ [**h**]\ [**v**] Sets distance units used. Append **h** to indicate that both horizontal distances are in km [m], - and append **z** to indicate vertical distances are in km [m]. + and append **z** to indicate vertical distances are in km [m]. If selected, we will internally + convert any affected distance provided by data input or command line options to meters. **Note**: + Any output will retain the original units. .. _-N: @@ -157,12 +160,12 @@ Optional Arguments .. _-S: **-S**\ *height* - Set the grid with the full seamount height, either to make prisms from or as required by **-H**. + Give name of grid with the full seamount heights, either for making prisms or as required by **-H**. .. _-T: **-T**\ *top* - Set the top grid surface for a layer we wish to make prisms for, or give a constant level. + Give name of the top surface grid for a layer we wish to fill with prisms, or give a constant *z*-level. .. |Add_-V| replace:: |Add_-V_links| .. include:: /explain_-V.rst_ @@ -172,8 +175,8 @@ Optional Arguments .. _-W: **-W**\ *avedens* - Give the name of a grid with spatially varying, vertically-averaged prism densities (see :ref:`Grid File Formats - `). Requires **-C** and must be co-registered with such grids. + Give name of a grid with spatially varying, vertically-averaged prism densities. Requires **-C** and + the grid must be co-registered with the grid provided by **-S**. .. _-Z: diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index b818e7c3845..8826fe78558 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -403,6 +403,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, -2, "Change distance units used, via one or two directives:"); GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); GMT_Usage (API, 3, "z: All z-distances are given in km [meters]."); + GMT_Usage (API, -2, "Note: All horixontal and/or vertical quantities will be affected by a factor of 1000"); GMT_Usage (API, 1, "\n-N"); GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally z may be read as well. No grid " "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); @@ -493,7 +494,7 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub g221 = -(dz2 * atan ((dx1 * dy2) / (dz2 * R221)) - dx1 * log (R221 + dy2) - dy2 * log (R221 + dx1)); g222 = +(dz2 * atan ((dx2 * dy2) / (dz2 * R222)) - dx2 * log (R222 + dy2) - dy2 * log (R222 + dx2)); - g = rho * GRAVITATIONAL_CONST_MGAL * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); + g = -rho * GRAVITATIONAL_CONST_MGAL * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); return (g); } @@ -579,7 +580,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { bool flat_earth = false; char *uname[2] = {"meter", "km"}, *kind[3] = {"FAA", "VGG", "GEOID"}, remark[GMT_LEN64] = {""}; - double z_level, rho, dx, dy, dz, lat, G0, *prism[7]; + double z_level, rho, dx, dy, dz, lat, G0, scl_xy, scl_z, i_scl_xy, i_scl_z, *prism[7]; double (*eval) (double, double, double, uint64_t, double **, double); struct GRAVPRISMS_CTRL *Ctrl = NULL; @@ -610,6 +611,13 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_enable_threads (GMT); /* Set number of active threads, if supported */ + /* Handle any -M[hv] settings so that internally we have meters */ + scl_xy = (Ctrl->M.active[GRAVPRISMS_HOR]) ? METERS_IN_A_KM : 1.0; + scl_z = (Ctrl->M.active[GRAVPRISMS_VER]) ? METERS_IN_A_KM : 1.0; + if (Ctrl->A.active) scl_z = -scl_z; + i_scl_xy = 1.0 / scl_xy; /* Scale use for output horizontal distances */ + i_scl_z = 1.0 / scl_z; /* Scale use for output vertical distances */ + if (Ctrl->C.active) { /* Need to create prisms from two surfaces first */ struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *Rho = NULL; double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid; @@ -666,9 +674,10 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_alloc, double); } - prism[0][n_prisms] = H->x[col] - dx; prism[1][n_prisms] = H->x[col] + dx; - prism[2][n_prisms] = H->y[row] - dy; prism[3][n_prisms] = H->y[row] + dy; - prism[4][n_prisms] = z_prev; prism[5][n_prisms] = z_next; + /* Here we ensure prism is in meters */ + prism[0][n_prisms] = scl_xy * (H->x[col] - dx); prism[1][n_prisms] = scl_xy * (H->x[col] + dx); + prism[2][n_prisms] = scl_xy * (H->y[row] - dy); prism[3][n_prisms] = scl_xy * (H->y[row] + dy); + prism[4][n_prisms] = scl_z * z_prev; prism[5][n_prisms] = scl_z * z_next; prism[6][n_prisms] = rho; n_prisms++; z_prev = z_next; /* The the top of this prism be the bottom of the next */ @@ -680,31 +689,35 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms constructed: %" PRIu64 "\n", n_prisms); if (!Ctrl->T.active && T) T = NULL; /* Undo what we set earlier */ if (B && GMT_Destroy_Data (API, &B) != GMT_NOERROR) { - Return (GMT_MEMORY_ERROR); + error = GMT_MEMORY_ERROR; + goto end_it_all; } if (T && GMT_Destroy_Data (API, &T) != GMT_NOERROR) { - Return (GMT_MEMORY_ERROR); + error = GMT_MEMORY_ERROR; + goto end_it_all; } if (H && GMT_Destroy_Data (API, &H) != GMT_NOERROR) { - Return (GMT_MEMORY_ERROR); + error = GMT_MEMORY_ERROR; + goto end_it_all; } if (Ctrl->C.file) { /* Wish to write prisms to a data table */ uint64_t dim[GMT_DIM_SIZE] = {1, 1, n_prisms, 7}; if ((P = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) Return (API->error); S = P->table[0]->segment[0]; /* The only segment in the table */ - for (k = 0; k < n_prisms; k++) { /* Unscramble edges to coordinates and dimensions */ - S->data[0][k] = 0.5 * (prism[1][k] + prism[0][k]); /* Get x */ - S->data[1][k] = 0.5 * (prism[3][k] + prism[2][k]); /* Get y */ - S->data[2][k] = 0.5 * (prism[5][k] + prism[4][k]); /* Get z */ - S->data[3][k] = prism[1][k] - prism[0][k]; /* Get dx */ - S->data[4][k] = prism[3][k] - prism[2][k]; /* Get dy */ - S->data[5][k] = prism[5][k] - prism[4][k]; /* Get dz */ + for (k = 0; k < n_prisms; k++) { /* Unscramble edges to coordinates and dimensions and possibly convert back to km */ + S->data[0][k] = 0.5 * (prism[1][k] + prism[0][k]) * i_scl_xy; /* Get x */ + S->data[1][k] = 0.5 * (prism[3][k] + prism[2][k]) * i_scl_xy; /* Get y */ + S->data[2][k] = 0.5 * (prism[5][k] + prism[4][k]) * i_scl_xy; /* Get z */ + S->data[3][k] = (prism[1][k] - prism[0][k]) * i_scl_xy; /* Get dx */ + S->data[4][k] = (prism[3][k] - prism[2][k]) * i_scl_xy; /* Get dy */ + S->data[5][k] = (prism[5][k] - prism[4][k]) * i_scl_xy; /* Get dz */ } gmt_M_memcpy (S->data[6], prism[6], n_prisms, double); /* Copy over densities */ if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->C.file, P) != GMT_NOERROR) { GMT_Report (API, GMT_MSG_ERROR, "Unable to write prisms to file %s\n", Ctrl->C.file); - Return (GMT_RUNTIME_ERROR); + error = GMT_RUNTIME_ERROR; + goto end_it_all; } if (Ctrl->C.quit) goto end_it_all; } @@ -726,7 +739,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { Return (API->error); } /* To avoid the need to loop over tables and segments, we extrac the data into a separate array and compute prism edges instead */ - for (k = 0; k < n_prisms; k++) + for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free the dataset */ @@ -743,12 +756,13 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { dz = 0.5 * S->data[5][row]; } if (!Ctrl->D.active) rho = S->data[col][row]; - prism[0][k] = S->data[GMT_X][row] - dx; /* Store the x-coordinates of the x-edges of prism */ - prism[1][k] = S->data[GMT_X][row] + dx; - prism[2][k] = S->data[GMT_Y][row] - dy; /* Store the y-coordinates of the y-edges of prism */ - prism[3][k] = S->data[GMT_Y][row] + dy; - prism[4][k] = S->data[GMT_Z][row] - dz; /* Store the z-coordinates of the z-edges of prism */ - prism[5][k] = S->data[GMT_Z][row] + dz; + /* Here we ensure prism is in meters */ + prism[0][k] = (S->data[GMT_X][row] - dx) * scl_xy; /* Store the x-coordinates of the x-edges of prism */ + prism[1][k] = (S->data[GMT_X][row] + dx) * scl_xy; + prism[2][k] = (S->data[GMT_Y][row] - dy) * scl_xy; /* Store the y-coordinates of the y-edges of prism */ + prism[3][k] = (S->data[GMT_Y][row] + dy) * scl_xy; + prism[4][k] = (S->data[GMT_Z][row] - dz) * scl_z; /* Store the z-coordinates of the z-edges of prism */ + prism[5][k] = (S->data[GMT_Z][row] + dz) * scl_z; prism[6][k] = rho; /* Store the fixed or variable density of prism */ } } @@ -756,34 +770,48 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { n_prisms = k; /* Free the data set */ if (GMT_Destroy_Data (API, &P) != GMT_NOERROR) { - Return (GMT_MEMORY_ERROR); + error = GMT_MEMORY_ERROR; + goto end_it_all; } GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms read: %" PRIu64 "\n", n_prisms); } if (n_prisms == 0) { GMT_Report (API, GMT_MSG_WARNING, "No prisms found - exiting\n"); - Return (GMT_NOERROR); + goto end_it_all; } if (Ctrl->Z.mode == 1) { /* Got grid with observation levels which also sets output locations; it could also set -fg so do this first */ - if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Z.file, NULL)) == NULL) - Return (API->error); + if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->Z.file, NULL)) == NULL) { + error = API->error; + goto end_it_all; + } if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); /* Mid-latitude needed if geoid is selected */ } else if (GMT->common.R.active[RSET]) { /* Gave -R -I [-r] and possibly -fg indirectly via geographic coordinates in -R */ if ((G = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, - GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) - Return (API->error); + GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) { + error = API->error; + goto end_it_all; + } if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); /* Mid-latitude needed if geoid is selected */ } else { /* Got a dataset with output locations via -N */ + unsigned int n_expected = 4; /* Max number of columns */ + if (Ctrl->Z.active) n_expected--; gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from the -N file */ - if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) - Return (API->error); + if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { + error = API->error; + goto end_it_all; + } + if ((D = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_READ_NORMAL, NULL, Ctrl->N.file, NULL)) == NULL) { + error = API->error; + goto end_it_all; + } if (D->n_columns < 2) { GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->N.file, (int)D->n_columns); - Return (GMT_DIM_TOO_SMALL); + error = GMT_DIM_TOO_SMALL; + goto end_it_all; } gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */ if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (D->min[GMT_Y] + D->max[GMT_Y]); /* Mid-latitude needed if geoid is selected */ @@ -793,7 +821,8 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (flat_earth && Ctrl->M.active[GRAVPRISMS_HOR]) { GMT_Report (API, GMT_MSG_ERROR, "Option -M: Cannot specify both geographic coordinates (degrees) AND -Mh\n"); - Return (GMT_RUNTIME_ERROR); + error = GMT_RUNTIME_ERROR; + goto end_it_all; } if (Ctrl->A.active) Ctrl->Z.level = -Ctrl->Z.level; @@ -825,28 +854,31 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (Ctrl->N.active) { /* Single loop over specified output locations */ unsigned int wmode = GMT_ADD_DEFAULT; - double scl = (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) ? (1.0 / METERS_IN_A_KM) : 1.0; /* Perhaps convert to km */ - double out[4]; + double out[4]; /* x, y, z, g|n|v */ struct GMT_RECORD *Rec = gmt_new_record (GMT, out, NULL); /* Must register Ctrl->G.file first since we are going to writing rec-by-rec */ if (Ctrl->G.active) { int out_ID; if ((out_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_OUT, NULL, Ctrl->G.file)) == GMT_NOTSET) { - Return (API->error); + error = API->error; + goto end_it_all; } wmode = GMT_ADD_EXISTING; } if ((error = GMT_Set_Columns (API, GMT_OUT, 4, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { - Return (error); + goto end_it_all; } if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, wmode, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */ - Return (API->error); + error = API->error; + goto end_it_all; } if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */ - Return (API->error); + error = API->error; + goto end_it_all; } if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */ - Return (API->error); + error = API->error; + goto end_it_all; } if (D->n_segments > 1) gmt_set_segmentheader (GMT, GMT_OUT, true); for (tbl = 0; tbl < D->n_tables; tbl++) { @@ -857,13 +889,13 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_prep_tmp_arrays (GMT, GMT_OUT, S->n_rows, 1); /* Init or reallocate tmp vector */ #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl,n_prisms,prism,flat_earth,G0) +#pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl_xy,scl_z,n_prisms,prism,flat_earth,G0) #endif /* Separate the calculation from the output in two separate row-loops since cannot do rec-by-rec output * with OpenMP due to race condiations that would mess up the output order */ for (row = 0; row < (int64_t)S->n_rows; row++) { /* Calculate attraction at all output locations for this segment */ z_level = (S->n_columns == 3 && !Ctrl->Z.active) ? S->data[GMT_Z][row] : Ctrl->Z.level; /* Default observation z level unless provided in input file */ - GMT->hidden.mem_coord[GMT_X][row] = eval (S->data[GMT_X][row] * scl, S->data[GMT_Y][row] * scl, z_level, n_prisms, prism, G0); + GMT->hidden.mem_coord[GMT_X][row] = eval (S->data[GMT_X][row] * scl_xy, S->data[GMT_Y][row] * scl_xy, z_level * scl_z, n_prisms, prism, G0); } /* This loop is not under OpenMP */ out[GMT_Z] = Ctrl->Z.level; /* Default observation z level unless provided in input file */ @@ -878,23 +910,22 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { } gmt_M_free (GMT, Rec); if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */ - Return (API->error); + error = API->error; + goto end_it_all; } } else { /* Dealing with a grid */ openmp_int row, col, n_columns = (openmp_int)G->header->n_columns, n_rows = (openmp_int)G->header->n_rows; /* To shut up compiler warnings */ double y_obs, *x_obs = gmt_M_memory (GMT, NULL, G->header->n_columns, double); for (col = 0; col < n_columns; col++) { - x_obs[col] = gmt_M_grd_col_to_x (GMT, col, G->header); - if (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) x_obs[col] /= METERS_IN_A_KM; /* Convert to km */ + x_obs[col] = scl_xy * gmt_M_grd_col_to_x (GMT, col, G->header); } #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,y_obs,col,node,z_level) shared(n_rows,GMT,G,flat_earth,Ctrl,n_columns,eval,x_obs,n_prisms,prism,G0) +#pragma omp parallel for private(row,y_obs,col,node,z_level) shared(n_rows,scl_xy,GMT,G,flat_earth,Ctrl,n_columns,eval,x_obs,scl_z,n_prisms,prism,G0) #endif for (row = 0; row < n_rows; row++) { /* Do row-by-row and report on progress if -V */ - y_obs = gmt_M_grd_row_to_y (GMT, row, G->header); - if (!(flat_earth || Ctrl->M.active[GRAVPRISMS_HOR])) y_obs /= METERS_IN_A_KM; /* Convert to km */ + y_obs = scl_xy * gmt_M_grd_row_to_y (GMT, row, G->header); #ifndef _OPENMP GMT_Report (API, GMT_MSG_INFORMATION, "Finished row %5d\n", row); #endif @@ -902,20 +933,23 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { /* Loop over cols; always save the next level before we update the array at that col */ node = gmt_M_ijp (G->header, row, col); z_level = (Ctrl->Z.mode == 1) ? G->data[node] : Ctrl->Z.level; /* Default observation z level unless provided in input grid */ - G->data[node] = (gmt_grdfloat) eval (x_obs[col], y_obs, z_level, n_prisms, prism, G0); + G->data[node] = (gmt_grdfloat) eval (x_obs[col], y_obs, z_level * scl_z, n_prisms, prism, G0); } } gmt_M_free (GMT, x_obs); GMT_Report (API, GMT_MSG_INFORMATION, "Create %s\n", Ctrl->G.file); sprintf (remark, "Calculated 3-D %s", kind[Ctrl->F.mode]); if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_REMARK, remark, G)) { - Return (API->error); + error = API->error; + goto end_it_all; } if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G)) { - Return (API->error); + error = API->error; + goto end_it_all; } if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, G) != GMT_NOERROR) { - Return (API->error); + error = API->error; + goto end_it_all; } } From e2d578eff9b1b5d1579ba48130c777999af0ae61 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 22:21:50 +0000 Subject: [PATCH 10/27] Update prism file format --- .../supplements/potential/gravprisms.rst | 9 ++-- src/potential/gravprisms.c | 44 ++++++++++--------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index b51834a364b..bd6d657b377 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -16,7 +16,7 @@ Synopsis [ |-A| ] [ |-C|\ [**+q**][**+w**\ *file*][**+z**\ *dz*] ] [ |-D|\ *density* ] ] -[ |-E|\ *dx*\ /*dy*\ /*dz* ] +[ |-E|\ *dx*\ /*dy* ] [ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] [ |-G|\ *outfile* ] [ |-H|\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] ] @@ -47,8 +47,9 @@ Description **gravprisms** will compute the geopotential field over vertically oriented, rectangular prisms. We either read the multi-segment *table* from file (or standard input), which may contain up to -7 colums: The first three are the center *x, y, z* coordinates of the prism, while the next -three are the dimensions *dx dy dz* of each prism (see **-E** if all prisms have the same dimesions. +7 columns: The first four are *x y z_low z_high*, i.e., the center *x, y* coordinates and the +vertical range of the prism from *zlow* to *z_high*, while the next two columns hold the dimensions +*dx dy* of each prism (see **-E** if all prisms have the same *x*- and *y*-dimensions. Last column may contain individual prism densities (may be overridden by a fixed density contrast given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate the entire feature (**-S**) or just the volume between two surfaces (one of which may be a constant) @@ -67,7 +68,7 @@ Required Arguments ------------------ *table* - The file describing the prisms with record format *x y z* [ *dx dy dz* ] [ *rho* ], + The file describing the prisms with record format *x y z_lo z_hi* [ *dx dy* ] [ *rho* ], where the optional items are controlled by options **-E** and **-D**, respectively. Any density contrast can be given in kg/m^3 of g/cm^3. diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 8826fe78558..bdf6d03229b 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -76,9 +76,9 @@ struct GRAVPRISMS_CTRL { char *file; double dz; } C; - struct GRAVPRISMS_E { /* -E/ fixed prism dimensions [read from file] */ + struct GRAVPRISMS_E { /* -E fixed prism x/y dimensions [read from file] */ bool active; - double dx, dy, dz; + double dx, dy; } E; struct GRAVPRISMS_F { /* -F[f|n[]|v] */ bool active, lset; @@ -139,6 +139,8 @@ static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new /* Initialize values whose defaults are not 0/false/NULL */ C->F.lat = 45.0; /* So we compute normal gravity at 45 */ + C->H.p = 1.0; /* Linear density increase */ + C->H.densify = 0.0; /* No water-driven compaction on flanks */ return (C); } @@ -202,10 +204,10 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->D.active = true; Ctrl->D.rho = atof (opt->arg); break; - case 'E': /* Set fixed prism parameters instead of reading from prismfile */ + case 'E': /* Set fixed prism dx, dy parameters instead of reading from prismfile */ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); Ctrl->E.active = true; - sscanf (opt->arg, "%lg/%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy, &Ctrl->E.dz); + sscanf (opt->arg, "%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy); break; case 'F': /* Seldct geopotential type */ n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active); @@ -362,9 +364,9 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-C[+q][+w][+z]] [-D] [-Ff|n[]|v] [-G] " - "[-H//[+d][+p]] [-L] [%s] [-M[hz]] [-N] [%s] [-S] " - "[-T] [%s] [-W] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", + GMT_Usage (API, 0, "usage: %s [-A] [-C[+q][+w][+z]] [-D] [-D/ [-Ff|n[]|v] " + "[-G] [-H//[+d][+p]] [-L] [%s] [-M[hz]] [-N] [%s] " + "[-S] [-T] [%s] [-W] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_x_OPT, GMT_PAR_OPT); @@ -383,6 +385,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "+z Set increment for discretization of rho(r,z) when -H is used."); GMT_Usage (API, 1, "\n-D"); GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column or computes it via -H]."); + GMT_Usage (API, 1, "\n-E/"); + GMT_Usage (API, -2, "Set fixed x- and y-dimensions for all prisms [Default reads it from columns 4-5]."); GMT_Usage (API, 1, "\n-Ff|n[]|v]"); GMT_Usage (API, -2, "Specify desired geopotential field component:"); GMT_Usage (API, 3, "f: Free-air anomalies (mGal) [Default]."); @@ -580,7 +584,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { bool flat_earth = false; char *uname[2] = {"meter", "km"}, *kind[3] = {"FAA", "VGG", "GEOID"}, remark[GMT_LEN64] = {""}; - double z_level, rho, dx, dy, dz, lat, G0, scl_xy, scl_z, i_scl_xy, i_scl_z, *prism[7]; + double z_level, rho, dx, dy, lat, G0, scl_xy, scl_z, i_scl_xy, i_scl_z, *prism[7]; double (*eval) (double, double, double, uint64_t, double **, double); struct GRAVPRISMS_CTRL *Ctrl = NULL; @@ -708,10 +712,10 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { for (k = 0; k < n_prisms; k++) { /* Unscramble edges to coordinates and dimensions and possibly convert back to km */ S->data[0][k] = 0.5 * (prism[1][k] + prism[0][k]) * i_scl_xy; /* Get x */ S->data[1][k] = 0.5 * (prism[3][k] + prism[2][k]) * i_scl_xy; /* Get y */ - S->data[2][k] = 0.5 * (prism[5][k] + prism[4][k]) * i_scl_xy; /* Get z */ - S->data[3][k] = (prism[1][k] - prism[0][k]) * i_scl_xy; /* Get dx */ - S->data[4][k] = (prism[3][k] - prism[2][k]) * i_scl_xy; /* Get dy */ - S->data[5][k] = (prism[5][k] - prism[4][k]) * i_scl_xy; /* Get dz */ + S->data[2][k] = prism[4][k] * i_scl_xy; /* Get z_low */ + S->data[3][k] = prism[5][k] * i_scl_xy; /* Get z_high */ + S->data[4][k] = (prism[1][k] - prism[0][k]) * i_scl_xy; /* Get dx */ + S->data[5][k] = (prism[3][k] - prism[2][k]) * i_scl_xy; /* Get dy */ } gmt_M_memcpy (S->data[6], prism[6], n_prisms, double); /* Copy over densities */ if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_SET, NULL, Ctrl->C.file, P) != GMT_NOERROR) { @@ -743,17 +747,17 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { prism[k] = gmt_M_memory (GMT, NULL, P->n_records, double); /* Fill in prism arrays from data set, including optional choices for fixed or variable dimension and density, then free the dataset */ - dx = 0.5 * Ctrl->E.dx; dy = 0.5 * Ctrl->E.dy; dz = 0.5 * Ctrl->E.dz; /* Distances from prism center to respective edges */ + dx = 0.5 * Ctrl->E.dx; dy = 0.5 * Ctrl->E.dy; /* Distances from prism center to respective edges */ rho = Ctrl->D.rho; - col = (Ctrl->E.active) ? 3 : 6; /* Input column for density, if -D is not set */ + col = (Ctrl->E.active) ? 4 : 6; /* Input column for density, if -D is not set */ for (tbl = k = 0; tbl < P->n_tables; tbl++) { for (seg = 0; seg < P->table[tbl]->n_segments; seg++) { S = P->table[tbl]->segment[seg]; for (row = 0; row < S->n_rows; row++, k++) { - if (!Ctrl->E.active) { /* Update the half-dimensions of the sides */ - dx = 0.5 * S->data[3][row]; - dy = 0.5 * S->data[4][row]; - dz = 0.5 * S->data[5][row]; + /* x y z_lo z_hi [dx dy] [rho] */ + if (!Ctrl->E.active) { /* Update the half-dimensions of the x/y-sides */ + dx = 0.5 * S->data[4][row]; + dy = 0.5 * S->data[5][row]; } if (!Ctrl->D.active) rho = S->data[col][row]; /* Here we ensure prism is in meters */ @@ -761,8 +765,8 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { prism[1][k] = (S->data[GMT_X][row] + dx) * scl_xy; prism[2][k] = (S->data[GMT_Y][row] - dy) * scl_xy; /* Store the y-coordinates of the y-edges of prism */ prism[3][k] = (S->data[GMT_Y][row] + dy) * scl_xy; - prism[4][k] = (S->data[GMT_Z][row] - dz) * scl_z; /* Store the z-coordinates of the z-edges of prism */ - prism[5][k] = (S->data[GMT_Z][row] + dz) * scl_z; + prism[4][k] = S->data[2][row] * scl_z; /* Store the z-coordinates of the z-edges of prism */ + prism[5][k] = S->data[3][row] * scl_z; prism[6][k] = rho; /* Store the fixed or variable density of prism */ } } From ef670a8dad10aff3db05083842fff275f1d7e128 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Sun, 13 Mar 2022 22:35:00 +0000 Subject: [PATCH 11/27] Update gmt_api.c --- src/gmt_api.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/gmt_api.c b/src/gmt_api.c index c57c1caabd6..216adcb9797 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -13138,6 +13138,10 @@ struct GMT_RESOURCE * GMT_Encode_Options (void *V_API, const char *module_name, deactivate_output = true; /* Turn off implicit output since none is in effect, only secondary -D output */ type = (GMT_Find_Option (API, 'I', *head)) ? 'I' : 'G'; /* Giving -I means we are reading an image */ } + /* 1y. Check if this is the gravprisms module, where primary dataset input should be turned off if -C is used */ + else if (!strncmp (module, "gravprisms", 10U) && (opt = GMT_Find_Option (API, 'C', *head))) { + deactivate_input = true; /* Turn off implicit input since none is in effect */ + } /* 2a. Get the option key array for this module */ key = gmtapi_process_keys (API, keys, type, *head, n_per_family, &n_keys); /* This is the array of keys for this module, e.g., " Date: Sun, 13 Mar 2022 22:50:02 +0000 Subject: [PATCH 12/27] Fix docs --- .../supplements/potential/gravprisms.rst | 2 +- src/potential/gravprisms.c | 20 +++++++++++-------- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index bd6d657b377..3592b1bbbf2 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -248,7 +248,7 @@ References ---------- -Grant, F. S. and West, G. F., 1965, *Interpretation theory in applied geophysics*, +Grant, F. S. and West, G. F., 1965, *Interpretation Theory in Applied Geophysics*, 583 pp., McGraw-Hill. Kim, S.-S., and P. Wessel, 2016, New analytic solutions for modeling vertical diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index bdf6d03229b..096ff2796fd 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -22,8 +22,12 @@ * Calculates geopotential anomalies due to any number of individual * vertical prisms that may have their own unique dimensions (or all * are constant) and individual density contrasts (or all the same). - * For vertically varying density we break prisms into a stack of - * short sub-prisms with constant density. + * Instead of reading in prisms we can create prisms based on grid + * limits (e.g., top and base surfaces) and then use these prisms to + * compute the geopotential anomalies. + * For vertically varying density we break the prisms into a stack of + * short sub-prisms with constant density that approximates the actual + * continuously-varying densities. * * Accelerated with OpenMP; see -x. */ @@ -375,13 +379,13 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n"); GMT_Usage (API, 1, "\n"); GMT_Usage (API, -2, "One or more multiple-segment ASCII data files. If no files are given, standard " - "input is read. Contains (x,y,z) center coordinates of prisms with optional [ dx dy dz] [rho] if not set via -D and ."); + "input is read. Contains (x,y,z_lo,z_hi), i.e., center coordinates of prisms and z-range, with optional [ dx dy] [rho] if not set via -E and -D."); GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n"); GMT_Usage (API, 1, "\n-A The z-axis is positive upwards [Default is positive down]."); GMT_Usage (API, 1, "\n-C[+q][+w][+z]"); GMT_Usage (API, -2, "Create prisms from the (-L) level to (-T) level, or for the full seamount (-S). No will be read. Modifiers are:"); GMT_Usage (API, 3, "+q Quit execution once prism file has been written (see +w). No geopotential calculations are preformed."); - GMT_Usage (API, 3, "+w Write the created prisms to "); + GMT_Usage (API, 3, "+w Write the created prisms to ."); GMT_Usage (API, 3, "+z Set increment for discretization of rho(r,z) when -H is used."); GMT_Usage (API, 1, "\n-D"); GMT_Usage (API, -2, "Set fixed density contrast (in kg/m^3) [Default reads it from last numerical column or computes it via -H]."); @@ -402,7 +406,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "+p Exponential coefficient (> 0) for density change with burial depth [1 (linear)]."); GMT_Option (API, "I"); GMT_Usage (API, 1, "\n-L"); - GMT_Usage (API, -2, "Set the lower (base) surface grid or constant of a layer to create prisms for; requires -C, -S, and -T [0]"); + GMT_Usage (API, -2, "Set the lower (base) surface grid or constant of a layer to create prisms for; requires -C and -T [0]"); GMT_Usage (API, 1, "\n-M[hz]"); GMT_Usage (API, -2, "Change distance units used, via one or two directives:"); GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); @@ -416,7 +420,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, -2, "Set the full surface grid height for making prisms of the entire feature (or if -H is set) (or use -L and -T to select a layer); requires -C"); GMT_Option (API, "V"); GMT_Usage (API, 1, "\n-T"); - GMT_Usage (API, -2, "Set the top surface grid or constant of a layer to create prisms for; requires -C, -L and -T"); + GMT_Usage (API, -2, "Set the top surface grid or constant of a layer to create prisms for; requires -C and -L"); GMT_Usage (API, 1, "\n-W"); GMT_Usage (API, -2, "Give grid with vertically-averaged spatially varying densities; requires co-registered grids with -C."); GMT_Usage (API, 1, "\n-Z"); @@ -469,7 +473,7 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou } GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { - /* Gravity anomaly from a single prism */ + /* Gravity anomaly from a single prism [Grant & West, 1965] */ double g, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double g111, g112, g121, g122, g211, g212, g221, g222; @@ -504,7 +508,7 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub } GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { - /* Vertical gravity gradient from a single prism [Kim & Wessel, 2016] */ + /* Vertical gravity gradient anomaly from a single prism [Kim & Wessel, 2016] */ double v, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double v111, v112, v121, v122, v211, v212, v221, v222; From 33e09a737d053b60753955627f6b62fdcb1b5705 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 09:47:14 +0000 Subject: [PATCH 13/27] Add test/doc script --- .../supplements/potential/gravprisms.rst | 9 +++++- doc/scripts/GMT_seamount_prisms.sh | 31 +++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) create mode 100755 doc/scripts/GMT_seamount_prisms.sh diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 3592b1bbbf2..ef75bafc328 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -63,6 +63,13 @@ output points specified via **-N**. Choose between free-air anomalies, vertical gravity gradient anomalies, or geoid anomalies. Options are available to control axes units and direction. +.. figure:: /_images/GMT_seamount_prisms.* + :width: 500 px + :align: center + + Three density models explored for a truncated Gaussian seamount via **-C**: (left) Constant + density (**-D**), (middle) Vertically-averaged density varying radially (**-W**), and + (right) density varies with *r* and *z* (**-H**). Required Arguments ------------------ @@ -70,7 +77,7 @@ Required Arguments *table* The file describing the prisms with record format *x y z_lo z_hi* [ *dx dy* ] [ *rho* ], where the optional items are controlled by options **-E** and **-D**, respectively. - Any density contrast can be given in kg/m^3 of g/cm^3. + Density contrasts can be given in kg/m^3 of g/cm^3. .. _-I: diff --git a/doc/scripts/GMT_seamount_prisms.sh b/doc/scripts/GMT_seamount_prisms.sh new file mode 100755 index 00000000000..ec5d87fdf5d --- /dev/null +++ b/doc/scripts/GMT_seamount_prisms.sh @@ -0,0 +1,31 @@ +#!/usr/bin/env bash +# Show the capability of gravprisms to generate the correct prisms +# We start with the grid of a Gaussian seamount (Truncated) and generate three +# prism sets: (1) constant density, (2) variable horizontal density, and (3) full +# variable density. + +# 1. Make the seamount, then remove one quadrant: +echo "0 0 30 6" | gmt grdseamount -R-30/30/-30/30 -I1 -r -Ggausssmt.grd -Cg -F0.2 -H7/2400/3030+p0.8 -Wgaussaverho.grd +gmt grdmath gausssmt.grd 1 X 0 LT Y 0 LT MUL SUB MUL = gausssmt.grd +rho_l=$(gmt grdinfo gaussaverho.grd | grep Remark | awk '{print $NF}') +# 2. Make the constant density prisms +gmt gravprisms -Sgausssmt.grd -D${rho_l} -C+wgaussprisms1.txt+q +# 3. Make variable constant prism densities +gmt gravprisms -Sgausssmt.grd -Wgaussaverho.grd -C+wgaussprisms2.txt+q +# 4. Make variable prism densities +gmt gravprisms -Sgausssmt.grd -H7/2400/3030+p0.8 -C+wgaussprisms3.txt+q+z0.2 +# Make plots +gmt begin GMT_seamount_prisms + gmt set FONT_TAG 12p,Times-Italic + gmt makecpt -Cturbo -I -T2400/2950 + gmt subplot begin 1x3 -Fs4.8c/2.75c -A + gmt subplot set 0 -A"@~r@~@-l@-" + gmt plot3d -R-30/30/-30/30/0/6 -JX5c -JZ1.75c -C -So1q+b gaussprisms1.txt -i0:2,6,3 -Baf -Bzaf -BWSrtZ -p200/20 -Wfaint + gmt subplot set 1 -A"@~r@~@-l@-(r)" + gmt plot3d -R-30/30/-30/30/0/6 -JX5c -JZ1.75c -C -So1q+b gaussprisms2.txt -i0:2,6,3 -Baf -BwSrt -p200/20 -Wfaint + gmt subplot set 2 -A"@~r@~@-l@-(r,z)" + gmt plot3d -R-30/30/-30/30/0/6 -JX5c -JZ1.75c -C -So1q+b gaussprisms3.txt -i0:2,6,3 -Baf -BwSrt -p200/20 -Wfaint + gmt subplot end + gmt colorbar -DJRM+o1.5c/0.5c -Bxaf -By+l"kg\267m@+-3@+" +gmt end show +rm -f gaussprisms?.txt gaussaverho.grd gausssmt.grd From 778bf855845dbe1f76b7c5405302cac7687ec2d2 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 09:48:44 +0000 Subject: [PATCH 14/27] Update images.dvc --- doc/scripts/images.dvc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/scripts/images.dvc b/doc/scripts/images.dvc index 06ffc05c23e..c5603260b2e 100644 --- a/doc/scripts/images.dvc +++ b/doc/scripts/images.dvc @@ -1,5 +1,5 @@ outs: -- md5: 1bff23b2d513c7b9c09693eafac399c1.dir - size: 18545710 - nfiles: 193 +- md5: 34a9efcc6aaabecfb74b746b0fbc0c5f.dir + size: 31115765 + nfiles: 194 path: images From b66439847c288f84f403a67bf73fc6c35151612e Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 10:08:31 +0000 Subject: [PATCH 15/27] Doc updates --- doc/rst/source/supplements/potential/gravprisms.rst | 3 ++- src/potential/gravprisms.c | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index ef75bafc328..211e775cbac 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -162,8 +162,9 @@ Optional Arguments **-N**\ *trackfile* Specifies individual (x, y[, z]) locations where we wish to compute the predicted value. When this option - is used there are no grids and the output data records are written to standard output (see **-bo** for binary output). + is used there are no grids involved and the output data records are written to standard output (see **-bo** for binary output). If *trackfile* has 3 columns we take the *z* value as our observation level; this level may be overridden via **-Z**. + **Note**: If **-G** is used to set an output file we will write the output table to that file instead of standard output. .. _-S: diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 096ff2796fd..d6ae7eacd12 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -805,7 +805,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (gmt_M_is_geographic (GMT, GMT_IN)) lat = 0.5 * (G->header->wesn[YLO] + G->header->wesn[YHI]); /* Mid-latitude needed if geoid is selected */ } else { /* Got a dataset with output locations via -N */ - unsigned int n_expected = 4; /* Max number of columns */ + unsigned int n_expected = 3; /* Max number of columns is 3 (x, y, z) but if -Z is set then just (x, y) */ if (Ctrl->Z.active) n_expected--; gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading from the -N file */ if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { @@ -816,8 +816,8 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { error = API->error; goto end_it_all; } - if (D->n_columns < 2) { - GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->N.file, (int)D->n_columns); + if (D->n_columns < n_expected) { + GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but %d are needed\n", Ctrl->N.file, (int)D->n_columns, n_expected); error = GMT_DIM_TOO_SMALL; goto end_it_all; } From 19ed61e6369c21d8e2e5a37a9f12a378b3f9a915 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 10:40:19 +0000 Subject: [PATCH 16/27] Update gravprisms.c --- src/potential/gravprisms.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index d6ae7eacd12..b42d40b8643 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -221,7 +221,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'n': Ctrl->F.mode = GRAVPRISMS_GEOID; if (opt->arg[1]) Ctrl->F.lat = atof (&opt->arg[1]), Ctrl->F.lset = true; break; - case 'g': Ctrl->F.mode = GRAVPRISMS_FAA; break; + case 'f': Ctrl->F.mode = GRAVPRISMS_FAA; break; default: GMT_Report (API, GMT_MSG_WARNING, "Option -F: Unrecognized field %c\n", opt->arg[0]); n_errors++; From 7698c4a8d0f708d894f48fdad7e956fa129f93fd Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 15:03:29 +0000 Subject: [PATCH 17/27] Add geo case --- src/potential/gravprisms.c | 77 +++++++++++++++++++++++++++++--------- src/potential/talwani.h | 3 +- 2 files changed, 61 insertions(+), 19 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index b42d40b8643..287ec2226ff 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -43,19 +43,10 @@ #define THIS_MODULE_NEEDS "r" #define THIS_MODULE_OPTIONS "-VRbdefhior" GMT_ADD_x_OPT -/* Geographic coordinates not implemented yet */ +/* Note: All calculations assume distances in meters so degrees -> meters and km -> meters first */ -#define DX_FROM_DLON(x1, x2, y1, y2) (((x1) - (x2)) * DEG_TO_KM * cos (0.5 * ((y1) + (y2)) * D2R)) -#define DY_FROM_DLAT(y1, y2) (((y1) - (y2)) * DEG_TO_KM) -#define gmt_M_is_zero2(x) (fabs (x) < 2e-5) /* Check for near-zero angles [used in geoid integral()]*/ - -#define GMT_GRAVPRISMS_N_DEPTHS GMT_BUFSIZ /* Max depths allowed due to OpenMP needing stack array */ - -#if 0 -/* this is for various debug purposes and will eventually be purged */ -bool dump = false; -FILE *fp = NULL; -#endif +#define DX_FROM_DLON(x1, x2, y1, y2) (((x1) - (x2)) * DEG_TO_M * cos (0.5 * ((y1) + (y2)) * D2R)) +#define DY_FROM_DLAT(y1, y2) (((y1) - (y2)) * DEG_TO_M) enum gravprisms_fields { GRAVPRISMS_FAA = 0, @@ -89,7 +80,7 @@ struct GRAVPRISMS_CTRL { unsigned int mode; double lat; } F; - struct GRAVPRISMS_G { /* Output file */ + struct GRAVPRISMS_G { /* Output grid or profile */ bool active; char *file; } G; @@ -551,6 +542,23 @@ GMT_LOCAL double gravprisms_get_one_g_output (double x, double y, double z, uint return (g); } +GMT_LOCAL double gravprisms_get_one_g_output_geo (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { + /* (x, y, z) is the observation point */ + double g = 0.0; + double dx1, dx2, dy1, dy2, ym; + gmt_M_unused (unused); + for (uint64_t k = 0; k < n_prisms; k++) { + /* Got lon, lat and must convert to Flat-Earth km */ + ym = 0.5 * (P[2][k] + P[3][k]); /* Prism mid-y-value */ + dx1 = DX_FROM_DLON (P[0][k], x, ym, y); + dx2 = DX_FROM_DLON (P[1][k], x, ym, y); + dy1 = DY_FROM_DLAT (P[2][k], y); + dy2 = DY_FROM_DLAT (P[3][k], y); + g += gravprism (dx1, dx2, dy1, dy2, P[4][k]-z, P[5][k]-z, P[6][k]); + } + return (g); +} + GMT_LOCAL double gravprisms_get_one_n_output (double x, double y, double z, uint64_t n_prisms, double **P, double constant) { /* (x, y, z) is the observation point */ double n = 0.0; @@ -559,6 +567,22 @@ GMT_LOCAL double gravprisms_get_one_n_output (double x, double y, double z, uint return (n * constant * 0.01); /* To get geoid in meter */ } +GMT_LOCAL double gravprisms_get_one_n_output_geo (double x, double y, double z, uint64_t n_prisms, double **P, double constant) { + /* (x, y, z) is the observation point */ + double n = 0.0; + double dx1, dx2, dy1, dy2, ym; + for (uint64_t k = 0; k < n_prisms; k++) { + /* Got lon, lat and must convert to Flat-Earth km */ + ym = 0.5 * (P[2][k] + P[3][k]); /* Prism mid-y-value */ + dx1 = DX_FROM_DLON (P[0][k], x, ym, y); + dx2 = DX_FROM_DLON (P[1][k], x, ym, y); + dy1 = DY_FROM_DLAT (P[2][k], y); + dy2 = DY_FROM_DLAT (P[3][k], y); + n += geoidprism (dx1, dx2, dy1, dy2, P[4][k]-z, P[5][k]-z, P[6][k]); + } + return (n * constant * 0.01); /* To get geoid in meter */ +} + GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { /* (x, y, z) is the observation point */ double v = 0.0; @@ -568,6 +592,23 @@ GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ } +GMT_LOCAL double gravprisms_get_one_v_output_geo (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { + /* (x, y, z) is the observation point */ + double v = 0.0; + double dx1, dx2, dy1, dy2, ym; + gmt_M_unused (unused); + for (uint64_t k = 0; k < n_prisms; k++) { + /* Got lon, lat and must convert to Flat-Earth km */ + ym = 0.5 * (P[2][k] + P[3][k]); /* Prism mid-y-value */ + dx1 = DX_FROM_DLON (P[0][k], x, ym, y); + dx2 = DX_FROM_DLON (P[1][k], x, ym, y); + dy1 = DY_FROM_DLAT (P[2][k], y); + dy2 = DY_FROM_DLAT (P[3][k], y); + v += vggprism (dx1, dx2, dy1, dy2, P[4][k]-z, P[5][k]-z, P[6][k]); + } + return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ +} + GMT_LOCAL double gravprisms_rho (struct GRAVPRISMS_CTRL *Ctrl, double z, double h) { /* Passing in the current seamount's normalized height, h(r) and the current normalized radial point z(r). * We return the density at this point inside the seamount with reference height 1 */ @@ -845,15 +886,15 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { switch (Ctrl->F.mode) { /* Set pointer to chosen geopotential evaluation function */ case GRAVPRISMS_VGG: - eval = &gravprisms_get_one_v_output; + eval = (flat_earth) ? &gravprisms_get_one_v_output_geo : &gravprisms_get_one_v_output; break; case GRAVPRISMS_GEOID: - eval = &gravprisms_get_one_n_output; + eval = (flat_earth) ? &gravprisms_get_one_n_output_geo : &gravprisms_get_one_n_output; G0 = (Ctrl->F.lset) ? g_normal (Ctrl->F.lat) : g_normal (lat); G0 = 1.0 / G0; /* To avoid dividing in the loop */ break; case GRAVPRISMS_FAA: - eval = &gravprisms_get_one_g_output; + eval = (flat_earth) ? &gravprisms_get_one_g_output_geo : &gravprisms_get_one_g_output; break; default: /* Just for Coverity */ @@ -897,7 +938,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { gmt_prep_tmp_arrays (GMT, GMT_OUT, S->n_rows, 1); /* Init or reallocate tmp vector */ #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl_xy,scl_z,n_prisms,prism,flat_earth,G0) +#pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl_xy,scl_z,n_prisms,prism,G0) #endif /* Separate the calculation from the output in two separate row-loops since cannot do rec-by-rec output * with OpenMP due to race condiations that would mess up the output order */ @@ -930,7 +971,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { } #ifdef _OPENMP /* Spread calculation over selected cores */ -#pragma omp parallel for private(row,y_obs,col,node,z_level) shared(n_rows,scl_xy,GMT,G,flat_earth,Ctrl,n_columns,eval,x_obs,scl_z,n_prisms,prism,G0) +#pragma omp parallel for private(row,y_obs,col,node,z_level) shared(n_rows,scl_xy,GMT,G,Ctrl,n_columns,eval,x_obs,scl_z,n_prisms,prism,G0) #endif for (row = 0; row < n_rows; row++) { /* Do row-by-row and report on progress if -V */ y_obs = scl_xy * gmt_M_grd_row_to_y (GMT, row, G->header); diff --git a/src/potential/talwani.h b/src/potential/talwani.h index 6bd7755d5d7..434f76eea9a 100644 --- a/src/potential/talwani.h +++ b/src/potential/talwani.h @@ -26,7 +26,8 @@ #define TOL 1.0e-7 /* Gotta leave a bit slack for these calculations */ #define SI_TO_MGAL 1.0e5 /* Convert m/s^2 to mGal */ #define SI_TO_EOTVOS 1.0e9 /* Convert (m/s^2)/m to Eotvos */ -#define DEG_TO_KM 111.319490793 /* For flat-Earth scaling of degrees to km on WGS-84 Equator */ +#define DEG_TO_KM 111.319490793 /* For flat-Earth scaling of degrees to km on WGS-84 Equator */ +#define DEG_TO_M 111319.490793 /* For flat-Earth scaling of degrees to m on WGS-84 Equator */ #define SI_GAMMA 6.673e-11 /* Gravitational constant (SI units) */ #define GAMMA 6.673 /* Gravitational constant for distances in km and mass in kg/m^3 */ From b1013fb10b5da8b4545e40e4398c66cb643dc4c8 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 21:46:49 +0000 Subject: [PATCH 18/27] Update gravprisms.rst --- .../supplements/potential/gravprisms.rst | 41 ++++++++++--------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 211e775cbac..76c1761c614 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -15,7 +15,7 @@ Synopsis **gmt gravprisms** [ *table* ] [ |-A| ] [ |-C|\ [**+q**][**+w**\ *file*][**+z**\ *dz*] ] -[ |-D|\ *density* ] ] +[ |-D|\ *density* ] [ |-E|\ *dx*\ /*dy* ] [ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] [ |-G|\ *outfile* ] @@ -48,10 +48,10 @@ Description **gravprisms** will compute the geopotential field over vertically oriented, rectangular prisms. We either read the multi-segment *table* from file (or standard input), which may contain up to 7 columns: The first four are *x y z_low z_high*, i.e., the center *x, y* coordinates and the -vertical range of the prism from *zlow* to *z_high*, while the next two columns hold the dimensions -*dx dy* of each prism (see **-E** if all prisms have the same *x*- and *y*-dimensions. -Last column may contain individual prism densities (may be overridden by a fixed density contrast -given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate +vertical range of the prism from *z_low* to *z_high*, while the next two columns hold the dimensions +*dx dy* of each prism (see **-E** if all prisms have the same *x*- and *y*-dimensions). +Last column may contain individual prism densities (but will be overridden by a fixed density contrast +if given via **-D**). Alternatively, we can use **-C** to create the prisms needed to approximate the entire feature (**-S**) or just the volume between two surfaces (one of which may be a constant) that define a layer (set via **-L** and **-T**). If a variable density model (**-H**) is selected then each vertical prism will be broken into constant-density, stacked sub-prisms using a prescribed @@ -67,9 +67,9 @@ axes units and direction. :width: 500 px :align: center - Three density models explored for a truncated Gaussian seamount via **-C**: (left) Constant - density (**-D**), (middle) Vertically-averaged density varying radially (**-W**), and - (right) density varies with *r* and *z* (**-H**). + Three density models modeled by prisms for a truncated Gaussian seamount via **-C**: (left) Constant + density (**-D**), (middle) vertically-averaged density varying radially (**-W**), and + (right) density varies with *r* and *z* (**-H**), requiring a stack of prisms. Required Arguments ------------------ @@ -94,8 +94,7 @@ Optional Arguments .. _-A: **-A** - The *z*-axis should be positive upwards [Default is down]. All internal *z*-values will have their sign - changed. **Note**: Output *z*-values are not affected. + The *z*-axis should be positive upwards [Default is down]. .. _-C: @@ -117,13 +116,13 @@ Optional Arguments .. _-E: -**-E**\ *dx*\ /*dy*\ /*dz* - If all prisms in *table* have constant dimensions then they can be set here. In that case *table* - only contains the centers of each prism (and optionally *density*; see **-D**). +**-E**\ *dx*\ /*dy* + If all prisms in *table* have constant x/y-dimensions then they can be set here. In that case *table* + only contains the centers of each prism and the *z* range (and optionally *density*; see **-D**). .. _-F: -**-F**\ **f**\|\ **n**\|\ **v** +**-F**\ **f**\|\ **n**\ [*lat*]\|\ **v** Specify desired gravitational field component. Choose between **f** (free-air anomaly) [Default], **n** (geoid; optionally append average latitude for normal gravity reference value [Default is mid-grid (or mid-profile if **-N**)]) or **v** (vertical gravity gradient). @@ -133,7 +132,8 @@ Optional Arguments **-G**\ *outfile* Specify the name of the output data (for grids, see :ref:`Grid File Formats `). Required when an equidistant grid is implied for output. - If **-N** is used then output is written to standard output unless **-G** specifies an output file. + If **-N** is used then output is written to standard output unless **-G** + specifies an output file. .. _-H: @@ -148,7 +148,8 @@ Optional Arguments .. _-L: **-L**\ *base* - Give name of the base grid surface for a layer we wish to fill with prisms, or give a constant *z*-level [0]. + Give name of the base surface grid for a layer we wish to approximate with prisms, or give + a constant *z*-level [0]. .. _-M: @@ -163,7 +164,7 @@ Optional Arguments **-N**\ *trackfile* Specifies individual (x, y[, z]) locations where we wish to compute the predicted value. When this option is used there are no grids involved and the output data records are written to standard output (see **-bo** for binary output). - If *trackfile* has 3 columns we take the *z* value as our observation level; this level may be overridden via **-Z**. + If **-Z** is not set then *trackfile* must have 3 columns and we take the *z* value as our observation level; otherwise the level must be set via **-Z**. **Note**: If **-G** is used to set an output file we will write the output table to that file instead of standard output. .. _-S: @@ -174,7 +175,7 @@ Optional Arguments .. _-T: **-T**\ *top* - Give name of the top surface grid for a layer we wish to fill with prisms, or give a constant *z*-level. + Give name of the top surface grid for a layer we wish to approximate with prisms, or give a constant *z*-level. .. |Add_-V| replace:: |Add_-V_links| .. include:: /explain_-V.rst_ @@ -184,8 +185,8 @@ Optional Arguments .. _-W: **-W**\ *avedens* - Give name of a grid with spatially varying, vertically-averaged prism densities. Requires **-C** and - the grid must be co-registered with the grid provided by **-S**. + Give name of an input grid with spatially varying, vertically-averaged prism densities. Requires **-C** and + the grid must be co-registered with the grid provided by **-S** (or **L** and **-T**). .. _-Z: From 35f240cd19d536f9173614bad126d7f716a095e8 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Mon, 14 Mar 2022 21:48:56 +0000 Subject: [PATCH 19/27] Update gravprisms.c --- src/potential/gravprisms.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 287ec2226ff..5b58cbbe8b6 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -309,13 +309,13 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'Z': /* Observation level(s) */ n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active); Ctrl->Z.active = true; - if (!gmt_access (GMT, opt->arg, F_OK)) { /* File with z-levels exists */ + if (opt->arg[0] && !gmt_access (GMT, opt->arg, F_OK)) { /* File with z-levels exists */ Ctrl->Z.file = strdup (opt->arg); Ctrl->Z.mode = 1; } else { /* Got a constant z-level */ Ctrl->Z.mode = 0; - Ctrl->Z.level = atof (opt->arg); + Ctrl->Z.level = (opt->arg[0]) ? atof (opt->arg) : 0.0; } break; default: From be6ca76cf59a25e305f806dad82161700355bfab Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 15 Mar 2022 09:18:28 +0000 Subject: [PATCH 20/27] Add numerical test for gravity --- .../supplements/potential/gravprisms.rst | 6 ++++-- src/potential/gravprisms.c | 10 ++++++---- test/potential/prism_check_faa.sh | 20 +++++++++++++++++++ 3 files changed, 30 insertions(+), 6 deletions(-) create mode 100755 test/potential/prism_check_faa.sh diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 76c1761c614..326d74df4f5 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -16,7 +16,7 @@ Synopsis [ |-A| ] [ |-C|\ [**+q**][**+w**\ *file*][**+z**\ *dz*] ] [ |-D|\ *density* ] -[ |-E|\ *dx*\ /*dy* ] +[ |-E|\ *dx*\ [/*dy*] ] [ |-F|\ **f**\|\ **n**\ [*lat*]\|\ **v** ] [ |-G|\ *outfile* ] [ |-H|\ *H*/*rho_l*/*rho_h*\ [**+d**\ *densify*][**+p**\ *power*] ] @@ -116,9 +116,11 @@ Optional Arguments .. _-E: -**-E**\ *dx*\ /*dy* +**-E**\ *dx*\ [/*dy*] If all prisms in *table* have constant x/y-dimensions then they can be set here. In that case *table* only contains the centers of each prism and the *z* range (and optionally *density*; see **-D**). + If only *dx* is given then we set *dy = dx*. **Note**: For geographic coordinates the *dx* dimension is + in geographic delta longitude and hence the physical width of the prism will decrease with latitude if *dx* stays the same. .. _-F: diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 5b58cbbe8b6..a0574fc83c9 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -71,7 +71,7 @@ struct GRAVPRISMS_CTRL { char *file; double dz; } C; - struct GRAVPRISMS_E { /* -E fixed prism x/y dimensions [read from file] */ + struct GRAVPRISMS_E { /* -E[/] fixed prism x/y dimensions [read from file] */ bool active; double dx, dy; } E; @@ -154,6 +154,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *C) { /* Dea static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT_OPTION *options) { unsigned int k, n_errors = 0; + int ns; char *c = NULL; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; @@ -202,7 +203,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT case 'E': /* Set fixed prism dx, dy parameters instead of reading from prismfile */ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); Ctrl->E.active = true; - sscanf (opt->arg, "%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy); + ns = sscanf (opt->arg, "%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy); + if (ns == 1) Ctrl->E.dy = Ctrl->E.dx; break; case 'F': /* Seldct geopotential type */ n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active); @@ -359,7 +361,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT 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] [-C[+q][+w][+z]] [-D] [-D/ [-Ff|n[]|v] " + GMT_Usage (API, 0, "usage: %s [-A] [-C[+q][+w][+z]] [-D] [-E[/]] [-Ff|n[]|v] " "[-G] [-H//[+d][+p]] [-L] [%s] [-M[hz]] [-N] [%s] " "[-S] [-T] [%s] [-W] [-Z] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, @@ -775,7 +777,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { unsigned int n_expected = 7; /* Max number of columns */ /* Specify input expected columns to be at least 3 */ if (Ctrl->D.active) n_expected--; /* Not reading density from records */ - if (Ctrl->E.active) n_expected -= 3; /* Not reading dx dy dz from records */ + if (Ctrl->E.active) n_expected -= 2; /* Not reading dx dy from records */ if ((error = GMT_Set_Columns (API, GMT_IN, n_expected, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) { Return (error); } diff --git a/test/potential/prism_check_faa.sh b/test/potential/prism_check_faa.sh new file mode 100755 index 00000000000..66c7f2c4058 --- /dev/null +++ b/test/potential/prism_check_faa.sh @@ -0,0 +1,20 @@ +#!/usr/bin/env bash +# Check that the calculation for gravity over a prism is correct +# Results validated by hand with separate Matlab code +# We test both Cartesian (in km) and geographic (degrees) + +cat <<- EOF > answer.txt +10.000000 10.000000 7000.000000 5.186362 +0.100000 0.100000 7000.000000 4.861885 +EOF +# Prism centered at (0,0) ranging from z = 0-5000 m +echo "0 0 0 5000" > p.txt +# Cartesian observation point at (10, 10) km at 7000 m elevation +echo " 10 10 7000" > c.txt +# Geopgraphic observation point at lon/lat (0.1, 0.1) at 7000 m elevation +echo "0.1 0.1 7000" > g.txt +# Cartesian prism has dimensions 10x10 km +gmt gravprisms p.txt -E10 -Nc.txt -Mh -D1000 --FORMAT_FLOAT_OUT=%.6lf > result.txt +# Geographic prism has dimension 0.1 degrees by 0.1 degrees in longitude/latitude +gmt gravprisms p.txt -E0.1 -Ng.txt -fg -D1000 --FORMAT_FLOAT_OUT=%.6lf >> result.txt +diff result.txt answer.txt --strip-trailing-cr > fail From 86e809cf3043f2f5c4cabfc5bb160fb5f995caa7 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 15 Mar 2022 15:16:25 +0000 Subject: [PATCH 21/27] Update gravprisms.c --- src/potential/gravprisms.c | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index a0574fc83c9..0e5ed285eeb 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -427,8 +427,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) { return (GMT_MODULE_USAGE); } -#define GRAVITATIONAL_CONST 6.674e-11 -#define GRAVITATIONAL_CONST_MGAL 6.674e-6 +#define GRAVITATIONAL_CONST_FAA 6.674e-6 /* To convert m/s^2 to mGal requires 1e5 */ +#define GRAVITATIONAL_CONST_GEOID 6.674e-11 /* This gives geoid in meter once we divide by g0 */ +#define GRAVITATIONAL_CONST_VGG 6.674e-2 /* To convert mGal/m to 0.1 mGal/km requires 1e4 */ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { /* Geoid anomaly from a single prism [Nagy et al, 2000] */ @@ -450,17 +451,17 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); /* Evaluate at dz1 */ - n111 = -(0.5 * (dx1_sq * atan ((dy1 * dz1) / (dz1 * R111)) + dy1_sq * atan ((dx1 * dz1) / (dz1 * R111)) + dz1_sq * atan ((dx1 * dy1) / (dz1 * R111))) - dx1 * dz1 * log(R111 + dy1) - dy1 * dz1 * log (R111 + dx1) - dx1 * dy1 * log (R111 + dz1)); - n112 = +(0.5 * (dx2_sq * atan ((dy1 * dz1) / (dz1 * R112)) + dy1_sq * atan ((dx2 * dz1) / (dz1 * R112)) + dz1_sq * atan ((dx2 * dy1) / (dz1 * R112))) - dx2 * dz1 * log(R112 + dy1) - dy1 * dz1 * log (R112 + dx2) - dx2 * dy1 * log (R112 + dz1)); - n121 = +(0.5 * (dx1_sq * atan ((dy2 * dz1) / (dz1 * R121)) + dy2_sq * atan ((dx1 * dz1) / (dz1 * R121)) + dz1_sq * atan ((dx1 * dy2) / (dz1 * R121))) - dx1 * dz1 * log(R121 + dy2) - dy2 * dz1 * log (R121 + dx1) - dx1 * dy2 * log (R121 + dz1)); - n122 = -(0.5 * (dx2_sq * atan ((dy2 * dz1) / (dz1 * R122)) + dy2_sq * atan ((dx2 * dz1) / (dz1 * R122)) + dz1_sq * atan ((dx2 * dy2) / (dz1 * R122))) - dx2 * dz1 * log(R122 + dy2) - dy2 * dz1 * log (R122 + dx2) - dx2 * dy2 * log (R122 + dz1)); + n111 = -(0.5 * (dx1_sq * atan ((dy1 * dz1) / (dz1 * R111)) + dy1_sq * atan ((dx1 * dz1) / (dz1 * R111)) + dz1_sq * atan ((dx1 * dy1) / (dz1 * R111))) - dx1 * dz1 * log (R111 + dy1) - dy1 * dz1 * log (R111 + dx1) - dx1 * dy1 * log (R111 + dz1)); + n112 = +(0.5 * (dx2_sq * atan ((dy1 * dz1) / (dz1 * R112)) + dy1_sq * atan ((dx2 * dz1) / (dz1 * R112)) + dz1_sq * atan ((dx2 * dy1) / (dz1 * R112))) - dx2 * dz1 * log (R112 + dy1) - dy1 * dz1 * log (R112 + dx2) - dx2 * dy1 * log (R112 + dz1)); + n121 = +(0.5 * (dx1_sq * atan ((dy2 * dz1) / (dz1 * R121)) + dy2_sq * atan ((dx1 * dz1) / (dz1 * R121)) + dz1_sq * atan ((dx1 * dy2) / (dz1 * R121))) - dx1 * dz1 * log (R121 + dy2) - dy2 * dz1 * log (R121 + dx1) - dx1 * dy2 * log (R121 + dz1)); + n122 = -(0.5 * (dx2_sq * atan ((dy2 * dz1) / (dz1 * R122)) + dy2_sq * atan ((dx2 * dz1) / (dz1 * R122)) + dz1_sq * atan ((dx2 * dy2) / (dz1 * R122))) - dx2 * dz1 * log (R122 + dy2) - dy2 * dz1 * log (R122 + dx2) - dx2 * dy2 * log (R122 + dz1)); /* Evaluate at dz2 */ - n211 = +(0.5 * (dx1_sq * atan ((dy1 * dz2) / (dz2 * R211)) + dy1_sq * atan ((dx1 * dz2) / (dz2 * R211)) + dz2_sq * atan ((dx1 * dy1) / (dz2 * R211))) - dx1 * dz2 * log (R211 + dy1) - dy1 * dz2 * log (R211 + dx1) - dx1 * dy1 * log (R111 + dz2)); - n212 = -(0.5 * (dx2_sq * atan ((dy1 * dz2) / (dz2 * R212)) + dy1_sq * atan ((dx2 * dz2) / (dz2 * R212)) + dz2_sq * atan ((dx2 * dy1) / (dz2 * R212))) - dx2 * dz2 * log (R212 + dy1) - dy1 * dz2 * log (R212 + dx2) - dx2 * dy1 * log (R112 + dz2)); - n221 = -(0.5 * (dx1_sq * atan ((dy2 * dz2) / (dz2 * R221)) + dy2_sq * atan ((dx1 * dz2) / (dz2 * R221)) + dz2_sq * atan ((dx1 * dy2) / (dz2 * R221))) - dx1 * dz2 * log (R221 + dy2) - dy2 * dz2 * log (R221 + dx1) - dx1 * dy2 * log (R121 + dz2)); - n222 = +(0.5 * (dx2_sq * atan ((dy2 * dz2) / (dz2 * R222)) + dy2_sq * atan ((dx2 * dz2) / (dz2 * R222)) + dz2_sq * atan ((dx2 * dy2) / (dz2 * R222))) - dx2 * dz2 * log (R222 + dy2) - dy2 * dz2 * log (R222 + dx2) - dx2 * dy2 * log (R122 + dz2)); + n211 = +(0.5 * (dx1_sq * atan ((dy1 * dz2) / (dz2 * R211)) + dy1_sq * atan ((dx1 * dz2) / (dz2 * R211)) + dz2_sq * atan ((dx1 * dy1) / (dz2 * R211))) - dx1 * dz2 * log (R211 + dy1) - dy1 * dz2 * log (R211 + dx1) - dx1 * dy1 * log (R211 + dz2)); + n212 = -(0.5 * (dx2_sq * atan ((dy1 * dz2) / (dz2 * R212)) + dy1_sq * atan ((dx2 * dz2) / (dz2 * R212)) + dz2_sq * atan ((dx2 * dy1) / (dz2 * R212))) - dx2 * dz2 * log (R212 + dy1) - dy1 * dz2 * log (R212 + dx2) - dx2 * dy1 * log (R212 + dz2)); + n221 = -(0.5 * (dx1_sq * atan ((dy2 * dz2) / (dz2 * R221)) + dy2_sq * atan ((dx1 * dz2) / (dz2 * R221)) + dz2_sq * atan ((dx1 * dy2) / (dz2 * R221))) - dx1 * dz2 * log (R221 + dy2) - dy2 * dz2 * log (R221 + dx1) - dx1 * dy2 * log (R221 + dz2)); + n222 = +(0.5 * (dx2_sq * atan ((dy2 * dz2) / (dz2 * R222)) + dy2_sq * atan ((dx2 * dz2) / (dz2 * R222)) + dz2_sq * atan ((dx2 * dy2) / (dz2 * R222))) - dx2 * dz2 * log (R222 + dy2) - dy2 * dz2 * log (R222 + dx2) - dx2 * dy2 * log (R222 + dz2)); - n = rho * GRAVITATIONAL_CONST_MGAL * (n111 + n112 + n121 + n122 + n211 + n212 + n221 + n222); + n = -rho * GRAVITATIONAL_CONST_GEOID * (n111 + n112 + n121 + n122 + n211 + n212 + n221 + n222); return (n); } @@ -495,7 +496,7 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub g221 = -(dz2 * atan ((dx1 * dy2) / (dz2 * R221)) - dx1 * log (R221 + dy2) - dy2 * log (R221 + dx1)); g222 = +(dz2 * atan ((dx2 * dy2) / (dz2 * R222)) - dx2 * log (R222 + dy2) - dy2 * log (R222 + dx2)); - g = -rho * GRAVITATIONAL_CONST_MGAL * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); + g = -rho * GRAVITATIONAL_CONST_FAA * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); return (g); } @@ -530,7 +531,7 @@ GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, doubl v221 = -atan ((dx1 * dy2) / (dz2 * R221)); v222 = +atan ((dx2 * dy2) / (dz2 * R222)); - v = rho * GRAVITATIONAL_CONST_MGAL * (v111 + v112 + v121 + v122 + v211 + v212 + v221 + v222); + v = -rho * GRAVITATIONAL_CONST_VGG * (v111 + v112 + v121 + v122 + v211 + v212 + v221 + v222); return (v); } @@ -566,7 +567,7 @@ GMT_LOCAL double gravprisms_get_one_n_output (double x, double y, double z, uint double n = 0.0; for (uint64_t k = 0; k < n_prisms; k++) n += geoidprism (P[0][k]-x, P[1][k]-x, P[2][k]-y, P[3][k]-y, P[4][k]-z, P[5][k]-z, P[6][k]); - return (n * constant * 0.01); /* To get geoid in meter */ + return (n * constant); /* To get geoid in meter */ } GMT_LOCAL double gravprisms_get_one_n_output_geo (double x, double y, double z, uint64_t n_prisms, double **P, double constant) { @@ -582,7 +583,7 @@ GMT_LOCAL double gravprisms_get_one_n_output_geo (double x, double y, double z, dy2 = DY_FROM_DLAT (P[3][k], y); n += geoidprism (dx1, dx2, dy1, dy2, P[4][k]-z, P[5][k]-z, P[6][k]); } - return (n * constant * 0.01); /* To get geoid in meter */ + return (n * constant); /* To get geoid in meter */ } GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { @@ -591,7 +592,7 @@ GMT_LOCAL double gravprisms_get_one_v_output (double x, double y, double z, uint gmt_M_unused (unused); for (uint64_t k = 0; k < n_prisms; k++) v += vggprism (P[0][k]-x, P[1][k]-x, P[2][k]-y, P[3][k]-y, P[4][k]-z, P[5][k]-z, P[6][k]); - return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ + return (v); /* Converted units from mGal/m to Eotvos = 0.1 mGal/km */ } GMT_LOCAL double gravprisms_get_one_v_output_geo (double x, double y, double z, uint64_t n_prisms, double **P, double unused) { @@ -608,7 +609,7 @@ GMT_LOCAL double gravprisms_get_one_v_output_geo (double x, double y, double z, dy2 = DY_FROM_DLAT (P[3][k], y); v += vggprism (dx1, dx2, dy1, dy2, P[4][k]-z, P[5][k]-z, P[6][k]); } - return (v * 0.0001); /* From mGal/m to Eotvos = 0.1 mGal/km */ + return (v); /* Converted units from mGal/m to Eotvos = 0.1 mGal/km */ } GMT_LOCAL double gravprisms_rho (struct GRAVPRISMS_CTRL *Ctrl, double z, double h) { From 3a192107546228f36af547ac5210002766d3ec84 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 15 Mar 2022 15:25:54 +0000 Subject: [PATCH 22/27] Spelling --- .../supplements/potential/gravprisms.rst | 22 +++++++++++-------- src/potential/gravprisms.c | 16 +++++++------- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 326d74df4f5..a956e59b385 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -119,8 +119,9 @@ Optional Arguments **-E**\ *dx*\ [/*dy*] If all prisms in *table* have constant x/y-dimensions then they can be set here. In that case *table* only contains the centers of each prism and the *z* range (and optionally *density*; see **-D**). - If only *dx* is given then we set *dy = dx*. **Note**: For geographic coordinates the *dx* dimension is - in geographic delta longitude and hence the physical width of the prism will decrease with latitude if *dx* stays the same. + If only *dx* is given then we set *dy = dx*. **Note**: For geographic coordinates the *dx* dimension + is in geographic delta longitude and hence the physical width of the prism will decrease with latitude + if *dx* stays the same. .. _-F: @@ -164,10 +165,12 @@ Optional Arguments .. _-N: **-N**\ *trackfile* - Specifies individual (x, y[, z]) locations where we wish to compute the predicted value. When this option - is used there are no grids involved and the output data records are written to standard output (see **-bo** for binary output). - If **-Z** is not set then *trackfile* must have 3 columns and we take the *z* value as our observation level; otherwise the level must be set via **-Z**. - **Note**: If **-G** is used to set an output file we will write the output table to that file instead of standard output. + Specifies individual (x, y[, z]) locations where we wish to compute the predicted value. + When this option is used there are no grids involved and the output data records are written + to standard output (see **-bo** for binary output). If **-Z** is not set then *trackfile* must + have 3 columns and we take the *z* value as our observation level; otherwise the level must be + set via **-Z**. **Note**: If **-G** is used to set an output file we will write the output table + to that file instead of standard output. .. _-S: @@ -177,7 +180,8 @@ Optional Arguments .. _-T: **-T**\ *top* - Give name of the top surface grid for a layer we wish to approximate with prisms, or give a constant *z*-level. + Give name of the top surface grid for a layer we wish to approximate with prisms, or give a + constant *z*-level. .. |Add_-V| replace:: |Add_-V_links| .. include:: /explain_-V.rst_ @@ -187,8 +191,8 @@ Optional Arguments .. _-W: **-W**\ *avedens* - Give name of an input grid with spatially varying, vertically-averaged prism densities. Requires **-C** and - the grid must be co-registered with the grid provided by **-S** (or **L** and **-T**). + Give name of an input grid with spatially varying, vertically-averaged prism densities. Requires + **-C** and the grid must be co-registered with the grid provided by **-S** (or **L** and **-T**). .. _-Z: diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 0e5ed285eeb..5b78988046f 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -206,7 +206,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT ns = sscanf (opt->arg, "%lg/%lg", &Ctrl->E.dx, &Ctrl->E.dy); if (ns == 1) Ctrl->E.dy = Ctrl->E.dx; break; - case 'F': /* Seldct geopotential type */ + case 'F': /* Select geopotential type */ n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active); Ctrl->F.active = true; switch (opt->arg[0]) { @@ -258,7 +258,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->H.p1 = Ctrl->H.p + 1; if (c) c[0] = '+'; /* Restore modifiers */ break; - case 'I': /* Grid iincrements */ + case 'I': /* Grid increments */ n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active); Ctrl->I.active = true; n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg); @@ -404,7 +404,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, -2, "Change distance units used, via one or two directives:"); GMT_Usage (API, 3, "h: All x- and y-distances are given in km [meters]."); GMT_Usage (API, 3, "z: All z-distances are given in km [meters]."); - GMT_Usage (API, -2, "Note: All horixontal and/or vertical quantities will be affected by a factor of 1000"); + GMT_Usage (API, -2, "Note: All horizontal and/or vertical quantities will be affected by a factor of 1000"); GMT_Usage (API, 1, "\n-N"); GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally z may be read as well. No grid " "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); @@ -432,7 +432,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) { #define GRAVITATIONAL_CONST_VGG 6.674e-2 /* To convert mGal/m to 0.1 mGal/km requires 1e4 */ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { - /* Geoid anomaly from a single prism [Nagy et al, 2000] */ + /* Geoid anomaly from a single vertical prism [Nagy et al, 2000] */ double n, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double n111, n112, n121, n122, n211, n212, n221, n222; @@ -467,7 +467,7 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou } GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { - /* Gravity anomaly from a single prism [Grant & West, 1965] */ + /* Gravity anomaly from a single vertical prism [Grant & West, 1965] */ double g, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double g111, g112, g121, g122, g211, g212, g221, g222; @@ -502,7 +502,7 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub } GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { - /* Vertical gravity gradient anomaly from a single prism [Kim & Wessel, 2016] */ + /* Vertical gravity gradient anomaly from a single vertical prism [Kim & Wessel, 2016] */ double v, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double v111, v112, v121, v122, v211, v212, v221, v222; @@ -675,7 +675,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid; size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; - if (Ctrl->L.active) { /* Specifed layer base */ + if (Ctrl->L.active) { /* Specified layer base */ if (access (Ctrl->L.file, F_OK)) /* No file, just a constant */ base = atof (Ctrl->L.file); else if ((B = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->L.file, NULL)) == NULL) @@ -944,7 +944,7 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { #pragma omp parallel for private(row,z_level) shared(S,Ctrl,GMT,eval,scl_xy,scl_z,n_prisms,prism,G0) #endif /* Separate the calculation from the output in two separate row-loops since cannot do rec-by-rec output - * with OpenMP due to race condiations that would mess up the output order */ + * with OpenMP due to race conditions that would mess up the output order */ for (row = 0; row < (int64_t)S->n_rows; row++) { /* Calculate attraction at all output locations for this segment */ z_level = (S->n_columns == 3 && !Ctrl->Z.active) ? S->data[GMT_Z][row] : Ctrl->Z.level; /* Default observation z level unless provided in input file */ GMT->hidden.mem_coord[GMT_X][row] = eval (S->data[GMT_X][row] * scl_xy, S->data[GMT_Y][row] * scl_xy, z_level * scl_z, n_prisms, prism, G0); From b678bf1160d3ccb8549bf02ce3fe77e5f4621c9b Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 15 Mar 2022 16:07:55 +0000 Subject: [PATCH 23/27] Add single test for all three outputs Simplify terms by computing the cross products first. --- src/potential/gravprisms.c | 71 ++++++++++++++++++++----------- test/potential/prism_check.sh | 45 ++++++++++++++++++++ test/potential/prism_check_faa.sh | 20 --------- 3 files changed, 92 insertions(+), 44 deletions(-) create mode 100755 test/potential/prism_check.sh delete mode 100755 test/potential/prism_check_faa.sh diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 5b78988046f..ad4df213fd5 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -436,6 +436,9 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou double n, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double n111, n112, n121, n122, n211, n212, n221, n222; + double dx1dy1, dx2dy1, dx1dy2, dx2dy2; + double dx1dz1, dx2dz1, dx1dz2, dx2dz2; + double dy1dz1, dy2dz1, dy1dz2, dy2dz2; /* Square distances */ dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; @@ -450,16 +453,30 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Get cross-terms */ + dx1dy1 = dx1 * dy1; dx2dy1 = dx2 * dy1; dx1dy2 = dx1 * dy2; dx2dy2 = dx2 * dy2; + dx1dz1 = dx1 * dz1; dx2dz1 = dx2 * dz1; dx1dz2 = dx1 * dz2; dx2dz2 = dx2 * dz2; + dy1dz1 = dy1 * dz1; dy2dz1 = dy2 * dz1; dy1dz2 = dy1 * dz2; dy2dz2 = dy2 * dz2; /* Evaluate at dz1 */ - n111 = -(0.5 * (dx1_sq * atan ((dy1 * dz1) / (dz1 * R111)) + dy1_sq * atan ((dx1 * dz1) / (dz1 * R111)) + dz1_sq * atan ((dx1 * dy1) / (dz1 * R111))) - dx1 * dz1 * log (R111 + dy1) - dy1 * dz1 * log (R111 + dx1) - dx1 * dy1 * log (R111 + dz1)); - n112 = +(0.5 * (dx2_sq * atan ((dy1 * dz1) / (dz1 * R112)) + dy1_sq * atan ((dx2 * dz1) / (dz1 * R112)) + dz1_sq * atan ((dx2 * dy1) / (dz1 * R112))) - dx2 * dz1 * log (R112 + dy1) - dy1 * dz1 * log (R112 + dx2) - dx2 * dy1 * log (R112 + dz1)); - n121 = +(0.5 * (dx1_sq * atan ((dy2 * dz1) / (dz1 * R121)) + dy2_sq * atan ((dx1 * dz1) / (dz1 * R121)) + dz1_sq * atan ((dx1 * dy2) / (dz1 * R121))) - dx1 * dz1 * log (R121 + dy2) - dy2 * dz1 * log (R121 + dx1) - dx1 * dy2 * log (R121 + dz1)); - n122 = -(0.5 * (dx2_sq * atan ((dy2 * dz1) / (dz1 * R122)) + dy2_sq * atan ((dx2 * dz1) / (dz1 * R122)) + dz1_sq * atan ((dx2 * dy2) / (dz1 * R122))) - dx2 * dz1 * log (R122 + dy2) - dy2 * dz1 * log (R122 + dx2) - dx2 * dy2 * log (R122 + dz1)); + n111 = -(0.5 * (dx1_sq * atan (dy1dz1 / (dx1 * R111)) + dy1_sq * atan (dx1dz1 / (dy1 * R111)) + dz1_sq * atan (dx1dy1 / (dz1 * R111))) - dx1dz1 * log (R111 + dy1) - dy1dz1 * log (R111 + dx1) - dx1dy1 * log (R111 + dz1)); + n112 = +(0.5 * (dx2_sq * atan (dy1dz1 / (dx2 * R112)) + dy1_sq * atan (dx2dz1 / (dy1 * R112)) + dz1_sq * atan (dx2dy1 / (dz1 * R112))) - dx2dz1 * log (R112 + dy1) - dy1dz1 * log (R112 + dx2) - dx2dy1 * log (R112 + dz1)); + n121 = +(0.5 * (dx1_sq * atan (dy2dz1 / (dx1 * R121)) + dy2_sq * atan (dx1dz1 / (dy2 * R121)) + dz1_sq * atan (dx1dy2 / (dz1 * R121))) - dx1dz1 * log (R121 + dy2) - dy2dz1 * log (R121 + dx1) - dx1dy2 * log (R121 + dz1)); + n122 = -(0.5 * (dx2_sq * atan (dy2dz1 / (dx2 * R122)) + dy2_sq * atan (dx2dz1 / (dy2 * R122)) + dz1_sq * atan (dx2dy2 / (dz1 * R122))) - dx2dz1 * log (R122 + dy2) - dy2dz1 * log (R122 + dx2) - dx2dy2 * log (R122 + dz1)); /* Evaluate at dz2 */ - n211 = +(0.5 * (dx1_sq * atan ((dy1 * dz2) / (dz2 * R211)) + dy1_sq * atan ((dx1 * dz2) / (dz2 * R211)) + dz2_sq * atan ((dx1 * dy1) / (dz2 * R211))) - dx1 * dz2 * log (R211 + dy1) - dy1 * dz2 * log (R211 + dx1) - dx1 * dy1 * log (R211 + dz2)); - n212 = -(0.5 * (dx2_sq * atan ((dy1 * dz2) / (dz2 * R212)) + dy1_sq * atan ((dx2 * dz2) / (dz2 * R212)) + dz2_sq * atan ((dx2 * dy1) / (dz2 * R212))) - dx2 * dz2 * log (R212 + dy1) - dy1 * dz2 * log (R212 + dx2) - dx2 * dy1 * log (R212 + dz2)); - n221 = -(0.5 * (dx1_sq * atan ((dy2 * dz2) / (dz2 * R221)) + dy2_sq * atan ((dx1 * dz2) / (dz2 * R221)) + dz2_sq * atan ((dx1 * dy2) / (dz2 * R221))) - dx1 * dz2 * log (R221 + dy2) - dy2 * dz2 * log (R221 + dx1) - dx1 * dy2 * log (R221 + dz2)); - n222 = +(0.5 * (dx2_sq * atan ((dy2 * dz2) / (dz2 * R222)) + dy2_sq * atan ((dx2 * dz2) / (dz2 * R222)) + dz2_sq * atan ((dx2 * dy2) / (dz2 * R222))) - dx2 * dz2 * log (R222 + dy2) - dy2 * dz2 * log (R222 + dx2) - dx2 * dy2 * log (R222 + dz2)); + n211 = +(0.5 * (dx1_sq * atan (dy1dz2 / (dx1 * R211)) + dy1_sq * atan (dx1dz2 / (dy1 * R211)) + dz2_sq * atan (dx1dy1 / (dz2 * R211))) - dx1dz2 * log (R211 + dy1) - dy1dz2 * log (R211 + dx1) - dx1dy1 * log (R211 + dz2)); + n212 = -(0.5 * (dx2_sq * atan (dy1dz2 / (dx2 * R212)) + dy1_sq * atan (dx2dz2 / (dy1 * R212)) + dz2_sq * atan (dx2dy1 / (dz2 * R212))) - dx2dz2 * log (R212 + dy1) - dy1dz2 * log (R212 + dx2) - dx2dy1 * log (R212 + dz2)); + n221 = -(0.5 * (dx1_sq * atan (dy2dz2 / (dx1 * R221)) + dy2_sq * atan (dx1dz2 / (dy2 * R221)) + dz2_sq * atan (dx1dy2 / (dz2 * R221))) - dx1dz2 * log (R221 + dy2) - dy2dz2 * log (R221 + dx1) - dx1dy2 * log (R221 + dz2)); + n222 = +(0.5 * (dx2_sq * atan (dy2dz2 / (dx2 * R222)) + dy2_sq * atan (dx2dz2 / (dy2 * R222)) + dz2_sq * atan (dx2dy2 / (dz2 * R222))) - dx2dz2 * log (R222 + dy2) - dy2dz2 * log (R222 + dx2) - dx2dy2 * log (R222 + dz2)); + + /* As implemented this function is subject to cancellations and round-off. The n??? terms print the same each time but n differs in the 3rd decimal. */ +//fprintf (stderr, "n111 = %.20lg\n", n111); +//fprintf (stderr, "n112 = %.20lg\n", n112); +//fprintf (stderr, "n121 = %.20lg\n", n121); +//fprintf (stderr, "n122 = %.20lg\n", n122); +//fprintf (stderr, "n211 = %.20lg\n", n211); +//fprintf (stderr, "n212 = %.20lg\n", n212); +//fprintf (stderr, "n221 = %.20lg\n", n221); +//fprintf (stderr, "n222 = %.20lg\n", n222); n = -rho * GRAVITATIONAL_CONST_GEOID * (n111 + n112 + n121 + n122 + n211 + n212 + n221 + n222); @@ -471,6 +488,7 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub double g, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double g111, g112, g121, g122, g211, g212, g221, g222; + double dx1dy1, dx2dy1, dx1dy2, dx2dy2; /* Square distances */ dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; @@ -485,16 +503,18 @@ GMT_LOCAL double gravprism (double dx1, double dx2, double dy1, double dy2, doub R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Get cross-terms */ + dx1dy1 = dx1 * dy1; dx2dy1 = dx2 * dy1; dx1dy2 = dx1 * dy2; dx2dy2 = dx2 * dy2; /* Evaluate at dz1 */ - g111 = -(dz1 * atan ((dx1 * dy1) / (dz1 * R111)) - dx1 * log(R111 + dy1) - dy1 * log (R111 + dx1)); - g112 = +(dz1 * atan ((dx2 * dy1) / (dz1 * R112)) - dx2 * log(R112 + dy1) - dy1 * log (R112 + dx2)); - g121 = +(dz1 * atan ((dx1 * dy2) / (dz1 * R121)) - dx1 * log(R121 + dy2) - dy2 * log (R121 + dx1)); - g122 = -(dz1 * atan ((dx2 * dy2) / (dz1 * R122)) - dx2 * log(R122 + dy2) - dy2 * log (R122 + dx2)); + g111 = -(dz1 * atan (dx1dy1 / (dz1 * R111)) - dx1 * log (R111 + dy1) - dy1 * log (R111 + dx1)); + g112 = +(dz1 * atan (dx2dy1 / (dz1 * R112)) - dx2 * log (R112 + dy1) - dy1 * log (R112 + dx2)); + g121 = +(dz1 * atan (dx1dy2 / (dz1 * R121)) - dx1 * log (R121 + dy2) - dy2 * log (R121 + dx1)); + g122 = -(dz1 * atan (dx2dy2 / (dz1 * R122)) - dx2 * log (R122 + dy2) - dy2 * log (R122 + dx2)); /* Evaluate at dz2 */ - g211 = +(dz2 * atan ((dx1 * dy1) / (dz2 * R211)) - dx1 * log (R211 + dy1) - dy1 * log (R211 + dx1)); - g212 = -(dz2 * atan ((dx2 * dy1) / (dz2 * R212)) - dx2 * log (R212 + dy1) - dy1 * log (R212 + dx2)); - g221 = -(dz2 * atan ((dx1 * dy2) / (dz2 * R221)) - dx1 * log (R221 + dy2) - dy2 * log (R221 + dx1)); - g222 = +(dz2 * atan ((dx2 * dy2) / (dz2 * R222)) - dx2 * log (R222 + dy2) - dy2 * log (R222 + dx2)); + g211 = +(dz2 * atan (dx1dy1 / (dz2 * R211)) - dx1 * log (R211 + dy1) - dy1 * log (R211 + dx1)); + g212 = -(dz2 * atan (dx2dy1 / (dz2 * R212)) - dx2 * log (R212 + dy1) - dy1 * log (R212 + dx2)); + g221 = -(dz2 * atan (dx1dy2 / (dz2 * R221)) - dx1 * log (R221 + dy2) - dy2 * log (R221 + dx1)); + g222 = +(dz2 * atan (dx2dy2 / (dz2 * R222)) - dx2 * log (R222 + dy2) - dy2 * log (R222 + dx2)); g = -rho * GRAVITATIONAL_CONST_FAA * (g111 + g112 + g121 + g122 + g211 + g212 + g221 + g222); @@ -506,6 +526,7 @@ GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, doubl double v, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; double R111, R112, R121, R122, R211, R212, R221, R222; double v111, v112, v121, v122, v211, v212, v221, v222; + double dx1dy1, dx2dy1, dx1dy2, dx2dy2; /* Square distances */ dx1_sq = dx1 * dx1; dx2_sq = dx2 * dx2; @@ -520,16 +541,18 @@ GMT_LOCAL double vggprism (double dx1, double dx2, double dy1, double dy2, doubl R212 = sqrt (dx2_sq + dy1_sq + dz2_sq); R221 = sqrt (dx1_sq + dy2_sq + dz2_sq); R222 = sqrt (dx2_sq + dy2_sq + dz2_sq); + /* Get cross-terms */ + dx1dy1 = dx1 * dy1; dx2dy1 = dx2 * dy1; dx1dy2 = dx1 * dy2; dx2dy2 = dx2 * dy2; /* Evaluate at dz1 */ - v111 = -atan ((dx1 * dy1) / (dz1 * R111)); - v112 = +atan ((dx2 * dy1) / (dz1 * R112)); - v121 = +atan ((dx1 * dy2) / (dz1 * R121)); - v122 = -atan ((dx2 * dy2) / (dz1 * R122)); + v111 = -atan (dx1dy1 / (dz1 * R111)); + v112 = +atan (dx2dy1 / (dz1 * R112)); + v121 = +atan (dx1dy2 / (dz1 * R121)); + v122 = -atan (dx2dy2 / (dz1 * R122)); /* Evaluate at dz2 */ - v211 = +atan ((dx1 * dy1) / (dz2 * R211)); - v212 = -atan ((dx2 * dy1) / (dz2 * R212)); - v221 = -atan ((dx1 * dy2) / (dz2 * R221)); - v222 = +atan ((dx2 * dy2) / (dz2 * R222)); + v211 = +atan (dx1dy1 / (dz2 * R211)); + v212 = -atan (dx2dy1 / (dz2 * R212)); + v221 = -atan (dx1dy2 / (dz2 * R221)); + v222 = +atan (dx2dy2 / (dz2 * R222)); v = -rho * GRAVITATIONAL_CONST_VGG * (v111 + v112 + v121 + v122 + v211 + v212 + v221 + v222); diff --git a/test/potential/prism_check.sh b/test/potential/prism_check.sh new file mode 100755 index 00000000000..1add3ae8c45 --- /dev/null +++ b/test/potential/prism_check.sh @@ -0,0 +1,45 @@ +#!/usr/bin/env bash +# Check that the calculation for faa, vgg, and going over a prism is correct +# Results validated by hand with separate Matlab code as well as talwani3d +# We test both Cartesian (in km) and geographic (degrees) + +cat <<- EOF > answer.txt +10.000000 10.000000 7000.000000 5.186362 +10.000000 10.000000 7000.000000 -7.576966 +10.00 10.00 7000.00 0.23 +0.100000 0.100000 7000.000000 4.861885 +0.100000 0.100000 7000.000000 -7.670554 +0.10 0.10 7000.00 0.26 +EOF +# Prism centered at (0,0) ranging from z = 0-5000 m +echo "0 0 0 5000" > p.txt +# Cartesian observation point at (10, 10) km at 7000 m elevation +echo " 10 10 7000" > c.txt +# Geographic observation point at lon/lat (0.1, 0.1) at 7000 m elevation +echo "0.1 0.1 7000" > g.txt +# Cartesian prism has dimensions 10x10 km +gmt gravprisms p.txt -E10 -Nc.txt -Mh -D1000 -Ff --FORMAT_FLOAT_OUT=%.6lf > result.txt +gmt gravprisms p.txt -E10 -Nc.txt -Mh -D1000 -Fv --FORMAT_FLOAT_OUT=%.6lf >> result.txt +gmt gravprisms p.txt -E10 -Nc.txt -Mh -D1000 -Fn --FORMAT_FLOAT_OUT=%.2lf >> result.txt +# Geographic prism has dimension 0.1 degrees by 0.1 degrees in longitude/latitude +gmt gravprisms p.txt -E0.1 -Ng.txt -fg -D1000 -Ff --FORMAT_FLOAT_OUT=%.6lf >> result.txt +gmt gravprisms p.txt -E0.1 -Ng.txt -fg -D1000 -Fv --FORMAT_FLOAT_OUT=%.6lf >> result.txt +gmt gravprisms p.txt -E0.1 -Ng.txt -fg -D1000 -Fn --FORMAT_FLOAT_OUT=%.2lf >> result.txt +diff result.txt answer.txt --strip-trailing-cr > fail +## talwani3d solutions: +#cat <<- EOF > square.txt +#-5 -5 +#5 -5 +#5 5 +#-5 5 +#-5 -5 +#EOF +#gmt math -T0/5000/10 -o0 T = z.txt +#rm -f body.txt +#while read z; do +# echo "> $z 1000" >> body.txt +# cat square.txt >> body.txt +#done < z.txt +#gmt talwani3d body.txt -Fg -Mh -Z7 -Nc.txt +#gmt talwani3d body.txt -Fv -Mh -Z7 -Nc.txt +#gmt talwani3d body.txt -Fn -Mh -Z7 -Nc.txt diff --git a/test/potential/prism_check_faa.sh b/test/potential/prism_check_faa.sh deleted file mode 100755 index 66c7f2c4058..00000000000 --- a/test/potential/prism_check_faa.sh +++ /dev/null @@ -1,20 +0,0 @@ -#!/usr/bin/env bash -# Check that the calculation for gravity over a prism is correct -# Results validated by hand with separate Matlab code -# We test both Cartesian (in km) and geographic (degrees) - -cat <<- EOF > answer.txt -10.000000 10.000000 7000.000000 5.186362 -0.100000 0.100000 7000.000000 4.861885 -EOF -# Prism centered at (0,0) ranging from z = 0-5000 m -echo "0 0 0 5000" > p.txt -# Cartesian observation point at (10, 10) km at 7000 m elevation -echo " 10 10 7000" > c.txt -# Geopgraphic observation point at lon/lat (0.1, 0.1) at 7000 m elevation -echo "0.1 0.1 7000" > g.txt -# Cartesian prism has dimensions 10x10 km -gmt gravprisms p.txt -E10 -Nc.txt -Mh -D1000 --FORMAT_FLOAT_OUT=%.6lf > result.txt -# Geographic prism has dimension 0.1 degrees by 0.1 degrees in longitude/latitude -gmt gravprisms p.txt -E0.1 -Ng.txt -fg -D1000 --FORMAT_FLOAT_OUT=%.6lf >> result.txt -diff result.txt answer.txt --strip-trailing-cr > fail From ccb7dea0b6a523c806f39c95e640c5bcb2374a45 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 15 Mar 2022 21:39:01 +0000 Subject: [PATCH 24/27] Update gravprisms.c --- src/potential/gravprisms.c | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index ad4df213fd5..4344278f3c0 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -635,13 +635,16 @@ GMT_LOCAL double gravprisms_get_one_v_output_geo (double x, double y, double z, return (v); /* Converted units from mGal/m to Eotvos = 0.1 mGal/km */ } -GMT_LOCAL double gravprisms_rho (struct GRAVPRISMS_CTRL *Ctrl, double z, double h) { - /* Passing in the current seamount's normalized height, h(r) and the current normalized radial point z(r). - * We return the density at this point inside the seamount with reference height 1 */ - double u = (h - z) / Ctrl->H.H; /* Normalized depth into seamount. */ - double q = (Ctrl->H.H - h) / Ctrl->H.H; /* Water depth to flank point */ - - double rho = Ctrl->H.rho_l + Ctrl->H.densify * q + Ctrl->H.del_rho * pow (u, Ctrl->H.p); +GMT_LOCAL double gravprisms_mean_density (struct GRAVPRISMS_CTRL *Ctrl, double h, double z1, double z2) { + /* Passing in the current seamounts height h(r) and the two depths z1(r) and z2(r) defining a prism. + * When doing the whole seamount we pass z2 = 0 and z1 = h(r). + * The vertically averaged density for this radial position and range of z is returned */ + double u1 = (h - z1) / Ctrl->H.H; + double u2 = (h - z2) / Ctrl->H.H; + double q = (Ctrl->H.H - h) / Ctrl->H.H; + double dz = z2 - z1; /* Prism height */ + + double rho = Ctrl->H.rho_l + Ctrl->H.densify * q + Ctrl->H.del_rho * Ctrl->H.H * (pow (u1, Ctrl->H.p1) - pow (u2, Ctrl->H.p1)) / (dz * (Ctrl->H.p1)); return (rho); } @@ -738,7 +741,8 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if (z_next <= z_prev) z_next += Ctrl->C.dz; /* Can happen if z1 is a multiple of dz */ else if (z_next > z2) z_next = z2; /* At the top, clip to limit */ z_mid = 0.5 * (z_prev + z_next); /* Middle of prism - used to look up density */ - rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); + //rho = gravprisms_rho (Ctrl, z_mid, H->data[node]); + rho = gravprisms_mean_density (Ctrl, H->data[node], z_prev, z_next); } else { /* Constant density rho (set above via -D) or by Rho (via -W), just need a single prism per location */ z_next = z2; From eb35afc8a6904b2b1568238b5111d3549c57dd59 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Wed, 16 Mar 2022 10:03:04 +0000 Subject: [PATCH 25/27] Update gravprisms.c --- src/potential/gravprisms.c | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 4344278f3c0..fc5ece76a9a 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -431,6 +431,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) { #define GRAVITATIONAL_CONST_GEOID 6.674e-11 /* This gives geoid in meter once we divide by g0 */ #define GRAVITATIONAL_CONST_VGG 6.674e-2 /* To convert mGal/m to 0.1 mGal/km requires 1e4 */ +#define zatan(a,b) ((fabs(b) < GMT_CONV15_LIMIT) ? 0.0 : atan(a/b)) + GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { /* Geoid anomaly from a single vertical prism [Nagy et al, 2000] */ double n, dx1_sq, dx2_sq, dy1_sq, dy2_sq, dz1_sq, dz2_sq; @@ -458,15 +460,15 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou dx1dz1 = dx1 * dz1; dx2dz1 = dx2 * dz1; dx1dz2 = dx1 * dz2; dx2dz2 = dx2 * dz2; dy1dz1 = dy1 * dz1; dy2dz1 = dy2 * dz1; dy1dz2 = dy1 * dz2; dy2dz2 = dy2 * dz2; /* Evaluate at dz1 */ - n111 = -(0.5 * (dx1_sq * atan (dy1dz1 / (dx1 * R111)) + dy1_sq * atan (dx1dz1 / (dy1 * R111)) + dz1_sq * atan (dx1dy1 / (dz1 * R111))) - dx1dz1 * log (R111 + dy1) - dy1dz1 * log (R111 + dx1) - dx1dy1 * log (R111 + dz1)); - n112 = +(0.5 * (dx2_sq * atan (dy1dz1 / (dx2 * R112)) + dy1_sq * atan (dx2dz1 / (dy1 * R112)) + dz1_sq * atan (dx2dy1 / (dz1 * R112))) - dx2dz1 * log (R112 + dy1) - dy1dz1 * log (R112 + dx2) - dx2dy1 * log (R112 + dz1)); - n121 = +(0.5 * (dx1_sq * atan (dy2dz1 / (dx1 * R121)) + dy2_sq * atan (dx1dz1 / (dy2 * R121)) + dz1_sq * atan (dx1dy2 / (dz1 * R121))) - dx1dz1 * log (R121 + dy2) - dy2dz1 * log (R121 + dx1) - dx1dy2 * log (R121 + dz1)); - n122 = -(0.5 * (dx2_sq * atan (dy2dz1 / (dx2 * R122)) + dy2_sq * atan (dx2dz1 / (dy2 * R122)) + dz1_sq * atan (dx2dy2 / (dz1 * R122))) - dx2dz1 * log (R122 + dy2) - dy2dz1 * log (R122 + dx2) - dx2dy2 * log (R122 + dz1)); + n111 = -(0.5 * (dx1_sq * zatan (dy1dz1, (dx1 * R111)) + dy1_sq * zatan (dx1dz1, (dy1 * R111)) + dz1_sq * zatan (dx1dy1, (dz1 * R111))) - dx1dz1 * log (R111 + dy1) - dy1dz1 * log (R111 + dx1) - dx1dy1 * log (R111 + dz1)); + n112 = +(0.5 * (dx2_sq * zatan (dy1dz1, (dx2 * R112)) + dy1_sq * zatan (dx2dz1, (dy1 * R112)) + dz1_sq * zatan (dx2dy1, (dz1 * R112))) - dx2dz1 * log (R112 + dy1) - dy1dz1 * log (R112 + dx2) - dx2dy1 * log (R112 + dz1)); + n121 = +(0.5 * (dx1_sq * zatan (dy2dz1, (dx1 * R121)) + dy2_sq * zatan (dx1dz1, (dy2 * R121)) + dz1_sq * zatan (dx1dy2, (dz1 * R121))) - dx1dz1 * log (R121 + dy2) - dy2dz1 * log (R121 + dx1) - dx1dy2 * log (R121 + dz1)); + n122 = -(0.5 * (dx2_sq * zatan (dy2dz1, (dx2 * R122)) + dy2_sq * zatan (dx2dz1, (dy2 * R122)) + dz1_sq * zatan (dx2dy2, (dz1 * R122))) - dx2dz1 * log (R122 + dy2) - dy2dz1 * log (R122 + dx2) - dx2dy2 * log (R122 + dz1)); /* Evaluate at dz2 */ - n211 = +(0.5 * (dx1_sq * atan (dy1dz2 / (dx1 * R211)) + dy1_sq * atan (dx1dz2 / (dy1 * R211)) + dz2_sq * atan (dx1dy1 / (dz2 * R211))) - dx1dz2 * log (R211 + dy1) - dy1dz2 * log (R211 + dx1) - dx1dy1 * log (R211 + dz2)); - n212 = -(0.5 * (dx2_sq * atan (dy1dz2 / (dx2 * R212)) + dy1_sq * atan (dx2dz2 / (dy1 * R212)) + dz2_sq * atan (dx2dy1 / (dz2 * R212))) - dx2dz2 * log (R212 + dy1) - dy1dz2 * log (R212 + dx2) - dx2dy1 * log (R212 + dz2)); - n221 = -(0.5 * (dx1_sq * atan (dy2dz2 / (dx1 * R221)) + dy2_sq * atan (dx1dz2 / (dy2 * R221)) + dz2_sq * atan (dx1dy2 / (dz2 * R221))) - dx1dz2 * log (R221 + dy2) - dy2dz2 * log (R221 + dx1) - dx1dy2 * log (R221 + dz2)); - n222 = +(0.5 * (dx2_sq * atan (dy2dz2 / (dx2 * R222)) + dy2_sq * atan (dx2dz2 / (dy2 * R222)) + dz2_sq * atan (dx2dy2 / (dz2 * R222))) - dx2dz2 * log (R222 + dy2) - dy2dz2 * log (R222 + dx2) - dx2dy2 * log (R222 + dz2)); + n211 = +(0.5 * (dx1_sq * zatan (dy1dz2, (dx1 * R211)) + dy1_sq * zatan (dx1dz2, (dy1 * R211)) + dz2_sq * zatan (dx1dy1, (dz2 * R211))) - dx1dz2 * log (R211 + dy1) - dy1dz2 * log (R211 + dx1) - dx1dy1 * log (R211 + dz2)); + n212 = -(0.5 * (dx2_sq * zatan (dy1dz2, (dx2 * R212)) + dy1_sq * zatan (dx2dz2, (dy1 * R212)) + dz2_sq * zatan (dx2dy1, (dz2 * R212))) - dx2dz2 * log (R212 + dy1) - dy1dz2 * log (R212 + dx2) - dx2dy1 * log (R212 + dz2)); + n221 = -(0.5 * (dx1_sq * zatan (dy2dz2, (dx1 * R221)) + dy2_sq * zatan (dx1dz2, (dy2 * R221)) + dz2_sq * zatan (dx1dy2, (dz2 * R221))) - dx1dz2 * log (R221 + dy2) - dy2dz2 * log (R221 + dx1) - dx1dy2 * log (R221 + dz2)); + n222 = +(0.5 * (dx2_sq * zatan (dy2dz2, (dx2 * R222)) + dy2_sq * zatan (dx2dz2, (dy2 * R222)) + dz2_sq * zatan (dx2dy2, (dz2 * R222))) - dx2dz2 * log (R222 + dy2) - dy2dz2 * log (R222 + dx2) - dx2dy2 * log (R222 + dz2)); /* As implemented this function is subject to cancellations and round-off. The n??? terms print the same each time but n differs in the 3rd decimal. */ //fprintf (stderr, "n111 = %.20lg\n", n111); From 0fb95ed0392ffa6764e167deeac14b5979854270 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Wed, 16 Mar 2022 11:07:54 +0000 Subject: [PATCH 26/27] Update gravprisms.c --- src/potential/gravprisms.c | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index fc5ece76a9a..b4cfe8ef8ad 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -406,8 +406,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "z: All z-distances are given in km [meters]."); GMT_Usage (API, -2, "Note: All horizontal and/or vertical quantities will be affected by a factor of 1000"); GMT_Usage (API, 1, "\n-N"); - GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally z may be read as well. No grid " - "is produced and output (x,y,z,g) is written to standard output (see -bo for binary output)."); + GMT_Usage (API, -2, "File with output locations (x,y) where a calculation is requested. Optionally, z may be read as well if -Z not set. No grid " + "is produced and output (x,y,z,g) is written to standard output (but see -G; also see -bo for binary output)."); GMT_Option (API, "R"); GMT_Usage (API, 1, "\n-S"); GMT_Usage (API, -2, "Set the full surface grid height for making prisms of the entire feature (or if -H is set) (or use -L and -T to select a layer); requires -C"); @@ -427,11 +427,13 @@ static int usage (struct GMTAPI_CTRL *API, int level) { return (GMT_MODULE_USAGE); } +#define GRAVITATIONAL_CONST_GEOID 6.674e-11 /* To get geoid in meter we divide by g0 in gravprisms_get_one_n_output */ #define GRAVITATIONAL_CONST_FAA 6.674e-6 /* To convert m/s^2 to mGal requires 1e5 */ -#define GRAVITATIONAL_CONST_GEOID 6.674e-11 /* This gives geoid in meter once we divide by g0 */ -#define GRAVITATIONAL_CONST_VGG 6.674e-2 /* To convert mGal/m to 0.1 mGal/km requires 1e4 */ +#define GRAVITATIONAL_CONST_VGG 6.674e-2 /* To convert mGal/m to 0.1 mGal/km requires an additional 1e4 */ -#define zatan(a,b) ((fabs(b) < GMT_CONV15_LIMIT) ? 0.0 : atan(a/b)) +/* Geoid: Carefully checking terms to avoid divisions by zero in atan or log (zero) */ +#define zatan(a,b) ((fabs(b) < GMT_CONV15_LIMIT) ? 0.0 : atan(a/b)) /* For safe atan (a/b) */ +#define zlog(a,b) ((fabs(b) < GMT_CONV15_LIMIT) ? 0.0 : a * log(b)) /* For safe a * log(b) */ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, double dz1, double dz2, double rho) { /* Geoid anomaly from a single vertical prism [Nagy et al, 2000] */ @@ -460,25 +462,15 @@ GMT_LOCAL double geoidprism (double dx1, double dx2, double dy1, double dy2, dou dx1dz1 = dx1 * dz1; dx2dz1 = dx2 * dz1; dx1dz2 = dx1 * dz2; dx2dz2 = dx2 * dz2; dy1dz1 = dy1 * dz1; dy2dz1 = dy2 * dz1; dy1dz2 = dy1 * dz2; dy2dz2 = dy2 * dz2; /* Evaluate at dz1 */ - n111 = -(0.5 * (dx1_sq * zatan (dy1dz1, (dx1 * R111)) + dy1_sq * zatan (dx1dz1, (dy1 * R111)) + dz1_sq * zatan (dx1dy1, (dz1 * R111))) - dx1dz1 * log (R111 + dy1) - dy1dz1 * log (R111 + dx1) - dx1dy1 * log (R111 + dz1)); - n112 = +(0.5 * (dx2_sq * zatan (dy1dz1, (dx2 * R112)) + dy1_sq * zatan (dx2dz1, (dy1 * R112)) + dz1_sq * zatan (dx2dy1, (dz1 * R112))) - dx2dz1 * log (R112 + dy1) - dy1dz1 * log (R112 + dx2) - dx2dy1 * log (R112 + dz1)); - n121 = +(0.5 * (dx1_sq * zatan (dy2dz1, (dx1 * R121)) + dy2_sq * zatan (dx1dz1, (dy2 * R121)) + dz1_sq * zatan (dx1dy2, (dz1 * R121))) - dx1dz1 * log (R121 + dy2) - dy2dz1 * log (R121 + dx1) - dx1dy2 * log (R121 + dz1)); - n122 = -(0.5 * (dx2_sq * zatan (dy2dz1, (dx2 * R122)) + dy2_sq * zatan (dx2dz1, (dy2 * R122)) + dz1_sq * zatan (dx2dy2, (dz1 * R122))) - dx2dz1 * log (R122 + dy2) - dy2dz1 * log (R122 + dx2) - dx2dy2 * log (R122 + dz1)); + n111 = -(0.5 * (dx1_sq * zatan (dy1dz1, (dx1 * R111)) + dy1_sq * zatan (dx1dz1, (dy1 * R111)) + dz1_sq * zatan (dx1dy1, (dz1 * R111))) - zlog (dx1dz1, R111 + dy1) - zlog (dy1dz1, R111 + dx1) - zlog (dx1dy1, R111 + dz1)); + n112 = +(0.5 * (dx2_sq * zatan (dy1dz1, (dx2 * R112)) + dy1_sq * zatan (dx2dz1, (dy1 * R112)) + dz1_sq * zatan (dx2dy1, (dz1 * R112))) - zlog (dx2dz1, R112 + dy1) - zlog (dy1dz1, R112 + dx2) - zlog (dx2dy1, R112 + dz1)); + n121 = +(0.5 * (dx1_sq * zatan (dy2dz1, (dx1 * R121)) + dy2_sq * zatan (dx1dz1, (dy2 * R121)) + dz1_sq * zatan (dx1dy2, (dz1 * R121))) - zlog (dx1dz1, R121 + dy2) - zlog (dy2dz1, R121 + dx1) - zlog (dx1dy2, R121 + dz1)); + n122 = -(0.5 * (dx2_sq * zatan (dy2dz1, (dx2 * R122)) + dy2_sq * zatan (dx2dz1, (dy2 * R122)) + dz1_sq * zatan (dx2dy2, (dz1 * R122))) - zlog (dx2dz1, R122 + dy2) - zlog (dy2dz1, R122 + dx2) - zlog (dx2dy2, R122 + dz1)); /* Evaluate at dz2 */ - n211 = +(0.5 * (dx1_sq * zatan (dy1dz2, (dx1 * R211)) + dy1_sq * zatan (dx1dz2, (dy1 * R211)) + dz2_sq * zatan (dx1dy1, (dz2 * R211))) - dx1dz2 * log (R211 + dy1) - dy1dz2 * log (R211 + dx1) - dx1dy1 * log (R211 + dz2)); - n212 = -(0.5 * (dx2_sq * zatan (dy1dz2, (dx2 * R212)) + dy1_sq * zatan (dx2dz2, (dy1 * R212)) + dz2_sq * zatan (dx2dy1, (dz2 * R212))) - dx2dz2 * log (R212 + dy1) - dy1dz2 * log (R212 + dx2) - dx2dy1 * log (R212 + dz2)); - n221 = -(0.5 * (dx1_sq * zatan (dy2dz2, (dx1 * R221)) + dy2_sq * zatan (dx1dz2, (dy2 * R221)) + dz2_sq * zatan (dx1dy2, (dz2 * R221))) - dx1dz2 * log (R221 + dy2) - dy2dz2 * log (R221 + dx1) - dx1dy2 * log (R221 + dz2)); - n222 = +(0.5 * (dx2_sq * zatan (dy2dz2, (dx2 * R222)) + dy2_sq * zatan (dx2dz2, (dy2 * R222)) + dz2_sq * zatan (dx2dy2, (dz2 * R222))) - dx2dz2 * log (R222 + dy2) - dy2dz2 * log (R222 + dx2) - dx2dy2 * log (R222 + dz2)); - - /* As implemented this function is subject to cancellations and round-off. The n??? terms print the same each time but n differs in the 3rd decimal. */ -//fprintf (stderr, "n111 = %.20lg\n", n111); -//fprintf (stderr, "n112 = %.20lg\n", n112); -//fprintf (stderr, "n121 = %.20lg\n", n121); -//fprintf (stderr, "n122 = %.20lg\n", n122); -//fprintf (stderr, "n211 = %.20lg\n", n211); -//fprintf (stderr, "n212 = %.20lg\n", n212); -//fprintf (stderr, "n221 = %.20lg\n", n221); -//fprintf (stderr, "n222 = %.20lg\n", n222); + n211 = +(0.5 * (dx1_sq * zatan (dy1dz2, (dx1 * R211)) + dy1_sq * zatan (dx1dz2, (dy1 * R211)) + dz2_sq * zatan (dx1dy1, (dz2 * R211))) - zlog (dx1dz2, R211 + dy1) - zlog (dy1dz2, R211 + dx1) - zlog (dx1dy1, R211 + dz2)); + n212 = -(0.5 * (dx2_sq * zatan (dy1dz2, (dx2 * R212)) + dy1_sq * zatan (dx2dz2, (dy1 * R212)) + dz2_sq * zatan (dx2dy1, (dz2 * R212))) - zlog (dx2dz2, R212 + dy1) - zlog (dy1dz2, R212 + dx2) - zlog (dx2dy1, R212 + dz2)); + n221 = -(0.5 * (dx1_sq * zatan (dy2dz2, (dx1 * R221)) + dy2_sq * zatan (dx1dz2, (dy2 * R221)) + dz2_sq * zatan (dx1dy2, (dz2 * R221))) - zlog (dx1dz2, R221 + dy2) - zlog (dy2dz2, R221 + dx1) - zlog (dx1dy2, R221 + dz2)); + n222 = +(0.5 * (dx2_sq * zatan (dy2dz2, (dx2 * R222)) + dy2_sq * zatan (dx2dz2, (dy2 * R222)) + dz2_sq * zatan (dx2dy2, (dz2 * R222))) - zlog (dx2dz2, R222 + dy2) - zlog (dy2dz2, R222 + dx2) - zlog (dx2dy2, R222 + dz2)); n = -rho * GRAVITATIONAL_CONST_GEOID * (n111 + n112 + n121 + n122 + n211 + n212 + n221 + n222); From 9ed1062d3a7b032902c658d99048937d98983734 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Wed, 16 Mar 2022 12:44:28 +0000 Subject: [PATCH 27/27] Update gravprisms.rst --- doc/rst/source/supplements/potential/gravprisms.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index a956e59b385..ced53ac7373 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -259,6 +259,15 @@ for the same model (at 30S) is written to n_crossing.txt by gmt gravprisms -Ncrossing.txt -Mh prisms.txt -D1700 -Fn-30 -Gn_crossing.txt +Note +---- + +The analytical expression for the geoid over a vertical prism (Nagy et al., 2000) is +fairly involved and contains 48 terms. Due to various cancellations the end result +is more unstable than similar expressions for gravity and VGG. Be aware that the +result may have less significant digits that you may expect. + + References ----------