Skip to content

Commit

Permalink
Expand gmtgravmag3d to accept input from virtual files.
Browse files Browse the repository at this point in the history
Fixes for the newly introduced -T syntax.

Add option to force triangle handedness order swapping.

Fix bugs in the sphere and ellipsoidal bodies generation.
  • Loading branch information
joa-quim committed Jan 30, 2025
1 parent 9f10ccc commit aef2f7a
Show file tree
Hide file tree
Showing 6 changed files with 161 additions and 72 deletions.
27 changes: 20 additions & 7 deletions doc/rst/source/supplements/potential/gmtgravmag3d.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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| ]
Expand Down Expand Up @@ -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 </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;
Expand All @@ -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
------------------

Expand All @@ -119,28 +125,35 @@ 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.

.. _-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.
Expand Down
4 changes: 2 additions & 2 deletions src/gmt_api.c
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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 */
Expand Down
2 changes: 1 addition & 1 deletion src/gmt_grdio.c
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
Loading

0 comments on commit aef2f7a

Please sign in to comment.