diff --git a/doc/rst/source/supplements/potential/gmtgravmag3d.rst b/doc/rst/source/supplements/potential/gmtgravmag3d.rst index 3c6d50b8d12..7662b03ab94 100644 --- a/doc/rst/source/supplements/potential/gmtgravmag3d.rst +++ b/doc/rst/source/supplements/potential/gmtgravmag3d.rst @@ -12,13 +12,14 @@ Synopsis .. include:: ../../common_SYN_OPTs.rst_ -**gmt gravmag3d** *xyz_file* |-T|\ **+v**\ *vert_file* OR |-T|\ **+r\|+s**\ *raw_file* OR |-M|\ **+s**\ *body,params* +**gmt gravmag3d** *xyz_file* |-T|\ *vert_file*\ **+v**\ [**+n**]\ OR |-T|\ *raw_file*\ **+r\|+s**\ [**+n**]\ OR |-M|\ **+s**\ *body,params* [ |-C|\ *density* ] [ |-E|\ *thickness* ] [ |-F|\ *xy_file* ] [ |-G|\ *outgrid* ] [ |-H|\ *f_dec*/*f_dip*/*m_int*/*m_dec*/*m_dip* ] [ |-L|\ *z_observation* ] +[ |-Q| ] [ |-S|\ *radius* ] [ |-Z|\ *level* ] [ |SYN_OPT-V| ] @@ -97,8 +98,8 @@ Required Arguments (not all) .. _-T: -**-T+v**\ *vert_file* (must have when passing a *xyz_file*) OR **-T+r\|+s**\ *raw_file* - Gives names of a xyz and vertex (**-T+v**\ *vert_file*) files defining a closed surface. +**-T**\ *vert_file*\ **+v**\ [**+n**] (must have when passing a *xyz_file*) OR **-T**\ *raw_file*\ **+r\|+s**\ [**+n**]\ + Gives names of a xyz and vertex (**-T**\ *vert_file*\ **+v**) files defining a closed surface. The file formats correspond to the output of the :doc:`triangulate ` program. The *xyz* file can have 3, 4, 5, 6 or 8 columns. In first case (3 columns) the magnetization (or density) are assumed constant (controlled by |-C| or |-H|). Following cases are: 4 columns -> 4rth col magnetization intensity; @@ -108,6 +109,11 @@ Required Arguments (not all) vertex of each triangle. Alternatively, the **-T+s** option indicates that the surface file is in the ASCII STL (Stereo Lithographic) format. These two type of files are used to provide a closed surface. + The closed surface formats (_STL_ or _raw_) are assumed to be provided with the facets (triangles) following the + counter-clockwise order. If that is not the case, _i.e._ they are clockwise oriented, append the **+n** modifier to + bring them to the expected order. However, this order may not be easy to check. In case of doubt, compute the gravity + anomaly caused by the body and see if it has the expected signal. + Optional Arguments ------------------ @@ -119,20 +125,27 @@ Optional Arguments .. _-E: **-E**\ [*thickness*] - give layer thickness in m [Default = 0 m]. Use this option only when + Give layer thickness in m [Default = 0 m]. Use this option only when the triangles describe a non-closed surface and you want the anomaly of a constant thickness layer. .. _-L: **-L**\ [*z_observation*] - sets level of observation [Default = 0]. That is the height (z) at + Sets level of observation [Default = 0]. That is the height (z) at which anomalies are computed. +.. _-Q: + +**-Q** + Use this option if the indices file (**-T+v**) is 1-based instead of the default (C) 0-based. + Normally this option is useful specially when using GMT.jl or GMTMEX where indices are 1-based + and bodies may have been created in those environments. + .. _-S: **-S**\ *radius* - search radius in km. Triangle centroids that are further away than + Search radius in km. Triangle centroids that are further away than *radius* from current output point will not be taken into account. Use this option to speed up computation at expenses of a less accurate result. @@ -140,7 +153,7 @@ Optional Arguments .. _-Z: **-Z**\ [*level*] - level of reference plane [Default = 0]. Use this option when the + Level of reference plane [Default = 0]. Use this option when the triangles describe a non-closed surface and the volume is defined from each triangle and this reference level. An example will be the hater depth to compute a Bouguer anomaly. diff --git a/src/gmt_api.c b/src/gmt_api.c index 99e5e93505b..64ccd7f4e98 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -3780,7 +3780,7 @@ GMT_LOCAL bool gmtapi_vector_data_must_be_duplicated (struct GMTAPI_CTRL *API, s } /*! . */ -GMT_LOCAL struct GMT_DATASET * gmtapi_import_dataset (struct GMTAPI_CTRL *API, int object_ID, unsigned int mode) { +GMT_LOCAL struct GMT_DATASET *gmtapi_import_dataset(struct GMTAPI_CTRL *API, int object_ID, unsigned int mode) { /* Does the actual work of loading in the entire virtual data set (possibly via many sources) * If object_ID == GMT_NOTSET we get all registered input tables, otherwise we just get the one requested. * Note: Memory is allocated for the Dataset except for method GMT_IS_REFERENCE. @@ -5882,7 +5882,7 @@ GMT_LOCAL int gmtapi_export_grid (struct GMTAPI_CTRL *API, int object_ID, unsign case GMT_IS_REFERENCE: /* GMT grid and header in a GMT_GRID container object - just pass the reference */ if (S_obj->region) return (gmtlib_report_error (API, GMT_SUBSET_NOT_ALLOWED)); if (mode & GMT_CONTAINER_ONLY) return (gmtlib_report_error (API, GMT_NOT_A_VALID_MODE)); - GMT_Report (API, GMT_MSG_INFORMATION, "Referencing grid data to GMT_GRID memory location\n"); + GMT_Report(API, GMT_MSG_DEBUG, "Referencing grid data to GMT_GRID memory location\n"); gmt_grd_zminmax (GMT, G_obj->header, G_obj->data); /* Must set zmin/zmax since we are not writing to file */ gmt_BC_init (GMT, G_obj->header); /* Initialize grid interpolation and boundary condition parameters */ if (gmt_M_err_pass (GMT, gmt_grd_BC_set (GMT, G_obj, GMT_OUT), "Grid memory")) return (gmtlib_report_error (API, GMT_GRID_BC_ERROR)); /* Set boundary conditions */ diff --git a/src/gmt_grdio.c b/src/gmt_grdio.c index 129db7cff7a..9d81affaf4b 100644 --- a/src/gmt_grdio.c +++ b/src/gmt_grdio.c @@ -1392,7 +1392,7 @@ int gmtlib_get_grdtype (struct GMT_CTRL *GMT, unsigned int direction, struct GMT return (GMT_GRID_CARTESIAN); } } - GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Cartesian %s grid\n", dir[direction]); + GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Cartesian %s grid\n", dir[direction]); return (GMT_GRID_CARTESIAN); } diff --git a/src/potential/gmtgravmag3d.c b/src/potential/gmtgravmag3d.c index a3ff4315fb8..807a2508743 100644 --- a/src/potential/gmtgravmag3d.c +++ b/src/potential/gmtgravmag3d.c @@ -88,6 +88,7 @@ struct GMTGRAVMAG3D_CTRL { bool triangulate; bool raw; bool stl; + bool body; bool m_var, m_var1, m_var2, m_var3, m_var4; char *xyz_file; char *t_file; @@ -98,6 +99,8 @@ struct GMTGRAVMAG3D_CTRL { bool is_geog; double d_to_m, *mag_int, lon_0, lat_0; } box; + bool swap; + unsigned int zero_base; int n_triang, n_vert, n_raw_triang; int npts_circ; /* Number of points in which a circle is descretized. */ int n_slices; /* Spheres and Ellipsoids are made by 2*n_slices. Bells are made by n_slices. */ @@ -174,12 +177,13 @@ GMT_LOCAL void set_center (struct GMTGRAVMAG3D_CTRL *Ctrl); GMT_LOCAL int facet_triangulate (struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body_verts, unsigned int i, bool bat); GMT_LOCAL int facet_raw (struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body_verts, unsigned int i, bool geo); GMT_LOCAL int check_triang_cw (struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, unsigned int type); +GMT_LOCAL void swap_triang(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n_tri); 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 [] -T+v | -T+r|+s | -M+s/ [-C] [-E] " - "[-F] [-G%s] [-H///] [%s] [-L] [%s] " + GMT_Usage (API, 0, "usage: %s [
] -T+v|+f[+n] | -T+r|+s[+n] | -M+s, [-C] [-E] " + "[-F] [-G%s] [-H///] [%s] [-L] [%s] [-Q] " "[-S] [-Z] [%s] [-fg] [%s] [%s] [%s] [%s]%s[%s]\n", name, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_r_OPT, GMT_xg_OPT, GMT_PAR_OPT); @@ -188,15 +192,17 @@ 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 data files (in ASCII, binary, netCDF) with data; see -T for format. If no files are given, standard input is read"); - GMT_Usage (API, 1, "\n-T+v | -T+r|+s"); - GMT_Usage (API, -2, "Give names of xyz and vertex (-T+v) files defining a closed surface. " + GMT_Usage (API, 1, "\n-T+v|+f[+n] | -T+r|+s[+n]"); + GMT_Usage (API, -2, "Give names of xyz and vertex (-T+v|+f) files defining a surface that may or not be closed. " "If has more then 3 columns it means variable magnetization; see HTML docs for more details. " "The file formats correspond to the output of the triangulate program. " "Alternatively, use directives to indicate specific formats:"); + GMT_Usage (API, 3, "+f: This form assumes that the vertices define a closed body (the +v may or not be closed)."); GMT_Usage (API, 3, "+r: Append in raw triangle format (x1 y1 z1 x2 ... z3)."); GMT_Usage (API, 3, "+s: Append in STL format."); + GMT_Usage (API, 2, "+n: Closed bodies are expected to be ccw oriented. If not, append this modifier to swap triangle order."); GMT_Usage (API, 1, "\nOR"); - GMT_Usage (API, 1, "\n-M+s/"); + GMT_Usage (API, 1, "\n-M+s,"); GMT_Usage (API, -2, "Select among one or more of the following bodies and append , where x0 and y0 are the horizontal coordinates " "of the body center [default to 0,0], npts is the number of points that a circle is discretized and n_slices " "apply when bodies are made by a pile of slices. For example Spheres and Ellipsoids are made of 2*n_slices and " @@ -239,6 +245,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Option (API, "V"); GMT_Usage (API, 1, "\n-Z"); GMT_Usage (API, -2, "Set z level of reference plane [Default = 0]."); + GMT_Usage (API, 1, "\n-Q"); + GMT_Usage (API, -2, "Use this option if the indices file (-T+v|+f) is 1 based instead of the default (C) 0 based."); GMT_Option (API, "bi"); GMT_Usage (API, 1, "\n-fg Converts geographic grids to meters using a \"Flat Earth\" approximation."); GMT_Option (API, "h,i,o,r"); @@ -326,7 +334,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM } gmt_strtok(opt->arg, "+", &pos, p2); /* Get the string with the model parameters */ if (pos < strlen(opt->arg)) pos--; /* Need to receed 1 due to the (p[0] != '+') test */ - if (!strcmp(&p[2], "bell")) { + if (!strncmp(&p[2], "bell", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[BELL][nBELL][0], &Ctrl->M.params[BELL][nBELL][1], &Ctrl->M.params[BELL][nBELL][2], &Ctrl->M.params[BELL][nBELL][3], &Ctrl->M.params[BELL][nBELL][4], &Ctrl->M.params[BELL][nBELL][5], &Ctrl->M.params[BELL][nBELL][6], &Ctrl->M.params[BELL][nBELL][7], &Ctrl->M.params[BELL][nBELL][8]); if (n_par < 4) err_npar = 1; if (n_par < 7) Ctrl->M.params[BELL][nBELL][6] = Ctrl->n_sigmas; @@ -339,7 +347,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM Ctrl->M.type[BELL][nBELL] = BELL; nBELL++; } - else if (!strcmp(&p[2], "cylinder")) { + else if (!strncmp(&p[2], "cylinder", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[CYLINDER][nCIL][0], &Ctrl->M.params[CYLINDER][nCIL][1], &Ctrl->M.params[CYLINDER][nCIL][2], &Ctrl->M.params[CYLINDER][nCIL][3], &Ctrl->M.params[CYLINDER][nCIL][4], &Ctrl->M.params[CYLINDER][nCIL][5]); if (n_par < 3) err_npar = 1; if (n_par < 6) Ctrl->M.params[CYLINDER][nCIL][5] = Ctrl->npts_circ; @@ -350,14 +358,14 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM Ctrl->M.type[CYLINDER][nCIL] = CYLINDER; nCIL++; } - else if (!strcmp(&p[2], "cone")) { + else if (!strncmp(&p[2], "cone", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[CONE][nCONE][0], &Ctrl->M.params[CONE][nCONE][1], &Ctrl->M.params[CONE][nCONE][2], &Ctrl->M.params[CONE][nCONE][3], &Ctrl->M.params[CONE][nCONE][4]); if (n_par < 4) err_npar = 1; if (n_par == 4) Ctrl->M.params[CONE][nCONE][4] = Ctrl->npts_circ; Ctrl->M.type[CONE][nCONE] = CONE; nCONE++; } - else if (!strcmp(&p[2], "ellipsoid")) { + else if (!strncmp(&p[2], "ellipsoid", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[ELLIPSOID][nELL][0], &Ctrl->M.params[ELLIPSOID][nELL][1], &Ctrl->M.params[ELLIPSOID][nELL][2], &Ctrl->M.params[ELLIPSOID][nELL][3], &Ctrl->M.params[ELLIPSOID][nELL][4], &Ctrl->M.params[ELLIPSOID][nELL][5], &Ctrl->M.params[ELLIPSOID][nELL][6], &Ctrl->M.params[ELLIPSOID][nELL][7]); if (n_par < 4) err_npar = 1; if (n_par < 7) Ctrl->M.params[ELLIPSOID][nELL][6] = Ctrl->npts_circ; @@ -369,19 +377,19 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM Ctrl->M.type[ELLIPSOID][nELL] = ELLIPSOID; nELL++; } - else if (!strcmp(&p[2], "pyramid")) { + else if (!strncmp(&p[2], "pyramid", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[PYRAMID][nPIR][0], &Ctrl->M.params[PYRAMID][nPIR][1], &Ctrl->M.params[PYRAMID][nPIR][2], &Ctrl->M.params[PYRAMID][nPIR][3], &Ctrl->M.params[PYRAMID][nPIR][4], &Ctrl->M.params[PYRAMID][nPIR][5]); if (n_par < 4) err_npar = 1; Ctrl->M.type[PYRAMID][nPIR] = PYRAMID; nPIR++; } - else if (!strcmp(&p[2], "prism")) { + else if (!strncmp(&p[2], "prism", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[PRISM][nPRI][0], &Ctrl->M.params[PRISM][nPRI][1], &Ctrl->M.params[PRISM][nPRI][2], &Ctrl->M.params[PRISM][nPRI][3], &Ctrl->M.params[PRISM][nPRI][4], &Ctrl->M.params[PRISM][nPRI][5]); if (n_par < 4) err_npar = 1; Ctrl->M.type[PRISM][nPRI] = PRISM; nPRI++; } - else if (!strcmp(&p[2], "sphere")) { + else if (!strncmp(&p[2], "sphere", 2U)) { n_par = sscanf (p2, "%lg/%lg/%lg/%lg/%lg/%lg", &Ctrl->M.params[SPHERE][nSPHERE][0], &Ctrl->M.params[SPHERE][nSPHERE][1], &Ctrl->M.params[SPHERE][nSPHERE][2], &Ctrl->M.params[SPHERE][nSPHERE][3], &Ctrl->M.params[SPHERE][nSPHERE][4], &Ctrl->M.params[SPHERE][nSPHERE][5]); if (n_par < 2) err_npar = 1; if (n_par < 5) Ctrl->M.params[SPHERE][nSPHERE][4] = Ctrl->npts_circ; @@ -407,6 +415,9 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM n_errors += gmt_M_repeated_module_option(API, Ctrl->E.active); n_errors += gmt_get_required_double (GMT, opt->arg, opt->option, 0, &Ctrl->E.dz); break; + case 'Q': /* Use 1-based indexing for indices file */ + Ctrl->zero_base = 1; + break; case 'S': n_errors += gmt_M_repeated_module_option(API, Ctrl->S.active); n_errors += gmt_get_required_double (GMT, opt->arg, opt->option, 0, &Ctrl->S.radius); @@ -414,17 +425,21 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM break; case 'T': /* Selected input mesh format */ n_errors += gmt_M_repeated_module_option(API, Ctrl->T.active); - if ((c = strstr(opt->arg, "+v"))) { + if ((c = strstr(opt->arg, "+v")) || (c = strstr(opt->arg, "+f"))) { /* +v is for general triang and +f is for a closed body */ + Ctrl->swap = (strstr(c, "+n") != NULL) ? false : true; c[0] = '\0'; n_errors += gmt_get_required_file(GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->T.t_file)); - Ctrl->T.triangulate = true; + if (c[1] == 'v') Ctrl->T.triangulate = true; + else Ctrl->T.body = true; } else if ((c = strstr(opt->arg, "+r"))) { + Ctrl->swap = (strstr(c, "+n") != NULL) ? false : true; c[0] = '\0'; n_errors += gmt_get_required_file(GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->T.raw_file)); Ctrl->T.raw = true; } - else if ((c = strstr(opt->arg, "+s"))) { + else if ((c = strstr(opt->arg, "+s"))) { /* This option accepts only a true (e.g. not virtual) file name */ + Ctrl->swap = (strstr(c, "+n") != NULL) ? false : true; c[0] = '\0'; n_errors += gmt_get_required_file(GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->T.stl_file)); Ctrl->T.stl = true; @@ -448,6 +463,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM else if (opt->arg[0] == 'v') { Ctrl->T.triangulate = true; Ctrl->T.t_file = strdup(&opt->arg[1]); + Ctrl->swap = true; } else if (opt->arg[0] == 'r') { Ctrl->T.raw_file = strdup(&opt->arg[1]); @@ -481,7 +497,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM n_errors += gmt_M_check_condition(GMT, Ctrl->C.rho == 0.0 && !Ctrl->H.active && !Ctrl->T.m_var4 , "Must specify either -C or -H\n"); (void)gmt_M_check_condition(GMT, Ctrl->G.active && Ctrl->F.active, "Warning: -F overrides -G\n"); - if (gmt_M_check_condition(GMT, Ctrl->T.raw && Ctrl->S.active, "Warning: -Tr overrides -S\n")) + if (gmt_M_check_condition(GMT, Ctrl->T.raw && Ctrl->S.active, "Warning: -T+r overrides -S\n")) Ctrl->S.active = false; return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); @@ -496,9 +512,8 @@ GMT_LOCAL int read_xyz(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, str unsigned int k, n = 0; size_t n_alloc = 10 * GMT_CHUNK; char line[GMT_LEN256] = {""}; - double x1, x2, x3, x4, x5, x6, x7, x8; + double x1, x2, x3; struct GMT_RECORD *In = NULL; - FILE *fp = NULL; if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, 0, GMT_COL_VAR)) != GMT_NOERROR) return error; @@ -630,7 +645,7 @@ GMT_LOCAL int read_vertices(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl int n_skipped, error; uint64_t seg, row, np, k; unsigned int save_col_type[3]; - double d_n = (double)Ctrl->n_triang - 0.5; /* So we can use > in test near line 806 */ + double d_n = (double)Ctrl->n_triang + Ctrl->zero_base - 0.5; /* So we can use > in test near line 806 */ struct GMT_DATASET *Tin = NULL; struct GMT_DATATABLE *T = NULL; @@ -659,9 +674,9 @@ GMT_LOCAL int read_vertices(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl if (T->segment[seg]->data[0][row] > d_n || T->segment[seg]->data[1][row] > d_n || T->segment[seg]->data[2][row] > d_n) n_skipped++; /* Outside point range */ else { - Ctrl->vert[np].a = urint (T->segment[seg]->data[0][row]); - Ctrl->vert[np].b = urint (T->segment[seg]->data[1][row]); - Ctrl->vert[np].c = urint (T->segment[seg]->data[2][row]); + Ctrl->vert[np].a = urint(T->segment[seg]->data[0][row]) - Ctrl->zero_base; + Ctrl->vert[np].b = urint(T->segment[seg]->data[1][row]) - Ctrl->zero_base; + Ctrl->vert[np].c = urint(T->segment[seg]->data[2][row]) - Ctrl->zero_base; np++; } } @@ -679,7 +694,7 @@ GMT_LOCAL int read_vertices(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl /* -----------------------------------------------------------------*/ GMT_LOCAL int read_raw(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { - /* read a file with triangles in the raw format and returns nb of triangles */ + /* read a file with triangles in the raw format */ unsigned int row, seg, nt; struct GMT_DATASET *In = NULL; @@ -716,6 +731,31 @@ GMT_LOCAL int read_raw(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { return 0; } +/* -----------------------------------------------------------------*/ +GMT_LOCAL void closed_body(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { + /* Use the xyz and indices data read with read_xyz and read_raw to create a closed body */ + unsigned int row, nf = 0; + + Ctrl->raw_mesh = gmt_M_memory (GMT, NULL, Ctrl->n_vert, struct GMTGRAVMAG3D_RAW); /* n_vert is a bad name as it's = n_faces */ + + for (row = 0; row < Ctrl->n_vert; row++) { + Ctrl->raw_mesh[nf].t1[0] = Ctrl->triang[Ctrl->vert[nf].a].x; + Ctrl->raw_mesh[nf].t1[1] = Ctrl->triang[Ctrl->vert[nf].a].y; + Ctrl->raw_mesh[nf].t1[2] = Ctrl->triang[Ctrl->vert[nf].a].z; + + Ctrl->raw_mesh[nf].t2[0] = Ctrl->triang[Ctrl->vert[nf].b].x; + Ctrl->raw_mesh[nf].t2[1] = Ctrl->triang[Ctrl->vert[nf].b].y; + Ctrl->raw_mesh[nf].t2[2] = Ctrl->triang[Ctrl->vert[nf].b].z; + + Ctrl->raw_mesh[nf].t3[0] = Ctrl->triang[Ctrl->vert[nf].c].x; + Ctrl->raw_mesh[nf].t3[1] = Ctrl->triang[Ctrl->vert[nf].c].y; + Ctrl->raw_mesh[nf].t3[2] = Ctrl->triang[Ctrl->vert[nf].c].z; + nf++; + } + + Ctrl->n_raw_triang = Ctrl->n_vert; +} + #include "solids.c" /* -----------------------------------------------------------------*/ GMT_LOCAL void solids(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { @@ -842,29 +882,35 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { if ((error = read_xyz(GMT, Ctrl, options, &lon_0, &lat_0))) Return (error); - /* read vertex file */ - if ((error = read_vertices(GMT, Ctrl))) + if ((error = read_vertices(GMT, Ctrl))) /* read vertex file */ Return (error); vert = Ctrl->vert; - Ctrl->t_center = gmt_M_memory (GMT, NULL, Ctrl->n_vert, struct GMTGRAVMAG3D_XYZ); - /* compute approximate center of each triangle */ - n_swap = check_triang_cw (Ctrl, Ctrl->n_vert, 0); - set_center (Ctrl); + if (Ctrl->swap) n_swap = check_triang_cw(Ctrl, Ctrl->n_vert, 0); + Ctrl->t_center = gmt_M_memory(GMT, NULL, Ctrl->n_vert, struct GMTGRAVMAG3D_XYZ); + set_center(Ctrl); /* compute approximate center of each triangle */ } else if (Ctrl->T.stl) { /* Read STL file defining a closed volume */ if ((ndata_s = read_stl(GMT, Ctrl)) < 0) { GMT_Report (API, GMT_MSG_ERROR, "Cannot open file %s\n", Ctrl->T.stl_file); Return (GMT_ERROR_ON_FOPEN); } - /*n_swap = check_triang_cw (ndata_s, 1);*/ + if (Ctrl->swap) swap_triang(Ctrl, ndata_s); } else if (Ctrl->T.raw) { /* Read RAW file defining a closed volume */ if ((error = read_raw(GMT, Ctrl))) Return (error); + if (Ctrl->swap) swap_triang(Ctrl, Ctrl->n_raw_triang); + } + else if (Ctrl->T.body) { /* Read FV files defining a closed volume */ + if ((error = read_xyz(GMT, Ctrl, options, &lon_0, &lat_0))) + Return (error); + if ((error = read_vertices(GMT, Ctrl))) /* read vertex file */ + Return (error); - /*n_swap = check_triang_cw (Ctrl, Ctrl->n_raw_triang, 1);*/ + closed_body(GMT, Ctrl); + if (Ctrl->swap) swap_triang(Ctrl, Ctrl->n_vert); } else if (Ctrl->M.active) { solids(GMT, Ctrl); @@ -875,8 +921,8 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { /* --------------------------------------------------------------------------------------- */ if (Ctrl->G.active) { - if ((Gout = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, Ctrl->I.inc, \ - GMT_GRID_DEFAULT_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error); + if ((Gout = GMT_Create_Data(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, Ctrl->I.inc, \ + GMT_GRID_DEFAULT_REG, 0, NULL)) == NULL) Return (API->error); GMT_Report (API, GMT_MSG_INFORMATION, "Grid dimensions are n_columns = %d, n_rows = %d\n", Gout->header->n_columns, Gout->header->n_rows); @@ -928,8 +974,8 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { if (Ctrl->C.active) body_desc.n_f = 2; /* Number of prism facets that count */ } - else if (Ctrl->T.raw || Ctrl->T.stl || Ctrl->M.active) { - n_triang = (Ctrl->T.raw || Ctrl->M.active) ? Ctrl->n_raw_triang : ndata_s; + else if (Ctrl->T.raw || Ctrl->T.stl || Ctrl->T.body || Ctrl->M.active) { + n_triang = (Ctrl->T.raw || Ctrl->T.body || Ctrl->M.active) ? Ctrl->n_raw_triang : ndata_s; body_desc.n_f = 1; body_desc.n_v = gmt_M_memory (GMT, NULL, body_desc.n_f, unsigned int); body_desc.n_v[0] = 3; @@ -1017,7 +1063,7 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { /* ---------------> Now start computing <------------------------------------- */ gmt_M_tic(GMT); - one_100 = (gmt_grdfloat)(n_triang / 100.); + one_100 = (gmt_grdfloat)(n_triang / 100.0 * 2); /* Make it 2% */ s_rad2 = Ctrl->S.radius*Ctrl->S.radius; if (Ctrl->G.active) { /* Compute the cos(lat) vector only once */ @@ -1040,7 +1086,7 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { continue; if (Ctrl->T.triangulate) z_th = facet_triangulate (Ctrl, body_verts, i, bat); - else if (Ctrl->T.raw || Ctrl->T.stl || Ctrl->M.active) + else if (Ctrl->T.raw || Ctrl->T.stl || Ctrl->T.body || Ctrl->M.active) z_th = facet_raw (Ctrl, body_verts, i, Ctrl->box.is_geog); if (z_th) { if (Ctrl->G.active) { /* grid */ @@ -1056,7 +1102,7 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { if (!DO) continue; } - a = okabe (GMT, x_o, y_o, Ctrl->L.zobs, Ctrl->C.rho, Ctrl->C.active, body_desc, body_verts, km, pm, loc_or, okabe_mag_param, okabe_mag_var); + a = okabe(GMT, x_o, y_o, Ctrl->L.zobs, Ctrl->C.rho, Ctrl->C.active, body_desc, body_verts, km, pm, loc_or, okabe_mag_param, okabe_mag_var); Gout->data[ij] += (gmt_grdfloat)a; } } @@ -1069,8 +1115,8 @@ EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { DO = (DX*DX + DY*DY) < s_rad2; if (!DO) continue; } - a = okabe (GMT, x_obs[kk], y_obs[kk], Ctrl->L.zobs, Ctrl->C.rho, Ctrl->C.active, body_desc, body_verts, km, pm, loc_or, okabe_mag_param, okabe_mag_var); - g[kk] += (gmt_grdfloat)a; + a = okabe(GMT, x_obs[kk], y_obs[kk], Ctrl->L.zobs, Ctrl->C.rho, Ctrl->C.active, body_desc, body_verts, km, pm, loc_or, okabe_mag_param, okabe_mag_var); + g[kk] += a; } } } @@ -1390,15 +1436,19 @@ GMT_LOCAL int check_triang_cw(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, un x1 = Ctrl->triang[Ctrl->vert[i].a].x; y1 = Ctrl->triang[Ctrl->vert[i].a].y; x2 = Ctrl->triang[Ctrl->vert[i].b].x; y2 = Ctrl->triang[Ctrl->vert[i].b].y; x3 = Ctrl->triang[Ctrl->vert[i].c].x; y3 = Ctrl->triang[Ctrl->vert[i].c].y; + det = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1); } - else if (type == 1) { /* raw */ - x1 = Ctrl->raw_mesh[i].t1[0]; y1 = Ctrl->raw_mesh[i].t1[1]; - x2 = Ctrl->raw_mesh[i].t2[0]; y2 = Ctrl->raw_mesh[i].t2[1]; - x3 = Ctrl->raw_mesh[i].t3[0]; y3 = Ctrl->raw_mesh[i].t3[1]; + else if (type == 1) { /* Closed bodies. This is probably never used by now. */ + double AB[2], AC[2]; + AB[0] = Ctrl->raw_mesh[i].t2[0] - Ctrl->raw_mesh[i].t1[0]; + AB[1] = Ctrl->raw_mesh[i].t2[1] - Ctrl->raw_mesh[i].t1[1]; + //AB[2] = Ctrl->raw_mesh[i].t2[2] - Ctrl->raw_mesh[i].t1[2]; /* Not needed */ + AC[0] = Ctrl->raw_mesh[i].t3[0] - Ctrl->raw_mesh[i].t1[0]; + AC[1] = Ctrl->raw_mesh[i].t3[1] - Ctrl->raw_mesh[i].t1[1]; + //AC[2] = Ctrl->raw_mesh[i].t3[2] - Ctrl->raw_mesh[i].t1[2]; + det = (AB[0]*AC[1] - AB[1]*AC[0]); /* Not the det, the normal z component*/ } - det = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1); - if (det < 0.0) { /* counter clockwise triangle -> swap vertex order */ if (type == 0) { tmp = Ctrl->vert[i].b; @@ -1407,12 +1457,12 @@ GMT_LOCAL int check_triang_cw(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, un n_swaped++; } else if (type == 1) { - d_tmp[0] = Ctrl->raw_mesh[i].t2[0]; - d_tmp[1] = Ctrl->raw_mesh[i].t2[1]; - d_tmp[2] = Ctrl->raw_mesh[i].t2[2]; - Ctrl->raw_mesh[i].t2[0] = Ctrl->raw_mesh[i].t3[0]; - Ctrl->raw_mesh[i].t2[1] = Ctrl->raw_mesh[i].t3[1]; - Ctrl->raw_mesh[i].t2[2] = Ctrl->raw_mesh[i].t3[2]; + d_tmp[0] = Ctrl->raw_mesh[i].t1[0]; + d_tmp[1] = Ctrl->raw_mesh[i].t1[1]; + d_tmp[2] = Ctrl->raw_mesh[i].t1[2]; + Ctrl->raw_mesh[i].t1[0] = Ctrl->raw_mesh[i].t3[0]; + Ctrl->raw_mesh[i].t1[1] = Ctrl->raw_mesh[i].t3[1]; + Ctrl->raw_mesh[i].t1[2] = Ctrl->raw_mesh[i].t3[2]; Ctrl->raw_mesh[i].t3[0] = d_tmp[0]; Ctrl->raw_mesh[i].t3[1] = d_tmp[1]; Ctrl->raw_mesh[i].t3[2] = d_tmp[2]; @@ -1422,3 +1472,20 @@ GMT_LOCAL int check_triang_cw(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, un } return n_swaped; } + +GMT_LOCAL void swap_triang(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n_tri) { + /* Swaps vertex order in triangulation */ + unsigned int i; + double d_tmp[3]; + for (i = 0; i < n_tri; i++) { + d_tmp[0] = Ctrl->raw_mesh[i].t1[0]; + d_tmp[1] = Ctrl->raw_mesh[i].t1[1]; + d_tmp[2] = Ctrl->raw_mesh[i].t1[2]; + Ctrl->raw_mesh[i].t1[0] = Ctrl->raw_mesh[i].t3[0]; + Ctrl->raw_mesh[i].t1[1] = Ctrl->raw_mesh[i].t3[1]; + Ctrl->raw_mesh[i].t1[2] = Ctrl->raw_mesh[i].t3[2]; + Ctrl->raw_mesh[i].t3[0] = d_tmp[0]; + Ctrl->raw_mesh[i].t3[1] = d_tmp[1]; + Ctrl->raw_mesh[i].t3[2] = d_tmp[2]; + } +} diff --git a/src/potential/solids.c b/src/potential/solids.c index 4a2f43b6fbb..93323f96a02 100644 --- a/src/potential/solids.c +++ b/src/potential/solids.c @@ -151,7 +151,12 @@ int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_ty a = b; b = c; c = Ctrl->M.params[body_type][nb][0]; } - z_top = z_c + c; z_bot = z_c; + if (body_type == SPHERE || body_type == ELLIPSOID) { + z_top = z_c + c; z_bot = z_c - c; + } + else { + z_top = z_c + c; z_bot = z_c; + } n_tri = (hemi) ? 2 * npts_circ * n_slices : 2 * (npts_circ * (n_slices*2 - 1)); ellipse[0] = (struct GRAVMAG_XY *) calloc((size_t) (npts_circ+1), sizeof(struct GRAVMAG_XY)); @@ -165,7 +170,7 @@ int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_ty half_width_x = 0.5 * a; dx = n_sigmas * a / n_slices; dy = n_sigmas * b / n_slices; - for (j = 0; j < Ctrl->n_slices; j++) { + for (j = 0; j < n_slices; j++) { j1 = j + 1; if (cone || piram) { ai0 = j * a / n_slices; bi0 = j * b / n_slices; @@ -187,6 +192,7 @@ int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_ty ai1 = a*cos(M_PI_2-j1*d_tet); bi1 = b*cos(M_PI_2-j1*d_tet); zi1 = z_top - c * (1. - sqrt(1. - (ai1/a)*(ai1/a))); } + for (i = 0; i < npts_circ; i++) { /* compute slice j */ ellipse[0][i].x = x0 + ai0 * cos (i*dfi); ellipse[0][i].y = y0 + bi0 * sin (i*dfi); @@ -218,6 +224,7 @@ int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_ty } first = false; } + /* First half is ready. Now, either close it and return or construct the other half by simetry */ if (cone || piram || sino || hemi) { /* close the base and return */ if (sino && fabs ((zi1 - z_c) / z_c) > 0.01) { @@ -245,19 +252,21 @@ int five_psoid(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, int body_ty free((void *)ellipse[1]); return (k); } + n_tri = npts_circ * (n_slices*2 - 1); - for (j = n_tri-1, i = n_tri; j >= 0; j--, i++) { - Ctrl->raw_mesh[i+i_tri].t1[0] = Ctrl->raw_mesh[j+i_tri].t1[0]; - Ctrl->raw_mesh[i+i_tri].t1[1] = Ctrl->raw_mesh[j+i_tri].t1[1]; - Ctrl->raw_mesh[i+i_tri].t1[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t1[2]); + z_c *= -1; /* It was still positive up */ + for (j = n_tri-1, i = n_tri; j >= 0; j--, i++) { /* Have to use the sequence t3,t2,t1 here to keep the body ccw */ + Ctrl->raw_mesh[i+i_tri].t3[0] = Ctrl->raw_mesh[j+i_tri].t1[0]; + Ctrl->raw_mesh[i+i_tri].t3[1] = Ctrl->raw_mesh[j+i_tri].t1[1]; + Ctrl->raw_mesh[i+i_tri].t3[2] = (z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t1[2])); Ctrl->raw_mesh[i+i_tri].t2[0] = Ctrl->raw_mesh[j+i_tri].t2[0]; Ctrl->raw_mesh[i+i_tri].t2[1] = Ctrl->raw_mesh[j+i_tri].t2[1]; - Ctrl->raw_mesh[i+i_tri].t2[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t2[2]); + Ctrl->raw_mesh[i+i_tri].t2[2] = (z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t2[2])); - Ctrl->raw_mesh[i+i_tri].t3[0] = Ctrl->raw_mesh[j+i_tri].t3[0]; - Ctrl->raw_mesh[i+i_tri].t3[1] = Ctrl->raw_mesh[j+i_tri].t3[1]; - Ctrl->raw_mesh[i+i_tri].t3[2] = z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t3[2]); + Ctrl->raw_mesh[i+i_tri].t1[0] = Ctrl->raw_mesh[j+i_tri].t3[0]; + Ctrl->raw_mesh[i+i_tri].t1[1] = Ctrl->raw_mesh[j+i_tri].t3[1]; + Ctrl->raw_mesh[i+i_tri].t1[2] = (z_c + (z_c - Ctrl->raw_mesh[j+i_tri].t3[2])); } free ((void *)ellipse[0]); free ((void *)ellipse[1]); diff --git a/test/potential/spheres.sh b/test/potential/spheres.sh index 722cfee60d8..9f96c8977b8 100755 --- a/test/potential/spheres.sh +++ b/test/potential/spheres.sh @@ -10,7 +10,7 @@ echo -50 0 > li echo 50 0 >> li gmt sample1d li -Fl -I1 > li1.dat -gmt gmtgravmag3d -Tr@sphere_raw.txt -C$ro -Fli1.dat > ptodos_g.dat +gmt gmtgravmag3d -T@sphere_raw.txt+r+n -C$ro -Fli1.dat > ptodos_g.dat # xyzokb solution gmt psxy ptodos_g.dat -i0,2 -R-50/50/0/0.125 -JX14c/10c -Bx10f5 -By.01 -BWSne+t"Anomaly (mGal)" -W1p -P -K > $ps