From 504e6d226ba1da344137d610a5256042b7cd2638 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 22 Jan 2025 22:31:04 +0000 Subject: [PATCH] Several fixes to let gmtgravmag3d pass in memory data via modified -T option. (#8680) * Modify option -T to use modifiers. Let both input xyz file and vertices (new -T+v) be readable by externals. * Save number of columns read when requesrec reading a file with unknown n cols. * Change KEYs to allow -F from externals. Change syntax of -T option * Update the -T docs. --- .../supplements/potential/gmtgravmag3d.rst | 10 +- src/gmt_api.c | 11 +- src/potential/gmtgravmag3d.c | 207 +++++++++--------- 3 files changed, 120 insertions(+), 108 deletions(-) diff --git a/doc/rst/source/supplements/potential/gmtgravmag3d.rst b/doc/rst/source/supplements/potential/gmtgravmag3d.rst index b34624b10b4..3c6d50b8d12 100644 --- a/doc/rst/source/supplements/potential/gmtgravmag3d.rst +++ b/doc/rst/source/supplements/potential/gmtgravmag3d.rst @@ -12,7 +12,7 @@ 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|\ **+v**\ *vert_file* OR |-T|\ **+r\|+s**\ *raw_file* OR |-M|\ **+s**\ *body,params* [ |-C|\ *density* ] [ |-E|\ *thickness* ] [ |-F|\ *xy_file* ] @@ -97,15 +97,15 @@ Required Arguments (not all) .. _-T: -**-Tv**\ *vert_file* (must have when passing a *xyz_file*) OR **-Tr\|s**\ *raw_file* - Gives names of a xyz and vertex (**-Tv**\ *vert_file*) files defining a close surface. +**-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. 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; 5 columns: mag, mag dip; 6 columns: mag, mag dec, mag dip; 8 columns: field dec, field dip, mag, mag dec, mag dip. - When n columns > 3 the third argument of the |-H| option is ignored. A *raw* format (selected by the **-Tr** option) + When n columns > 3 the third argument of the |-H| option is ignored. A *raw* format (selected with the **-T+r** option) is a file with N rows (one per triangle) and 9 columns corresponding to the *x, y, z* coordinates of each of the three - vertex of each triangle. Alternatively, the **-Ts** option indicates that the surface file is in the ASCII STL + 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. Optional Arguments diff --git a/src/gmt_api.c b/src/gmt_api.c index 9afb1436385..99e5e93505b 100644 --- a/src/gmt_api.c +++ b/src/gmt_api.c @@ -10315,7 +10315,7 @@ GMT_LOCAL void gmtapi_get_record_init (struct GMTAPI_CTRL *API) { } } -void * GMT_Get_Record (void *V_API, unsigned int mode, int *retval) { +void *GMT_Get_Record(void *V_API, unsigned int mode, int *retval) { /* Retrieves the next data record from the virtual input source and * returns the number of columns found via *retval (unless retval == NULL). * If current record is a segment header then we return 0. @@ -10348,7 +10348,10 @@ void * GMT_Get_Record (void *V_API, unsigned int mode, int *retval) { record = API->api_get_record (API, mode, &n_fields); } while (API->get_next_record); - if (!(n_fields == EOF || n_fields == GMT_IO_NEXT_FILE)) API->current_rec[GMT_IN]++; /* Increase record count, unless EOF */ + if (!(n_fields == EOF || n_fields == GMT_IO_NEXT_FILE)) { /* Increase record count, unless EOF */ + API->current_rec[GMT_IN]++; + if (GMT->current.io.variable_in_columns) GMT->current.io.n_numerical_cols = (unsigned int)n_fields; /* Keep track of this */ + } if (retval) *retval = n_fields; /* Requested we return the number of fields found */ return (record); /* Return pointer to current record */ @@ -13461,6 +13464,10 @@ struct GMT_RESOURCE * GMT_Encode_Options (void *V_API, const char *module_name, 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 */ } + /* 1z. Check if gmtgravmag3d is producing grids or datasets */ + else if (!strncmp (module, "gmtgravmag3d", 12U)) { + //type = (opt = GMT_Find_Option (API, 'F', *head)) ? 'D' : 'G'; /* Giving -F means compute over a line, else grid */ + } /* 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., "D)" +#define THIS_MODULE_KEYS "DF" #define THIS_MODULE_NEEDS "R" #define THIS_MODULE_OPTIONS "-:RVfhior" GMT_ADD_xg_OPT @@ -178,7 +178,7 @@ GMT_LOCAL int check_triang_cw (struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, u 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 [] -Tv | -Tr|s | -M+s/ [-C] [-E] " + GMT_Usage (API, 0, "usage: %s [
] -T+v | -T+r|+s | -M+s/ [-C] [-E] " "[-F] [-G%s] [-H///] [%s] [-L] [%s] " "[-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,13 +188,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 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-Tv | -Tr|s"); - GMT_Usage (API, -2, "Give names of xyz and vertex (-Tv) files defining a closed surface. " - "If has more then 3 columns it means variable magnetization; see docs for more details. " + 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. " + "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, "r: Append in raw triangle format (x1 y1 z1 x2 ... z3)."); - GMT_Usage (API, 3, "s: Append in STL format."); + 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, 1, "\nOR"); 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 " @@ -263,7 +263,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM unsigned int pos = 0, n_errors = 0; int n_par, err_npar = 0, nBELL = 0, nCIL = 0, nPRI = 0, nCONE = 0, nELL = 0, nPIR = 0, nSPHERE = 0; - char p[GMT_LEN256] = {""}, p2[GMT_LEN256] = {""}; + char p[GMT_LEN256] = {""}, p2[GMT_LEN256] = {""}, *c = NULL; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; @@ -404,17 +404,32 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM } break; case 'E': - n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); + 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 'S': - n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active); + 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); Ctrl->S.radius *= 1000; /* Expect km, convert to meters */ break; case 'T': /* Selected input mesh format */ - n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active); - if (opt->arg[0] == 'p') { + n_errors += gmt_M_repeated_module_option(API, Ctrl->T.active); + if ((c = strstr(opt->arg, "+v"))) { + 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; + } + else if ((c = strstr(opt->arg, "+r"))) { + 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"))) { + 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; + } + else if (opt->arg[0] == 'p') { /* Non documented. And start of backward compat section. */ char *pch; Ctrl->T.xyz_file = strdup(&opt->arg[1]); Ctrl->T.triangulate = true; @@ -469,8 +484,6 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM if (gmt_M_check_condition(GMT, Ctrl->T.raw && Ctrl->S.active, "Warning: -Tr overrides -S\n")) Ctrl->S.active = false; - /*n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input file\n");*/ - return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR); } @@ -478,6 +491,7 @@ static int parse(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GM /* -------------------------------------------------------------------------*/ GMT_LOCAL int read_xyz(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, struct GMT_OPTION *options, double *lon_0, double *lat_0) { /* read xyz[m] file with point data coordinates */ + bool first_time = true; int n_cols = 0, error; unsigned int k, n = 0; size_t n_alloc = 10 * GMT_CHUNK; @@ -486,20 +500,7 @@ GMT_LOCAL int read_xyz(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, str struct GMT_RECORD *In = NULL; FILE *fp = NULL; - /* First, count number of columns */ - if ((fp = gmt_fopen (GMT, Ctrl->T.xyz_file, "r")) == NULL) return -1; - while (fgets (line, GMT_LEN256, fp)) { - if (line[0] == '#') continue; - n_cols = sscanf (line, "%lg %lg %lg %lg %lg %lg %lg %lg", &x1, &x2, &x3, &x4, &x5, &x6, &x7, &x8); - break; - } - fclose(fp); - if (n_cols < 3 || n_cols == 7 || n_cols > 8) { - GMT_Report (GMT->parent, GMT_MSG_ERROR, "Wrong number of columns (%d) in file %s\n", n_cols, Ctrl->T.xyz_file); - return -1; - } - - if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, (unsigned int)n_cols, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) + if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, 0, GMT_COL_VAR)) != GMT_NOERROR) return error; if (GMT_Init_IO (GMT->parent, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) /* Establishes data input */ return GMT->parent->error; @@ -507,55 +508,65 @@ GMT_LOCAL int read_xyz(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, str if (GMT_Begin_IO (GMT->parent, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) /* Enables data input and sets access mode */ return GMT->parent->error; - Ctrl->triang = gmt_M_memory (GMT, NULL, n_alloc, struct GMTGRAVMAG3D_XYZ); - Ctrl->T.m_var = (n_cols == 3) ? false : true; /* x,y,z */ - if (n_cols == 4) { - Ctrl->T.m_var1 = true; - Ctrl->box.mag_int = gmt_M_memory (GMT, NULL, n_alloc, double); - } - else if (n_cols == 5) { - Ctrl->T.m_var2 = true; - Ctrl->okabe_mag_var2 = gmt_M_memory (GMT, NULL, n_alloc, struct MAG_VAR2); - } - else if (n_cols == 6) { - Ctrl->T.m_var3 = true; - Ctrl->okabe_mag_var3 = gmt_M_memory (GMT, NULL, n_alloc, struct MAG_VAR3); - } - else if (n_cols == 8) { - Ctrl->T.m_var4 = true; - Ctrl->okabe_mag_var4 = gmt_M_memory (GMT, NULL, n_alloc, struct MAG_VAR4); - } - - if (n_cols > 3) { /* A bit ugly doing this here but only now we know enough */ - Ctrl->H.active = true; - Ctrl->C.active = false; - } - do { /* Keep returning records until we reach EOF */ - if ((In = GMT_Get_Record (GMT->parent, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */ - if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */ + if ((In = GMT_Get_Record(GMT->parent, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */ + if (gmt_M_rec_is_error(GMT)) /* Bail if there are any read errors */ return (GMT_RUNTIME_ERROR); - else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */ + else if (gmt_M_rec_is_eof(GMT)) /* Reached end of file */ break; continue; /* Go back and read the next record */ } if (In->data == NULL) { - gmt_quit_bad_record (GMT->parent, In); + gmt_quit_bad_record(GMT->parent, In); return (GMT->parent->error); } + if (first_time) { /* Only now we know the number of columns and can do some checks and memory allocs. */ + n_cols = GMT->current.io.n_numerical_cols; + + if (n_cols < 3 || n_cols == 7 || n_cols > 8) { + GMT_Report(GMT->parent, GMT_MSG_ERROR, "Wrong number of columns (%d) in file %s\n", n_cols, Ctrl->T.xyz_file); + return (GMT_RUNTIME_ERROR); + } + + Ctrl->triang = gmt_M_memory(GMT, NULL, n_alloc, struct GMTGRAVMAG3D_XYZ); + Ctrl->T.m_var = (n_cols == 3) ? false : true; /* x,y,z */ + if (n_cols == 4) { + Ctrl->T.m_var1 = true; + Ctrl->box.mag_int = gmt_M_memory(GMT, NULL, n_alloc, double); + } + else if (n_cols == 5) { + Ctrl->T.m_var2 = true; + Ctrl->okabe_mag_var2 = gmt_M_memory(GMT, NULL, n_alloc, struct MAG_VAR2); + } + else if (n_cols == 6) { + Ctrl->T.m_var3 = true; + Ctrl->okabe_mag_var3 = gmt_M_memory(GMT, NULL, n_alloc, struct MAG_VAR3); + } + else if (n_cols == 8) { + Ctrl->T.m_var4 = true; + Ctrl->okabe_mag_var4 = gmt_M_memory(GMT, NULL, n_alloc, struct MAG_VAR4); + } + + if (n_cols > 3) { + Ctrl->H.active = true; + Ctrl->C.active = false; + } + first_time = false; + } + if (n == n_alloc) { n_alloc = (size_t)(n_alloc * 1.7); - Ctrl->triang = gmt_M_memory (GMT, Ctrl->triang, n_alloc, struct GMTGRAVMAG3D_XYZ); + Ctrl->triang = gmt_M_memory(GMT, Ctrl->triang, n_alloc, struct GMTGRAVMAG3D_XYZ); if (Ctrl->T.m_var1) - Ctrl->box.mag_int = gmt_M_memory (GMT, Ctrl->box.mag_int, n_alloc, double); + Ctrl->box.mag_int = gmt_M_memory(GMT, Ctrl->box.mag_int, n_alloc, double); else if (Ctrl->T.m_var2) - Ctrl->okabe_mag_var2 = gmt_M_memory (GMT, Ctrl->okabe_mag_var2, n_alloc, struct MAG_VAR2); + Ctrl->okabe_mag_var2 = gmt_M_memory(GMT, Ctrl->okabe_mag_var2, n_alloc, struct MAG_VAR2); else if (Ctrl->T.m_var3) - Ctrl->okabe_mag_var3 = gmt_M_memory (GMT, Ctrl->okabe_mag_var3, n_alloc, struct MAG_VAR3); + Ctrl->okabe_mag_var3 = gmt_M_memory(GMT, Ctrl->okabe_mag_var3, n_alloc, struct MAG_VAR3); else - Ctrl->okabe_mag_var4 = gmt_M_memory (GMT, Ctrl->okabe_mag_var4, n_alloc, struct MAG_VAR4); + Ctrl->okabe_mag_var4 = gmt_M_memory(GMT, Ctrl->okabe_mag_var4, n_alloc, struct MAG_VAR4); } Ctrl->triang[n].x = In->data[0]; Ctrl->triang[n].y = -In->data[1]; /* - because y must be positive 'south'*/ @@ -581,14 +592,16 @@ GMT_LOCAL int read_xyz(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl, str n++; } while (true); - if (GMT_End_IO (GMT->parent, GMT_IN, 0) != GMT_NOERROR) /* Disables further data input */ + if (GMT->current.io.variable_in_columns) GMT->current.io.n_numerical_cols = 0; /* To undo what we did in GMT_Get_Record() */ + + if (GMT_End_IO(GMT->parent, GMT_IN, 0) != GMT_NOERROR) /* Disables further data input */ return GMT->parent->error; - Ctrl->triang = gmt_M_memory (GMT, Ctrl->triang, (size_t)n, struct GMTGRAVMAG3D_XYZ); - if (Ctrl->T.m_var1) Ctrl->box.mag_int = gmt_M_memory (GMT, Ctrl->box.mag_int, (size_t)n, double); - else if (Ctrl->T.m_var2) Ctrl->okabe_mag_var2 = gmt_M_memory (GMT, Ctrl->okabe_mag_var2, (size_t)n, struct MAG_VAR2); - else if (Ctrl->T.m_var3) Ctrl->okabe_mag_var3 = gmt_M_memory (GMT, Ctrl->okabe_mag_var3, (size_t)n, struct MAG_VAR3); - else Ctrl->okabe_mag_var4 = gmt_M_memory (GMT, Ctrl->okabe_mag_var4, (size_t)n, struct MAG_VAR4); + Ctrl->triang = gmt_M_memory(GMT, Ctrl->triang, (size_t)n, struct GMTGRAVMAG3D_XYZ); + if (Ctrl->T.m_var1) Ctrl->box.mag_int = gmt_M_memory(GMT, Ctrl->box.mag_int, (size_t)n, double); + else if (Ctrl->T.m_var2) Ctrl->okabe_mag_var2 = gmt_M_memory(GMT, Ctrl->okabe_mag_var2, (size_t)n, struct MAG_VAR2); + else if (Ctrl->T.m_var3) Ctrl->okabe_mag_var3 = gmt_M_memory(GMT, Ctrl->okabe_mag_var3, (size_t)n, struct MAG_VAR3); + else Ctrl->okabe_mag_var4 = gmt_M_memory(GMT, Ctrl->okabe_mag_var4, (size_t)n, struct MAG_VAR4); *lon_0 = 0.; *lat_0 = 0.; if (Ctrl->box.is_geog) { @@ -745,7 +758,7 @@ GMT_LOCAL void solids(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { #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_gmtgravmag3d (void *V_API, int mode, void *args) { +EXTERN_MSC int GMT_gmtgravmag3d(void *V_API, int mode, void *args) { bool bat = true, DO = true; unsigned int row, col, i, j, k, kk; @@ -857,14 +870,6 @@ EXTERN_MSC int GMT_gmtgravmag3d (void *V_API, int mode, void *args) { solids(GMT, Ctrl); } -#if 0 - for (i = 0; i < 24; i++) { - fprintf(stderr, "%.2f %.2f %.2f ", Ctrl->raw_mesh[i].t1[0], Ctrl->raw_mesh[i].t1[1], Ctrl->raw_mesh[i].t1[2]); - fprintf(stderr, "%.2f %.2f %.2f ", Ctrl->raw_mesh[i].t2[0], Ctrl->raw_mesh[i].t2[1], Ctrl->raw_mesh[i].t2[2]); - fprintf(stderr, "%.2f %.2f %.2f\n", Ctrl->raw_mesh[i].t3[0], Ctrl->raw_mesh[i].t3[1], Ctrl->raw_mesh[i].t3[2]); - } -#endif - if (n_swap > 0) GMT_Report (API, GMT_MSG_INFORMATION, "%d triangles had ccw order\n", n_swap); /* --------------------------------------------------------------------------------------- */ @@ -1071,26 +1076,7 @@ EXTERN_MSC int GMT_gmtgravmag3d (void *V_API, int mode, void *args) { } } - if (Ctrl->G.active) { - if (Ctrl->C.active) { - strcpy (Gout->header->title, "Gravity field"); - strcpy (Gout->header->z_units, "mGal"); - } - else { - strcpy (Gout->header->title, "Magnetic field"); - strcpy (Gout->header->z_units, "nT"); - } - - if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Gout)) { - error = API->error; - goto END; - } - if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Gout) != GMT_NOERROR) { - error = API->error; - goto END; - } - } - else { + if (Ctrl->F.active) { double out[3]; char save[GMT_LEN64] = {""}; struct GMT_RECORD *Out = gmt_new_record (GMT, out, NULL); @@ -1125,6 +1111,25 @@ EXTERN_MSC int GMT_gmtgravmag3d (void *V_API, int mode, void *args) { goto END; } } + else { + if (Ctrl->C.active) { + strcpy(Gout->header->title, "Gravity field"); + strcpy(Gout->header->z_units, "mGal"); + } + else { + strcpy(Gout->header->title, "Magnetic field"); + strcpy(Gout->header->z_units, "nT"); + } + + if (GMT_Set_Comment(API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Gout)) { + error = API->error; + goto END; + } + if (GMT_Write_Data(API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Gout) != GMT_NOERROR) { + error = API->error; + goto END; + } + } END: gmt_M_toc(GMT,""); /* Print total run time, but only if -Vt was set */ @@ -1154,7 +1159,7 @@ EXTERN_MSC int GMT_gmtgravmag3d (void *V_API, int mode, void *args) { } /* -----------------------------------------------------------------*/ -GMT_LOCAL int read_stl (struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { +GMT_LOCAL int read_stl(struct GMT_CTRL *GMT, struct GMTGRAVMAG3D_CTRL *Ctrl) { /* read a file with triagles in the stl format and returns nb of triangles */ unsigned int ndata_s; size_t n_alloc; @@ -1206,7 +1211,7 @@ GMT_LOCAL int read_stl (struct GMT_CTRL *GMT, 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_triangulate(struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body_verts, unsigned int i, bool bat) { /* Sets coordinates for the facet whose effect is being calculated */ double x_a, x_b, x_c, y_a, y_b, y_c, z_a, z_b, z_c; struct GMTGRAVMAG3D_XYZ *triang = Ctrl->triang; @@ -1286,7 +1291,7 @@ GMT_LOCAL int facet_triangulate (struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VER } /* -----------------------------------------------------------------*/ -GMT_LOCAL int facet_raw (struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body_verts, unsigned int i, bool geo) { +GMT_LOCAL int facet_raw(struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body_verts, unsigned int i, bool geo) { /* Sets coordinates for the facet in the RAW format */ double cos_a, cos_b, cos_c, x_a, x_b, x_c, y_a, y_b, y_c, z_a, z_b, z_c; @@ -1310,7 +1315,7 @@ GMT_LOCAL int facet_raw (struct GMTGRAVMAG3D_CTRL *Ctrl, struct BODY_VERTS *body } /* ---------------------------------------------------------------------- */ -GMT_LOCAL void set_center (struct GMTGRAVMAG3D_CTRL *Ctrl) { +GMT_LOCAL void set_center(struct GMTGRAVMAG3D_CTRL *Ctrl) { /* Calculates triangle center by an approximate (iterative) formula */ int i, j, k = 5; double x, y, z, xa[6], ya[6], xb[6], yb[6], xc[6], yc[6]; @@ -1367,7 +1372,7 @@ GMT_LOCAL void gmtgravmag3d_triang_norm (int n_triang) { } #endif -GMT_LOCAL int check_triang_cw (struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, unsigned int type) { +GMT_LOCAL int check_triang_cw(struct GMTGRAVMAG3D_CTRL *Ctrl, unsigned int n, unsigned int type) { /* Checks that triangles are given in the correct clock-wise order. If not swap them. This is a tricky issue. In the case of "classic" trihedron (x positive right; y positive "north" and z positive up),