/*-------------------------------------------------------------------- * $Id: psxyz.c,v 1.2.4.3 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING 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 General Public License as published by * the Free Software Foundation; version 2 of the License. * * 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 General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * psxyz will read triplets from GMT_stdin and plot * the positions in 3-D using symbols selected by the user. If no * symbol size is specified, psxyz expects as input, where * size is the symbol size. Several map projections are supported. * For linear projections a 3-D basemap is provided. All symbols are * projected onto the xy plane (so that circles becomes ellipses) except * COLUMN and CUBE which are fully 3-D. * PostScript code is then written to GMT_stdout. * * Author: Paul Wessel * Date: 22-AUG-2000 * Version: 3.4.1 * */ #include "gmt.h" #define NONE 99 #define UNSPECIFIED 999 #define LINE 0 #define BAR 1 #define CROSS 2 #define POINT 3 #define CIRCLE 4 #define SQUARE 5 #define TRIANGLE 6 #define DIAMOND 7 #define STAR 8 #define HEXAGON 9 #define ITRIANGLE 10 #define CUBE 11 #define COLUMN 12 #define ELLIPSE 13 #define VECTOR 14 #define VECTOR2 15 #define PIE 16 #define TEXT 17 #define FRONT -18 /* FRONT symbols */ #define F_FAULT 0 #define F_TRIANGLE 1 #define F_SLIP 2 #define F_CIRCLE 3 #define F_BOX 4 /* Direction of front symbols: */ #define F_LEFT +1 #define F_CENTERED 0 #define F_RIGHT -1 #define POINTSIZE 0.005 double *xp, *yp; struct FRONTLINE { /* A sub-symbol for symbols along a front */ double f_gap; /* Gap between front symbols in inches */ double f_len; /* Length of front symbols in inches */ double f_off; /* Offset of first symbol from start of front in inches */ int f_sense; /* Draw symbols to left (+1), centered (0), or right (-1) of line */ int f_symbol; /* Which symbol to draw along the front line */ }; struct PSXYZ_SYMBOL { int symbol; /* Symbol id */ int n_required; /* Number of additional columns necessary to decode chosen symbol */ int u; /* Measure unit id (0 = cm, 1 = inch, 2 = m, 3 = point */ double size_x; /* Current symbol size in x */ double size_y; /* Current symbol size in y */ double given_size_x; /* Symbol size read from file or command line */ double given_size_y; /* Symbol size read from file or command line */ BOOLEAN read_size; /* TRUE when we must read symbol size from file */ int font_no; /* font to use for the -Sl symbol */ /* These apply to bar symbols */ double base; /* From what level to draw the bar */ BOOLEAN user_unit; /* if TRUE */ /* These apply to vectors */ BOOLEAN convert_angles; /* If TRUE, convert azimuth to angle on map */ BOOLEAN read_vector; /* if TRUE must read vector attributes */ BOOLEAN shrink; /* If TRUE, shrink vector attributes for small lengths */ double v_norm; /* shrink when lengths are smaller than this */ double v_shrink; /* Required scale factor */ double v_width; /* Width of vector stem in inches */ double h_length; /* Length of vector head in inches */ double h_width; /* Width of vector head in inches */ char string[64]; /* Character code to plot (could be octal) */ struct FRONTLINE f; /* parameters needed for a front */ } S; struct DATA1 { double lon, lat; double x, y, z, value, dist; float x_size, y_size, z_size; int symbol, flag; char *string; } *data1; struct DATA2 { double x, y, z; } *data2; void column3D(double x, double y, double z, double base, double x_size, double y_size, int *rgb, int outline); void cube3D(double x, double y, double z, double x_size, double y_size, int *rgb, int outline); void cross3D(double x, double y, double z, double size); void bar3D(double x, double y, double z, double base, double size, int *rgb, int outline); void square3D(double x, double y, double z, double size, int *rgb, int outline); void circle3D(double x, double y, double z, double size, int *rgb, int outline); void triangle3D(double x, double y, double z, double size, int *rgb, int outline); void itriangle3D(double x, double y, double z, double size, int *rgb, int outline); void hexagon3D(double x, double y, double z, double size, int *rgb, int outline); void star3D(double x, double y, double z, double size, int *rgb, int outline); void diamond3D(double x, double y, double z, double size, int *rgb, int outline); void ellipse3D (double x, double y, double z, double direction, double major, double minor, int rgb[], int outline); void pie3D (double x, double y, double z, double size, double dir1, double dir2, int rgb[], int outline); void sort_on_distance(struct DATA1 *data, int n); int decode_symbol_option (char *text, struct PSXYZ_SYMBOL *p, BOOLEAN cmd); void GMT_plot_ellipse (double lon, double lat, double z, double major, double minor, double azimuth, struct GMT_FILL fill, int outline); void GMT_draw_fence (double x[], double y[], int n, struct FRONTLINE *f, struct GMT_FILL *g, BOOLEAN outline); main (int argc, char **argv) { int i, j, n, n_files = 0, fno, rgb[9], n_cols_start = 3, n_fields; int n_alloc = GMT_CHUNK, n_args, three, four, five, n_required = 3; int n_total_read = 0, font_size, n_expected = 0; BOOLEAN error = FALSE, nofile = TRUE, polygon = FALSE, done, close_polygon = FALSE; BOOLEAN fill_set = FALSE, outline = FALSE, get_rgb = FALSE, read_symbol, set_z = FALSE, clip_set = FALSE; BOOLEAN skip_if_outside = TRUE, default_outline = FALSE, default_polygon = FALSE, cpt_given = FALSE; double west = 0.0, east = 0.0, south = 0.0, north = 0.0, new_z_level = 0.0; double lux[3], tmp, *in, x_1, x_2, y_1, y_2, z, xxx, yyy; double x2, y2, v_w, h_l, h_w; char line[BUFSIZ], *cpt, *p, *symb_arg = CNULL; FILE *fp = NULL; struct GMT_PEN pen, default_pen; struct GMT_FILL fill, default_fill; argc = GMT_begin (argc, argv); GMT_3D_mode = 1; /* Only do background axis first; do foreground at end */ GMT_init_pen (&pen, GMT_PENWIDTH); GMT_init_fill (&fill, -1, -1, -1); /* Default is no fill */ for (i = 0; i < 9; i++) rgb[i] = -1; /* Initialize PSXYZ_SYMBOL structure */ S.symbol = S.n_required = 0; S.size_x = S.size_y = S.v_norm = S.v_shrink = 0.0; S.read_size = S.user_unit = S.convert_angles = S.shrink = S.read_vector = FALSE; S.v_width = 0.03; memset ((void *)S.string, 0, (size_t)(16 * sizeof (char))); S.h_length = 0.12; S.h_width = 0.1; S.base = GMT_d_NaN; if ((S.u = gmtdefs.measure_unit) == GMT_CM) { S.v_width = 0.075 / 2.54; S.h_length = 0.3 / 2.54; S.h_width = 0.25 / 2.54; } /* Check and interpret the command line arguments */ for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch(argv[i][1]) { /* Common parameters */ case 'B': case 'H': case 'J': case 'K': case 'O': case 'P': case 'R': case 'U': case 'V': case 'X': case 'x': case 'Y': case 'y': case 'c': case ':': case '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; /* Supplemental parameters */ case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'C': cpt = &argv[i][2]; cpt_given = TRUE; break; case 'E': sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation); break; case 'G': /* Set color for symbol or polygon */ if (!argv[i][2] || (argv[i][2] && GMT_getfill (&argv[i][2], &fill))) { GMT_fill_syntax ('G'); error++; } fill_set = TRUE; break; case 'L': /* Force closed polygons */ close_polygon = TRUE; break; case 'M': /* Multiple line segments */ GMT_multisegment (&argv[i][2]); break; case 'N': /* Do not clip to map */ skip_if_outside = FALSE; break; case 'S': /* Get symbol [and size] */ symb_arg = &argv[i][2]; error += decode_symbol_option (&argv[i][2], &S, TRUE); break; case 'W': /* Set line attributes */ if (argv[i][2] && GMT_getpen (&argv[i][2], &pen)) { GMT_pen_syntax ('W'); error++; } outline = TRUE; break; case 'Z': if (argv[i][2]) { new_z_level = atof (&argv[i][2]); set_z = TRUE; } break; default: /* Options not recognized */ error = TRUE; GMT_default_error (argv[i][1]); break; } } else n_files++; } if (S.symbol > 0 && cpt_given) get_rgb = TRUE; if (argc == 1 || GMT_quick) { /* Display usage */ fprintf (stderr,"psxyz %s - Plot lines, polygons, and symbols in 3-D\n\n", GMT_VERSION); fprintf(stderr,"usage: psxyz -R -J\n"); fprintf(stderr, "\t-Jz [-B] [-C] [-E] [-G] [-H[]] [-K] [-L]\n"); fprintf(stderr, "\t[-M[]] [-N] [-O] [-P] [-S[/size_y]] [-U] [-V] [-W[]] [-X]\n"); fprintf(stderr, "\t[-Y] [-:] [-c] [-bi[s][]]\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\t is one or more files. If none, read standard input\n"); GMT_explain_option ('j'); GMT_explain_option ('Z'); GMT_explain_option ('r'); fprintf (stderr, "\n\tOPTIONS:\n"); GMT_explain_option ('b'); fprintf (stderr, "\t-C Use cpt-file to assign symbol colors based on t-value in 4th column (requires -S)\n"); fprintf (stderr, "\t Without -S, psxyz expects -M polygons and looks for -Z strings in each multiheader\n"); fprintf (stderr, "\t-E set user viewpoint azimuth and elevation [180/90].\n"); fprintf (stderr, "\t-G Specify color (for symbols/polygons) or pattern (for polygons). fill can be either\n"); fprintf (stderr, "\t 1) (each 0-255) for color or (0-255) for gray-shade [0].\n"); fprintf (stderr, "\t 2) p[or P]/ for predefined patterns (1-90).\n"); fprintf (stderr, "\t If -G is specified but not -S, psxyz draws a filled polygon.\n"); fprintf (stderr, "\t Default is no fill (transparent symbols or polygons)\n"); GMT_explain_option ('H'); GMT_explain_option ('K'); fprintf (stderr, "\t-L Force closed polygons\n"); GMT_explain_option ('M'); fprintf (stderr, "\t-N Do Not skip symbols that fall outside map border [Default will ignore those outside]\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf (stderr, "\t-S to select symbol type and symbol size (in %s). Choose between\n", GMT_unit_names[gmtdefs.measure_unit]); fprintf (stderr, "\t st(a)r, (b)ar, (c)ircle, (d)iamond, (e)llipse, (f)ront, (h)exagon, (i)nvtriangle\n"); fprintf (stderr, "\t (l)etter, c(o)lumn, (p)oint, (s)quare, (t)riangle, c(u)be, (v)ector, (w)edge, (x)cross\n"); fprintf (stderr, "\t If no size is specified, psxyz expects the 4th column to have sizes.\n"); fprintf (stderr, "\t If no symbol is specified, psxyz expects the last column to have symbol code.\n"); fprintf (stderr, "\t [Note: if -C is selected then 4th means 5th column, etc.]\n"); fprintf (stderr, "\t column and cube are true 3-D objects (give size as xsize/ysize), the rest is 2-D perspective only.\n"); fprintf (stderr, "\t Symbols A, C, D, H, I, S, T are adjusted to have same area of a circle of given size\n"); fprintf (stderr, "\t Bar (or Column): Append b to give the y- (or z-) value of the base [Default = 0 (1 for log-scales)]\n"); fprintf (stderr, "\t Ellipses: the direction, major, and minor axis must be in input columns 4, 5 and 6.\n"); fprintf (stderr, "\t If -SE rather than -Se is selected, psxyz will expect azimuth, and axes in km\n"); fprintf (stderr, "\t and convert azimuths based on the chosen map projection\n"); fprintf (stderr, "\t Fronts: Give tickgap/ticklen[dir][type][:offset], where\n"); fprintf (stderr, "\t dir = Plot symbol to the l(eft) or r(ight) of front [Default is centered]\n"); fprintf (stderr, "\t type = b(ox), c(ircle), f(ault), s(lip), t(riangle) [Default is f]\n"); fprintf (stderr, "\t box : square when centered, half-square otherwise.\n"); fprintf (stderr, "\t circle : full when centered, half-circle otherwise.\n"); fprintf (stderr, "\t fault : centered cross-tick or tick only in the direction.\n"); fprintf (stderr, "\t slip : left-lateral or right-lateral strike-slip arrows (centered is not defined).\n"); fprintf (stderr, "\t triangle : diagonal square when centered, directed triangle otherwise.\n"); fprintf (stderr, "\t offset = start plotting symbols when along-track distance equals offset [0].\n"); fprintf (stderr, "\t Letter: append / after symbol size, and optionally %%\n"); fprintf (stderr, "\t Vectors: the direction and length must be in input columns 4 and 5.\n"); fprintf (stderr, "\t Furthermore, means arrowwidth/headlength/headwith [Default is %lg/%lg/%lg]\n", S.v_width, S.h_length, S.h_width); fprintf (stderr, "\t If -SV rather than -Sv is selected, psxyz will expect azimuth and length\n"); fprintf (stderr, "\t and convert azimuths based on the chosen map projection\n"); GMT_explain_option ('U'); fprintf (stderr, "\t Wedges: the start and stop directions of pie wedge must be in input columns 4 and 5.\n"); GMT_explain_option ('V'); fprintf (stderr, "\t-W sets pen attributes [width = %lgp, color = (%d/%d/%d), solid line].\n", pen.width, pen.rgb[0], pen.rgb[1], pen.rgb[2]); fprintf (stderr, "\t Implicitly draws symbol outline with this pen.\n"); GMT_explain_option ('X'); GMT_explain_option ('c'); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf (stderr, "\t Default is the required number of columns\n"); GMT_explain_option ('.'); exit (EXIT_FAILURE); } /* Check that the options selected are mutually consistant */ if (!project_info.region_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program); error++; } if (fill.use_pattern && S.symbol != LINE) { /* fill-pattern only for polygons */ error++; fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Fill-pattern only used with polygons\n", GMT_program); } if (cpt_given && S.symbol == 0 && !GMT_io.multi_segments) { error++; fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Must also specify a symbol (see -S) or give polygon file with -M\n", GMT_program); } if (GMT_io.binary[0] && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } if (GMT_io.binary[0] && S.symbol == UNSPECIFIED) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have symbol information\n", GMT_program); error++; } if (error) exit (EXIT_FAILURE); read_symbol = (S.symbol == UNSPECIFIED); if (S.symbol == 0 && (fill_set || close_polygon)) polygon = TRUE; if (S.symbol == FRONT) polygon = FALSE; default_pen = pen; default_fill = fill; default_outline = outline; default_polygon = polygon; if (get_rgb) n_cols_start++; three = (get_rgb) ? 4 : 3; four = (get_rgb) ? 5 : 4; five = (get_rgb) ? 6 : 5; n_required = n_cols_start + S.n_required; n_expected = (GMT_io.binary[0]) ? GMT_io.ncol[0] : n_required; if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0]= n_required; if (n_expected < n_required) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data must have at least %d columns\n", GMT_program, n_required); error++; } GMT_put_history (argc, argv); /* Update .gmtcommands */ if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } if (cpt_given) GMT_read_cpt (cpt); if (n_files > 0) nofile = FALSE; else n_files = 1; n_args = (argc > 1) ? argc : 2; GMT_map_setup (west, east, south, north); if (fill.rgb[0] >= 0 && (S.symbol == COLUMN || S.symbol == CUBE)) { /* Modify the color for each facet */ lux[0] = fabs (z_project.sin_az * z_project.cos_el); lux[1] = fabs (z_project.cos_az * z_project.cos_el); lux[2] = fabs (z_project.sin_el); tmp = MAX (lux[0], MAX (lux[1], lux[2])); for (i = 0; i < 3; i++) { lux[i] = (lux[i] / tmp) - 0.5; rgb[3*i] = fill.rgb[0]; rgb[3*i+1] = fill.rgb[1]; rgb[3*i+2] = fill.rgb[2]; GMT_illuminate (lux[i], &rgb[3*i]); } } ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin, gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0])); GMT_echo_command (argc, argv); if (gmtdefs.unix_time) GMT_timestamp (argc, argv); ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0); if (set_z) project_info.z_level = new_z_level; if (frame_info.plot) GMT_map_basemap (); if (project_info.z_pars[0] == 0.0 && skip_if_outside) { GMT_map_clip_on (GMT_no_rgb, 3); clip_set = TRUE; } if (S.symbol == ELLIPSE) skip_if_outside = FALSE; GMT_setpen (&pen); if (S.symbol == TEXT && fill_set && !outline) ps_setpaint (fill.rgb); if (S.symbol == TEXT) { font_size = irint (S.size_x * 72.0); /* To get points */ ps_setfont (S.font_no); /* Set the required font */ } done = FALSE; for (fno = 1; !done && fno < n_args; fno++) { /* Loop over all input files */ if (!nofile && argv[fno][0] == '-') continue; if (nofile) { fp = GMT_stdin; done = TRUE; #ifdef SET_IO_MODE GMT_setmode (0); #endif } else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]); continue; } if (!nofile && gmtdefs.verbose) { fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]); sprintf (line, "File: %s", argv[fno]); ps_comment (line); } if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp); if (S.symbol > 0) { /* symbol part */ GMT_world_map = (S.symbol == ELLIPSE && S.convert_angles) ? FALSE : TRUE; if (read_symbol) n_expected = BUFSIZ; data1 = (struct DATA1 *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (struct DATA1), GMT_program); n = 0; while ((n_fields = GMT_input (fp, &n_expected, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) { while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) { /* Skip segment headers */ if (gmtdefs.verbose) ps_comment (GMT_io.segment_header); if ((p = strstr (GMT_io.segment_header, "-G"))) { strcpy (line, &p[2]); for (i = 0; line[i]; i++) if (line[i] == ' ') line[i] = '\0'; if (GMT_getfill (line, &fill)) { fprintf (stderr, "%s: Trouble decoding fill info from header near line %d\n", GMT_program, n_total_read); } } else if (cpt_given && (p = strstr (GMT_io.segment_header, "-Z"))) { /* Set symbol r/g/b via cpt-lookup */ z = (p[2] == 'N') ? GMT_d_NaN : atof (&p[2]); GMT_get_rgb24 (z, fill.rgb); } else fill = default_fill; if ((p = strstr (GMT_io.segment_header, "-W"))) { strcpy (line, &p[2]); for (i = 0; line[i]; i++) if (line[i] == ' ') line[i] = '\0'; if (GMT_getpen (line, &pen)) { fprintf (stderr, "%s: Trouble decoding pen info from header near line %d\n", GMT_program, n_total_read); } else { GMT_setpen (&pen); outline = TRUE; } } else { GMT_setpen (&default_pen); outline = default_outline; } if (read_symbol) n_expected = BUFSIZ; n_fields = GMT_input (fp, &n_expected, &in); } if ((GMT_io.status & GMT_IO_EOF)) continue; /* At EOF */ if (GMT_io.status & GMT_IO_MISMATCH && !read_symbol) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected, n_total_read); exit (EXIT_FAILURE); } n_total_read++; S.size_x = S.size_y = 0.0; if (S.read_size) { S.size_x = in[three] * GMT_u2u[S.u][GMT_INCH]; S.size_y = (S.symbol == COLUMN) ? in[four] * GMT_u2u[S.u][GMT_INCH] : S.size_x; } if (read_symbol) { i = strlen (GMT_io.current_record) - 1; while (GMT_io.current_record[i] && !strchr (" ,\t", (int)GMT_io.current_record[i])) i--; decode_symbol_option (&GMT_io.current_record[i+1], &S, FALSE); n_required = n_cols_start + S.n_required; } else { decode_symbol_option (symb_arg, &S, FALSE); n_required = n_cols_start + S.n_required; } data1[n].flag = 0; data1[n].lon = in[0]; data1[n].lat = in[1]; if (S.read_vector) { data1[n].x_size = (float)in[three]; /* direction */ data1[n].y_size = (float)(in[four] * GMT_u2u[S.u][GMT_INCH]); /* length */ data1[n].z_size = (float)S.v_shrink; /* Normalization shrink (if given) */ if (S.convert_angles) data1[n].flag = 1; /* Azim vs Dir. */ if (data1[i].y_size < S.v_norm) data1[n].flag |= 2; } else if (S.symbol == ELLIPSE) { data1[n].x_size = (float)in[three]; /* direction */ data1[n].y_size = (float)in[four]; /* major */ data1[n].z_size = (float)in[five]; /* minor */ if (S.convert_angles) { /* Got axis in km */ data1[n].flag = 1; /* Azim vs Dir. */ } else { /* Got axes in current units, change to inches */ data1[n].y_size *= (float)GMT_u2u[S.u][GMT_INCH]; data1[n].z_size *= (float)GMT_u2u[S.u][GMT_INCH]; } } else if (S.symbol == PIE) { data1[n].x_size = (float)S.size_x; data1[n].y_size = (float)(D2R * in[three+S.read_size]); /* Stop direction in radians */ data1[n].z_size = (float)(D2R * in[four+S.read_size]); /* Stop direction in radians */ } else { data1[n].x_size = (float)S.size_x; data1[n].y_size = (float)S.size_y; } if (S.user_unit) data1[n].flag |= 4; /* Skip zero-size symbols */ if (S.symbol != POINT && S.symbol < ELLIPSE && S.size_x <= 0.0) { if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Symbol size <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read); continue; } if (S.read_vector && data1[n].y_size <= 0.0) { if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: vector length <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read); continue; } if (S.symbol == ELLIPSE && (data1[n].y_size <= 0.0 || data1[n].z_size <= 0.0)) { if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: ellipse axes <= 0.0 near line %d (skipped)\n", GMT_program, n_total_read); continue; } if (skip_if_outside) { if (in[2] < project_info.z_bottom || in[2] > project_info.z_top) continue; GMT_map_outside (in[0], in[1]); if ( abs (GMT_x_status_new) > 1 || abs (GMT_y_status_new) > 1) continue; } GMT_project3D (in[0], in[1], in[2], &data1[n].x, &data1[n].y, &data1[n].z); if (get_rgb) data1[n].value = in[3]; if (S.symbol == COLUMN) { GMT_z_to_zz (S.base, &tmp); data1[n].z_size = (float)tmp; } if (S.symbol == BAR) { GMT_y_to_yy (S.base, &tmp); data1[n].z_size = (float)tmp; } if (S.symbol == TEXT) { data1[n].string = (char *) GMT_memory (VNULL, (size_t) (strlen (S.string) + 1), sizeof (char), GMT_program); strcpy (data1[n].string, S.string); } data1[n].symbol = S.symbol; n++; if (n == n_alloc) { n_alloc += GMT_CHUNK; data1 = (struct DATA1 *) GMT_memory ((void *)data1, (size_t)n_alloc, sizeof (struct DATA1), GMT_program); } if (read_symbol) n_expected = BUFSIZ; } data1 = (struct DATA1 *) GMT_memory ((void *)data1, (size_t)n, sizeof (struct DATA1), GMT_program); /* Sort according to distance from viewer */ sort_on_distance (data1, n); /* Now plot these symbols one at the time */ for (i = 0; i < n; i++) { if (get_rgb) { GMT_get_rgb24 (data1[i].value, fill.rgb); if (data1[i].symbol == COLUMN || data1[i].symbol == CUBE) { for (j = 0; j < 3; j++) { rgb[3*j] = fill.rgb[0]; rgb[3*j+1] = fill.rgb[1]; rgb[3*j+2] = fill.rgb[2]; GMT_illuminate (lux[j], &rgb[3*j]); } } } if (GMT_cpt_skip) continue; /* Skip this because cpt says so */ switch (data1[i].symbol) { case NONE: break; case STAR: star3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case BAR: if (data1[i].flag & 4) { GMT_geo_to_xy (data1[i].x-data1[i].x_size, data1[i].y, &x_1, &y_1); GMT_geo_to_xy (data1[i].x+data1[i].x_size, data1[i].y, &x_2, &y_2); xxx = 0.5 * hypot (x_1 - x_2, y_1 - y_2); bar3D (data1[i].x, data1[i].y, data1[i].z, data1[i].z_size, xxx, fill.rgb, outline); } else bar3D (data1[i].x, data1[i].y, data1[i].z, data1[i].z_size, (double)data1[i].x_size, fill.rgb, outline); break; case COLUMN: if (data1[i].flag & 4) { GMT_geo_to_xy (data1[i].x-data1[i].x_size, data1[i].y, &x_1, &y_1); GMT_geo_to_xy (data1[i].x+data1[i].x_size, data1[i].y, &x_2, &y_2); xxx = 0.5 * hypot (x_1 - x_2, y_1 - y_2); GMT_geo_to_xy (data1[i].x, data1[i].y-data1[i].y_size, &x_1, &y_1); GMT_geo_to_xy (data1[i].x, data1[i].y+data1[i].y_size, &x_2, &y_2); yyy = 0.5 * hypot (x_1 - x_2, y_1 - y_2); column3D (data1[i].x, data1[i].y, data1[i].z, data1[i].z_size, xxx, yyy, rgb, outline); } else column3D (data1[i].x, data1[i].y, data1[i].z, data1[i].z_size, (double)data1[i].x_size, (double)data1[i].y_size, rgb, outline); break; case CROSS: cross3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size); break; case POINT: cross3D (data1[i].x, data1[i].y, data1[i].z, POINTSIZE); break; case CIRCLE: circle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case SQUARE: square3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case HEXAGON: hexagon3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case TRIANGLE: triangle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case ITRIANGLE: itriangle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case DIAMOND: diamond3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.rgb, outline); break; case CUBE: cube3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, (double)data1[i].y_size, rgb, outline); break; case ELLIPSE: if (data1[i].flag) GMT_plot_ellipse (data1[i].lon, data1[i].lat, data1[i].z, data1[i].y_size, data1[i].z_size, data1[i].x_size, fill, outline); else ellipse3D (data1[i].x, data1[i].y, data1[i].z, data1[i].x_size, data1[i].y_size, data1[i].z_size, fill.rgb, outline); break; case TEXT: font_size = irint (2.0 * data1[i].x_size * 72.0); /* To get points */ if (outline && fill_set) { ps_setpaint (fill.rgb); GMT_text3d (data1[i].x, data1[i].y, data1[i].z, font_size, S.font_no, data1[i].string, 0.0, 6, FALSE); ps_setpaint (pen.rgb); GMT_text3d (data1[i].x, data1[i].y, data1[i].z, font_size, S.font_no, data1[i].string, 0.0, 6, TRUE); } else if (fill_set) GMT_text3d (data1[i].x, data1[i].y, data1[i].z, font_size, S.font_no, data1[i].string, 0.0, 6, FALSE); else GMT_text3d (data1[i].x, data1[i].y, data1[i].z, font_size, S.font_no, data1[i].string, 0.0, 6, TRUE); break; case VECTOR: if (data1[i].y_size <= 0.0) continue; if (data1[i].flag & 1) { GMT_azim_to_angle (data1[i].lon, data1[i].lat, 0.1, data1[i].x_size, &tmp); data1[i].x_size = (float)tmp; } x2 = data1[i].x + data1[i].y_size * cosd (data1[i].x_size); y2 = data1[i].y + data1[i].y_size * sind (data1[i].x_size); GMT_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, S.v_width, S.h_length, S.h_width, gmtdefs.vector_shape, fill.rgb, outline); break; case VECTOR2: if (data1[i].y_size <= 0.0) continue; if (data1[i].flag & 1) { GMT_azim_to_angle (data1[i].lon, data1[i].lat, 0.1, data1[i].x_size, &tmp); data1[i].x_size = (float)tmp; } x2 = data1[i].x + data1[i].y_size * cosd (data1[i].x_size); y2 = data1[i].y + data1[i].y_size * sind (data1[i].x_size); if (data1[i].flag & 2) { /* Scale arrow attributes down with length */ v_w = S.v_width * data1[i].y_size * data1[i].z_size; h_l = S.h_length * data1[i].y_size * data1[i].z_size; h_w = S.h_width * data1[i].y_size * data1[i].z_size; GMT_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, v_w, h_l, h_w, gmtdefs.vector_shape, fill.rgb, outline); } else /* Leave as specified */ GMT_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, S.v_width, S.h_length, S.h_width, gmtdefs.vector_shape, fill.rgb, outline); break; case PIE: pie3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, (double)data1[i].y_size, (double)data1[i].z_size, fill.rgb, outline); break; } } GMT_free ((void *)data1); } else { /* Line/polygon part */ data2 = (struct DATA2 *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (struct DATA2), GMT_program); n_required = 3; n_fields = GMT_input (fp, &n_expected, &in); while (! (GMT_io.status & GMT_IO_EOF)) { while (GMT_io.status & GMT_IO_SEGMENT_HEADER) { if (gmtdefs.verbose) ps_comment (GMT_io.segment_header); if ((p = strstr (GMT_io.segment_header, "-G"))) { strcpy (line, &p[2]); for (i = 0; line[i]; i++) if (line[i] == ' ') line[i] = '\0'; if (GMT_getfill (line, &fill)) { fprintf (stderr, "%s: Trouble decoding fill info from header near line %d\n", GMT_program, n_total_read); } polygon = TRUE; } else if (cpt_given && (p = strstr (GMT_io.segment_header, "-Z"))) { /* Set polygon r/g/b via cpt-lookup */ z = (p[2] == 'N') ? GMT_d_NaN : atof (&p[2]); GMT_get_rgb24 (z, fill.rgb); fill.use_pattern = FALSE; polygon = TRUE; } else { polygon = default_polygon; fill = default_fill; } if ((p = strstr (GMT_io.segment_header, "-W"))) { strcpy (line, &p[2]); for (i = 0; line[i]; i++) if (line[i] == ' ') line[i] = '\0'; if (GMT_getpen (line, &pen)) { fprintf (stderr, "%s: Trouble decoding pen info from header near line %d\n", GMT_program, n_total_read); } else { GMT_setpen (&pen); outline = TRUE; } } else { GMT_setpen (&default_pen); outline = default_outline; } n_fields = GMT_input (fp, &n_expected, &in); n_total_read++; } if ((GMT_io.status & GMT_IO_EOF)) continue; /* At EOF */ n = 0; while (! (GMT_io.status & (GMT_IO_SEGMENT_HEADER | GMT_IO_EOF))) { /* Keep going until FALSE or = 2 segment header */ if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected, n_total_read); exit (EXIT_FAILURE); } data2[n].x = in[0]; data2[n].y = in[1]; data2[n].z = in[2]; n++; if (n == n_alloc) { n_alloc += GMT_CHUNK; data2 = (struct DATA2 *) GMT_memory ((void *)data2, (size_t)n_alloc, sizeof (struct DATA2), GMT_program); } n_fields = GMT_input (fp, &n_expected, &in); } if (polygon) { /* Explicitly close polygon */ data2[n].x = data2[0].x; data2[n].y = data2[0].y; data2[n].z = data2[0].z; n++; } data2 = (struct DATA2 *) GMT_memory ((void *)data2, (size_t)n, sizeof (struct DATA2), GMT_program); n_alloc = n; if (!GMT_cpt_skip) { xp = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program); yp = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program); for (i = 0; i < n; i++) GMT_geoz_to_xy (data2[i].x, data2[i].y, data2[i].z, &xp[i], &yp[i]); if (polygon) GMT_fill (xp, yp, n, &fill, outline); else ps_line (xp, yp, n, 3, FALSE, TRUE); if (S.symbol == FRONT) { for (i = 0; i < n; i++) GMT_geo_to_xy (data2[i].x, data2[i].y, &xp[i], &yp[i]); GMT_draw_fence (xp, yp, n, &S.f, &fill, outline); /* Must draw fault crossbars */ } GMT_free ((void *)xp); GMT_free ((void *)yp); } n = 0; } GMT_free ((void *)data2); } if (fp != GMT_stdin) GMT_fclose (fp); } if (clip_set) GMT_map_clip_off(); if (pen.texture) ps_setdash (CNULL, 0); if (project_info.three_D) GMT_vertical_axis (2); /* Draw background axis */ ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0); ps_plotend (gmtdefs.last_page); GMT_end (argc, argv); } void column3D (double x, double y, double z, double base, double x_size, double y_size, int *rgb, int outline) { int i, j, k; double zz, xp[4], yp[4], zp[4], plot_x[4], plot_y[4], top, sign; x_size *= 0.5; y_size *= 0.5; top = z; if (top < base) d_swap (top, base); for (i = 0; i < 3; i++) { sign = -1.0; zz = base; switch (z_project.face[i]) { case 0: /* yz plane positive side */ sign = 1.0; case 1: /* negative side */ xp[0] = xp[1] = xp[2] = xp[3] = x + sign * x_size; yp[0] = yp[3] = y - y_size; yp[1] = yp[2] = y + y_size; zp[0] = zp[1] = base; zp[2] = zp[3] = top; break; case 2: /* xz plane positive side */ sign = 1.0; case 3: /* negative side */ xp[0] = xp[3] = x - x_size; xp[1] = xp[2] = x + x_size; yp[0] = yp[1] = yp[2] = yp[3] = y + sign * y_size; zp[0] = zp[1] = base; zp[2] = zp[3] = top; break; case 4: /* xy plane positive side */ zz = top; case 5: /* negative side */ xp[0] = xp[3] = x - x_size; yp[0] = yp[1] = y - y_size; xp[1] = xp[2] = x + x_size; yp[2] = yp[3] = y + y_size; zp[0] = zp[1] = zp[2] = zp[3] = zz; break; } k = z_project.face[i] / 2; for (j = 0; j < 4; j++) GMT_xyz_to_xy (xp[j], yp[j], zp[j], &plot_x[j], &plot_y[j]); ps_patch (plot_x, plot_y, 4, &rgb[3*k], outline); } } void cube3D (double x, double y, double z, double x_size, double y_size, int *rgb, int outline) { int i, j, k; double xp[4], yp[4], zp[4], plot_x[4], plot_y[4], sign; x_size *= 0.5; y_size *= 0.5; for (i = 0; i < 3; i++) { sign = -1.0; switch (z_project.face[i]) { case 4: /* xy plane positive side */ sign = 1.0; case 5: /* negative side */ xp[0] = xp[3] = x - x_size; yp[0] = yp[1] = y - y_size; xp[1] = xp[2] = x + x_size; yp[2] = yp[3] = y + y_size; zp[0] = zp[1] = zp[2] = zp[3] = z + sign * x_size; break; case 2: /* xz plane positive side */ sign = 1.0; case 3: /* negative side */ xp[0] = xp[3] = x - x_size; xp[1] = xp[2] = x + x_size; yp[0] = yp[1] = yp[2] = yp[3] = y + sign * y_size; zp[0] = zp[1] = z - x_size; zp[2] = zp[3] = z + x_size; break; case 0: /* yz plane positive side */ sign = 1.0; case 1: /* negative side */ xp[0] = xp[1] = xp[2] = xp[3] = x + sign * x_size; yp[0] = yp[3] = y - y_size; yp[1] = yp[2] = y + y_size; zp[0] = zp[1] = z - x_size; zp[2] = zp[3] = z + x_size; break; } k = z_project.face[i] / 2; for (j = 0; j < 4; j++) GMT_xyz_to_xy (xp[j], yp[j], zp[j], &plot_x[j], &plot_y[j]); ps_patch (plot_x, plot_y, 4, &rgb[3*k], outline); } } void cross3D (double x, double y, double z, double size) { double xp[2], yp[2], plot_x, plot_y; size *= 0.5; xp[0] = x - size; xp[1] = x + size; yp[0] = y - size; yp[1] = y + size; GMT_xyz_to_xy (xp[0], yp[0], z, &plot_x, &plot_y); ps_plot (plot_x, plot_y, 3); GMT_xyz_to_xy (xp[1], yp[1], z, &plot_x, &plot_y); ps_plot (plot_x, plot_y, 2); GMT_xyz_to_xy (xp[1], yp[0], z, &plot_x, &plot_y); ps_plot (plot_x, plot_y, 3); GMT_xyz_to_xy (xp[0], yp[1], z, &plot_x, &plot_y); ps_plot (plot_x, plot_y, 2); } void bar3D (double x, double y, double z, double base, double size, int rgb[], int outline) { int i; double xp[4], yp[4], plot_x[4], plot_y[4]; size *= 0.5; xp[0] = xp[3] = x - size; xp[1] = xp[2] = x + size; yp[0] = yp[1] = base; yp[2] = yp[3] = y; for (i = 0; i < 4; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 4, rgb, outline); } void square3D (double x, double y, double z, double size, int rgb[], int outline) { int i; double xp[4], yp[4], plot_x[4], plot_y[4]; size *= 0.3535533906; xp[0] = xp[3] = x - size; xp[1] = xp[2] = x + size; yp[0] = yp[1] = y - size; yp[2] = yp[3] = y + size; for (i = 0; i < 4; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 4, rgb, outline); } void circle3D (double x, double y, double z, double size, int rgb[], int outline) { /* Must plot a squashed circle */ int i; double xx, yy, a, da, s, c, plot_x[51], plot_y[51]; da = TWO_PI / 50.0; size *= 0.5; for (i = 0; i <= 50; i++) { a = i * da; sincos (a, &s, &c); xx = x + size * c; yy = y + size * s; GMT_xyz_to_xy (xx, yy, z, &plot_x[i], &plot_y[i]); } ps_polygon (plot_x, plot_y, 51, rgb, outline); } void ellipse3D (double x, double y, double z, double direction, double major, double minor, int rgb[], int outline) { /* Must plot a squashed circle */ int i; double dx, dy, a, da, s, c, sin_direction, cos_direction, x_prime, y_prime, plot_x[51], plot_y[51]; sincos (direction * D2R, &sin_direction, &cos_direction); da = TWO_PI / 50.0; for (i = 0; i <= 50; i++) { a = i * da; sincos (a, &s, &c); dx = major * c; dy = minor * s; x_prime = x + dx * cos_direction - dy * sin_direction; y_prime = y + dx * sin_direction + dy * cos_direction; GMT_xyz_to_xy (x_prime, y_prime, z, &plot_x[i], &plot_y[i]); } ps_polygon (plot_x, plot_y, 51, rgb, outline); } void pie3D (double x, double y, double z, double size, double dir1, double dir2, int rgb[], int outline) { /* Must plot a squashed pie wedge */ int i, j, n; double arc, xx, yy, a, da, s, c, plot_x[52], plot_y[52]; arc = (dir2 - dir1); while (arc > TWO_PI) arc -= TWO_PI; da = TWO_PI / 50.0; /* Standard step length for full circle */ n = irint (arc / da); /* But we are not doing a full circle so use less than 50 points */ da = arc / n; /* Step length for this pie wedge */ size *= 0.5; GMT_xyz_to_xy (x, y, z, &plot_x[0], &plot_y[0]); /* Start from center */ for (i = 0, j = 1; i <= n; i++, j++) { a = dir1 + i * da; sincos (a, &s, &c); xx = x + size * c; yy = y + size * s; GMT_xyz_to_xy (xx, yy, z, &plot_x[j], &plot_y[j]); } ps_polygon (plot_x, plot_y, j, rgb, outline); } void triangle3D (double x, double y, double z, double size, int rgb[], int outline) { int i; double xp[3], yp[3], plot_x[3], plot_y[3]; xp[0] = x - 0.433012701892*size; yp[0] = yp[1] = y - 0.25 * size; xp[1] = x + 0.433012701892*size; xp[2] = x; yp[2] = y + 0.5 * size; for (i = 0; i < 3; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 3, rgb, outline); } void itriangle3D (double x, double y, double z, double size, int rgb[], int outline) { int i; double xp[3], yp[3], plot_x[3], plot_y[3]; xp[0] = x - 0.433012701892*size; yp[0] = yp[1] = y + 0.25 * size; xp[1] = x + 0.433012701892*size; xp[2] = x; yp[2] = y - 0.5 * size; for (i = 0; i < 3; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 3, rgb, outline); } void diamond3D (double x, double y, double z, double size, int rgb[], int outline) { int i; double xp[4], yp[4], plot_x[4], plot_y[4]; size *= 0.5; xp[0] = xp[2] = x; xp[1] = x - size; xp[3] = x + size; yp[0] = y - size; yp[1] = yp[3] = y; yp[2] = y + size; for (i = 0; i < 4; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 4, rgb, outline); } void hexagon3D (double x, double y, double z, double size, int rgb[], int outline) { int i; double xp[6], yp[6], plot_x[6], plot_y[6], sx, sy; size *= 0.5; sx = 0.5 * size; sy = 0.8660254038 * size; xp[0] = x + size; yp[0] = y; xp[1] = x + sx; yp[1] = y + sy; xp[2] = x - sx; yp[2] = yp[1]; xp[3] = x - size; yp[3] = y; xp[4] = xp[2]; yp[4] = y - sy; xp[5] = xp[1]; yp[5] = yp[4]; for (i = 0; i < 6; i++) GMT_xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]); ps_patch (plot_x, plot_y, 6, rgb, outline); } void star3D (double x, double y, double z, double size, int rgb[], int outline) { int i, k; double xx, yy, plot_x[10], plot_y[10], a, s2; size *= 0.5; s2 = 0.38196601125 * size; for (i = k = 0; i < 5; i++) { a = -54.0 + i * 72.0; xx = x + size * cosd (a); yy = y + size * sind (a); GMT_xyz_to_xy (xx, yy, z, &plot_x[k], &plot_y[k]); k++; a += 36.0; xx = x + s2 * cosd (a); yy = y + s2 * sind (a); GMT_xyz_to_xy (xx, yy, z, &plot_x[k], &plot_y[k]); k++; } ps_patch (plot_x, plot_y, 10, rgb, outline); } void sort_on_distance (struct DATA1 *data, int n) { /* This function sorts the data array such that points farthest away are plotted first */ int i; int dist_compare(const void *point_1, const void *point_2); double dx, dy, x0, y0, x, y, dr, a, b, c; x0 = 0.5 * (project_info.xmin + project_info.xmax); y0 = 0.5 * (project_info.ymin + project_info.ymax); dx = 0.5 * (project_info.xmax - project_info.xmin); dy = 0.5 * (project_info.ymax - project_info.ymin); dr = hypot (dx, dy); x = x0 - dr * z_project.sin_az; y = y0 - dr * z_project.cos_az; if (z_project.cos_az == 0.0) { a = 1.0; b = 0.0; c = x; } else { a = -tan (z_project.view_azimuth); b = -1.0; c = y - x * a; } for (i = 0; i < n; i++) data[i].dist = fabs (a * data[i].x + b * data[i].y + c); qsort ((void *)data, (size_t)n, sizeof (struct DATA1), dist_compare); } int dist_compare (const void *point_1, const void *point_2) { int first; if ( ((struct DATA1 *)point_1)->dist > ((struct DATA1 *)point_2)->dist) return (-1); else if ( ((struct DATA1 *)point_1)->dist < ((struct DATA1 *)point_2)->dist) return (1); else { first = ( ((struct DATA1 *)point_1)->z < ((struct DATA1 *)point_2)->z); if (first && z_project.view_elevation >= 0.0) return (-1); else if (first && z_project.view_elevation < 0.0) return (1); else return (0); } } int decode_symbol_option (char *text, struct PSXYZ_SYMBOL *p, BOOLEAN cmd) { int decode_error = 0, bset = 0, j, n, k, len, slash = 0; BOOLEAN check = TRUE, old_style; char symbol_type, txt_a[32], txt_b[32], txt_c[32], text_cp[128], *c; p->n_required = 0; if (!text[0]) { /* No symbol or size given */ p->size_x = p->size_y = 0.0; symbol_type = '+'; } else if (isdigit ((int)text[0]) || text[0] == '.') { /* Size, but no symbol given */ n = sscanf (text, "%[^/]/%s", txt_a, txt_b); p->size_x = p->given_size_x = GMT_convert_units (txt_a, GMT_INCH); if (n == 2) p->size_y = p->given_size_y = GMT_convert_units (txt_b, GMT_INCH); else if (n == 1) p->size_y = p->given_size_y = p->size_x; else decode_error = TRUE; symbol_type = '+'; } else if (text[0] == 'l') { /* Letter symbol is special case */ strcpy (text_cp, text); if ((c = strchr (text_cp, '%'))) { /* Gave font name or number, too */ *c = ' '; /* Make the % a space */ c++; /* Go to next character */ if (c[0] >= '0' && c[0] <= '9') /* Gave a font # */ p->font_no = atoi (c); else p->font_no = GMT_key_lookup (c, GMT_font_name, N_FONTS); if (p->font_no >= N_FONTS) { fprintf (stderr, "%s: -Sl contains bad font (set to %s (0))\n", GMT_program, GMT_font_name[gmtdefs.anot_font]); p->font_no = gmtdefs.anot_font; } } if (text[1] == '/') { /* No size given */ symbol_type = 'l'; if (p->size_x == 0.0) p->size_x = p->given_size_x; if (p->size_y == 0.0) p->size_y = p->given_size_y; } else { n = sscanf (text_cp, "%c%[^/]/%s", &symbol_type, txt_a, p->string); p->size_x = p->given_size_x = GMT_convert_units (txt_a, GMT_INCH); decode_error = (n != 3); } } else if (strchr ("aAbCcDdeEfHhIipoSsTtuVvwx", (int) text[0]) && strchr ("CcIiMmPp", (int) text[1])) { /* Symbol, but no size given (size assumed given on command line), only unit information */ n = sscanf (text, "%c", &symbol_type); if (p->size_x == 0.0) p->size_x = p->given_size_x; if (p->size_y == 0.0) p->size_y = p->given_size_y; if (text[1] && (p->u = GMT_get_unit (text[1])) < 0) decode_error = TRUE; } else if (strchr ("aAbCcDdeEfHhIipoSsTtuVvwx", (int) text[0]) && (text[1] == '\n' || text[1] == '\0')) { /* Symbol, but no size given (size assumed given on command line) */ n = sscanf (text, "%c", &symbol_type); if (p->size_x == 0.0) p->size_x = p->given_size_x; if (p->size_y == 0.0) p->size_y = p->given_size_y; } else if (strchr ("bou", (int) text[0])) { /* Bar, column, cube with size */ /* Bar: -Sb[c|i|m|p|u][b] */ /* Column: -So[c|i|m|p][/[c|i|m|p]][u][b] */ /* Cube: -Su[c|i|m|p|u] */ for (j = 0; text[j]; j++) { if (text[j] == '/') slash = j; if (text[j] == 'b') bset = j; } strcpy (text_cp, text); if (bset) text_cp[bset] = 0; /* Chop off the b from copy */ if ((bset && text_cp[bset-1] == 'u') || (j && text[j-1] == 'u')) p->user_unit = TRUE; if (slash) { /* Separate x/y sizes */ n = sscanf (text_cp, "%c%[^/]/%s", &symbol_type, txt_a, txt_b); decode_error = (n != 3); if (p->user_unit) { p->size_x = p->given_size_x = atof (txt_a); p->size_y = p->given_size_y = atof (txt_b); } else { p->size_x = p->given_size_x = GMT_convert_units (txt_a, GMT_INCH); p->size_y = p->given_size_y = GMT_convert_units (txt_b, GMT_INCH); } } else { /* Only a single x = y size */ n = sscanf (text_cp, "%c%s", &symbol_type, txt_a); if (n == 2) { if (p->user_unit) { p->size_x = p->given_size_x = atof (txt_a); p->size_y = p->given_size_y = p->size_x; } else { p->size_x = p->given_size_x = GMT_convert_units (txt_a, GMT_INCH); p->size_y = p->given_size_y = p->size_x; } } else { if (p->size_x == 0.0) p->size_x = p->given_size_x; if (p->size_y == 0.0) p->size_y = p->given_size_y; } } } else { n = sscanf (text, "%c%[^/]/%s", &symbol_type, txt_a, txt_b); p->size_x = p->given_size_x = GMT_convert_units (txt_a, GMT_INCH); if (n == 3) p->size_y = p->given_size_y = GMT_convert_units (txt_b, GMT_INCH); else if (n == 2) p->size_y = p->given_size_y = p->size_x; else decode_error = TRUE; } check = TRUE; switch (symbol_type) { case '+': p->symbol = UNSPECIFIED; break; case 'A': p->size_x *= 1.6728932614; /* To equal area of circle with same size */ case 'a': p->symbol = STAR; break; case 'b': p->symbol = BAR; if (bset) p->base = atof (&text[bset+1]); break; case 'C': case 'c': p->symbol = CIRCLE; break; case 'D': p->size_x *= 1.2533141373; /* To equal area of circle with same size */ case 'd': p->symbol = DIAMOND; break; case 'E': /* Expect axis in km to be scaled based on -J */ p->convert_angles = TRUE; case 'e': p->symbol = ELLIPSE; p->n_required = 3; check = FALSE; break; case 'f': /* Fronts: -Sf/[dir][type][:] */ p->symbol = FRONT; p->f.f_off = 0.0; strcpy (text_cp, text); if ((c = strchr (text_cp, ':'))) { /* Gave :, set it and strip it off */ c++; /* SKip over the colon */ p->f.f_off = GMT_convert_units (c, GMT_INCH); c--; /* Go back to colon */ *c = 0; /* Effectively chops off the offset modifier */ } len = strlen (text_cp) - 1; old_style = FALSE; switch (text_cp[len]) { case 'f': /* Fault front */ p->f.f_symbol = F_FAULT; len--; break; case 't': /* Triangle front */ p->f.f_symbol = F_TRIANGLE; len--; break; case 's': /* Strike-slip front */ p->f.f_symbol = F_SLIP; len--; break; case 'c': /* [half-]circle front */ p->f.f_symbol = F_CIRCLE; len--; break; case 'b': /* [half-]square front */ p->f.f_symbol = F_BOX; len--; break; /* Old style (backward compatibility) */ case 'L': /* Left triangle */ p->f.f_symbol = F_TRIANGLE; case 'l': /* Left ticked fault */ p->f.f_sense = F_LEFT; old_style = TRUE; break; case 'R': /* Right triangle */ p->f.f_symbol = F_TRIANGLE; case 'r': /* Right ticked fault */ p->f.f_sense = F_RIGHT; old_style = TRUE; break; default: p->f.f_sense = F_CENTERED; break; } if (!old_style) { switch (text_cp[len]) { /* Get sense - default is centered */ case 'l': p->f.f_sense = F_LEFT; break; case 'r': p->f.f_sense = F_RIGHT; break; default: len++; p->f.f_sense = F_CENTERED; if (p->f.f_symbol == F_SLIP) { fprintf (stderr, "%s: Error in Option -Sf: Must specify (l)eft-lateral or (r)ight-lateral slip\n", GMT_program); exit (EXIT_FAILURE); } break; } } text_cp[len] = 0; /* Gets rid of the [dir][type] flags, if present */ /* Pull out and get spacing and size */ sscanf (&text_cp[1], "%[^/]/%s", txt_a, txt_b); p->f.f_gap = GMT_convert_units (txt_a, GMT_INCH); p->f.f_len = GMT_convert_units (txt_b, GMT_INCH); break; case 'H': p->size_x *= 1.09963611078; /* To equal area of circle with same size */ case 'h': p->symbol = HEXAGON; break; case 'I': p->size_x *= 1.55512030155; /* To equal area of circle with same size */ case 'i': p->symbol = ITRIANGLE; break; case 'l': p->symbol = TEXT; /* Look for a slash that separates size and string: */ for (j = 1, slash = 0; text_cp[j] && !slash; j++) if (text_cp[j] == '/') slash = j; /* Set j to the first char in the string: */ j = slash + 1; /* Copy string characters */ k = 0; while (text_cp[j] && text_cp[j] != ' ' && k < 63) p->string[k++] = text_cp[j++]; if (!k) { fprintf (stderr, "%s: GMT SYNTAX ERROR -Sl option: No string given\n", GMT_program); decode_error++; } p->string[k] = 0; break; case 'o': p->symbol = COLUMN; if (bset) p->base = atof (&text[bset+1]); break; case 'p': p->symbol = POINT; check = FALSE; break; case 'S': p->size_x *= 1.25331413731; /* To equal area of circle with same size */ case 's': p->symbol = SQUARE; break; case 'T': p->size_x *= 1.55512030155; /* To equal area of circle with same size */ case 't': p->symbol = TRIANGLE; break; case 'u': p->symbol = CUBE; break; case 'V': p->convert_angles = TRUE; case 'v': p->symbol = VECTOR; for (j = 1; text[j] && text[j] != 'n'; j++); len = strlen(text) - 1; if (text[j]) { /* Normalize option used */ k = GMT_get_unit (text[len]); if (k >= 0) p->u = k; p->v_norm = atof (&text[j+1]); if (p->v_norm > 0.0) { p->v_shrink = 1.0 / p->v_norm; p->symbol = VECTOR2; } text[j] = 0; /* Chop off the shrink part */ } if (text[1]) { /* It is possible that the user have appended a unit modifier after * the argument (which here are vector attributes). We use that * to set the unit, but only if the vector attributes themselves have * units. (If not we would override MEASURE_UNIT without cause). * So, -SV0.1i/0.2i/0.3ic will expect 4th colum to have length in cm * while SV0.1i/0.2i/0.3i expects data units in MEASURE_UNIT */ if (isalpha ((int)text[len]) && isalpha ((int)text[len-1])) { p->u = GMT_get_unit (text[len]); text[len] = 0; } sscanf (&text[1], "%[^/]/%[^/]/%s", txt_a, txt_b, txt_c); p->v_width = GMT_convert_units (txt_a, GMT_INCH); p->h_length = GMT_convert_units (txt_b, GMT_INCH); p->h_width = GMT_convert_units (txt_c, GMT_INCH); } if (p->symbol == VECTOR2) text[j] = 'n'; /* Put back the n part */ p->read_vector = TRUE; p->n_required = 2; check = FALSE; break; case 'w': p->symbol = PIE; p->n_required = 2; break; case 'x': p->symbol = CROSS; break; default: decode_error = TRUE; fprintf (stderr, "%s: GMT SYNTAX ERROR -S option: Unrecognized symbol type %c\n", GMT_program, symbol_type); break; } if (p->given_size_x == 0.0 && check) { p->read_size = TRUE; p->n_required++; } else p->read_size = FALSE; if (!cmd && p->symbol == COLUMN) { if (!bset) p->base = (project_info.xyz_projection[2] == LOG10) ? 1.0 : 0.0; } if (!cmd && p->symbol == BAR) { if (!bset) p->base = (project_info.xyz_projection[1] == LOG10) ? 1.0 : 0.0; } return (decode_error); } #define GMT_ELLIPSE_APPROX 72 void GMT_plot_ellipse (double lon, double lat, double z, double major, double minor, double azimuth, struct GMT_FILL fill, int outline) { /* GMT_plot_ellipse takes the location, axes (in km), and azimuth of an ellipse and draws the ellipse using the chosen map projection */ int i, plot_n, start, n_use; double delta_azimuth, sin_azimuth, cos_azimuth, sinp, cosp, angle, x, y, x_prime, y_prime, rho, c; double sin_c, cos_c, keep_z; double *px, *py, *xp, *yp; keep_z = project_info.z_level; project_info.z_level = z; px = (double *) GMT_memory (VNULL, (size_t)GMT_ELLIPSE_APPROX, sizeof (double), GMT_program); py = (double *) GMT_memory (VNULL, (size_t)GMT_ELLIPSE_APPROX, sizeof (double), GMT_program); delta_azimuth = 2.0 * M_PI / GMT_ELLIPSE_APPROX; major *= 1000.0; minor *= 1000.0; /* Convert to meters */ azimuth = 90.0 - azimuth; /* Because the code below originally used directions instead */ azimuth *= D2R; sincos (azimuth, &sin_azimuth, &cos_azimuth); sincos (lat * D2R, &sinp, &cosp); /* Set up azimuthal equidistant projection */ /* Approximate ellipse by a n-sided polygon */ for (i = 0; i < GMT_ELLIPSE_APPROX; i++) { angle = i * delta_azimuth; sincos (angle, &y, &x); x *= major; y *= minor; /* Get rotated coordinates in km */ x_prime = x * cos_azimuth - y * sin_azimuth; y_prime = x * sin_azimuth + y * cos_azimuth; /* Convert km back to lon lat */ rho = hypot (x_prime, y_prime); c = rho / project_info.EQ_RAD; sincos (c, &sin_c, &cos_c); py[i] = d_asin (cos_c * sinp + (y_prime * sin_c * cosp / rho)) * R2D; if (lat == 90.0) px[i] = lon + R2D * d_atan2 (x_prime, -y_prime); else if (lat == -90.0) px[i] = lon + R2D * d_atan2 (x_prime, y_prime); else px[i] = lon + R2D * d_atan2 (x_prime * sin_c, (rho * cosp * cos_c - y_prime * sinp * sin_c)); while ((px[i] - project_info.central_meridian) < -180.0) px[i] += 360.0; while ((px[i] - project_info.central_meridian) > 180.0) px[i] -= 360.0; } if ((plot_n = GMT_clip_to_map (px, py, GMT_ELLIPSE_APPROX, &xp, &yp)) == 0) return; while (plot_n > GMT_n_alloc) GMT_get_plot_array (); memcpy ((void *)GMT_x_plot, (void *)xp, (size_t)(plot_n * sizeof (double))); memcpy ((void *)GMT_y_plot, (void *)yp, (size_t)(plot_n * sizeof (double))); if ((*GMT_will_it_wrap) (xp, yp, plot_n, &start)) { /* Polygon wraps */ /* First truncate agains left border */ GMT_n_plot = (*GMT_truncate) (xp, yp, plot_n, start, -1); n_use = GMT_compact_line (GMT_x_plot, GMT_y_plot, GMT_n_plot, FALSE, 0); GMT_2D_to_3D (GMT_x_plot, GMT_y_plot, GMT_n_plot); GMT_fill (GMT_x_plot, GMT_y_plot, n_use, &fill, outline); /* Then truncate agains right border */ GMT_n_plot = (*GMT_truncate) (xp, yp, plot_n, start, +1); n_use = GMT_compact_line (GMT_x_plot, GMT_y_plot, GMT_n_plot, FALSE, 0); GMT_2D_to_3D (GMT_x_plot, GMT_y_plot, GMT_n_plot); GMT_fill (GMT_x_plot, GMT_y_plot, n_use, &fill, outline); } else { GMT_2D_to_3D (xp, yp, plot_n); GMT_fill (xp, yp, plot_n, &fill, outline); } project_info.z_level = keep_z; GMT_free ((void *)xp); GMT_free ((void *)yp); GMT_free ((void *)px); GMT_free ((void *)py); } void GMT_draw_fence (double x[], double y[], int n, struct FRONTLINE *f, struct GMT_FILL *g, BOOLEAN outline) { int i, ngap; BOOLEAN skip; double *s, xx[4], yy[4], dist, w, frac, dx, dy, angle, dir1, dir2; double gap, x0, y0, xp, yp, len2, len3, cosa, sina, sa, ca, offx, offy; if (n < 2) return; s = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program); for (i = 1, s[0] = 0.0; i < n; i++) { /* Watch out for longitude wraps */ dx = x[i] - x[i-1]; w = GMT_half_map_width (y[i]); if (GMT_world_map && dx > w) dx = copysign (2 * w - fabs (dx), -dx); s[i] = s[i-1] + hypot (dx, y[i] - y[i-1]); } if (f->f_gap > 0.0) { ngap = irint (s[n-1] / f->f_gap); gap = s[n-1] / (double)ngap; dist = f->f_off; ngap++; } else { ngap = (int) fabs (f->f_gap); gap = s[n-1] / (ngap - 1); if (ngap == 1) dist = 0.5 * s[n-1]; } len2 = 0.5 * f->f_len; len3 = 0.866025404 * f->f_len; if (f->f_sense == F_CENTERED) len3 = len2; i = 0; while (i < n) { while ((s[i] - dist) > -SMALL) { /* Time for tick */ if (i > 0) { dx = x[i] - x[i-1]; dy = y[i] - y[i-1]; } else { dx = x[1] - x[0]; dy = y[1] - y[0]; } if (fabs (dist - s[i]) < SMALL) { x0 = x[i]; y0 = y[i]; } else { frac = (s[i] - dist) / (s[i] - s[i-1]); x0 = x[i] - dx * frac; y0 = y[i] - dy * frac; } angle = d_atan2 (dy, dx); skip = (GMT_world_map && fabs (dx) > GMT_half_map_width (y[i])); /* Dont do ticks on jumps */ if (skip) { dist += gap; i++; continue; } switch (f->f_symbol) { case F_TRIANGLE: /* Triangle */ switch (f->f_sense) { case F_CENTERED: sincos (angle, &sina, &cosa); xx[0] = x0 + len2 * cosa; yy[0] = y0 + len2 * sina; xx[1] = x0 - len3 * sina; yy[1] = y0 + len3 * cosa; xx[2] = x0 - len2 * cosa; yy[2] = y0 - len2 * sina; xx[3] = x0 + len3 * sina; yy[3] = y0 - len3 * cosa; if (project_info.three_D) GMT_2D_to_3D (xx, yy, 4); ps_patch (xx, yy, 4, g->rgb, outline); break; case F_RIGHT: angle += M_PI; case F_LEFT: sincos (angle, &sina, &cosa); xx[0] = x0 + len2 * cosa; yy[0] = y0 + len2 * sina; xx[1] = x0 - len3 * sina; yy[1] = y0 + len3 * cosa; xx[2] = x0 - len2 * cosa; yy[2] = y0 - len2 * sina; if (project_info.three_D) GMT_2D_to_3D (xx, yy, 3); ps_patch (xx, yy, 3, g->rgb, outline); break; } break; case F_CIRCLE: /* Circles */ switch (f->f_sense) { case F_CENTERED: circle3D (x0, y0, 0.0, f->f_len, g->rgb, outline); break; case F_RIGHT: angle += M_PI; case F_LEFT: dir1 = angle; dir2 = angle + M_PI; if (dir1 > dir2) dir1 -= TWO_PI; pie3D (x0, y0, 0.0, f->f_len, dir1, dir2, g->rgb, outline); break; } break; case F_BOX: /* Squares */ switch (f->f_sense) { case F_CENTERED: /* Full square centered on line */ sincos (angle, &sina, &cosa); xx[0] = x0 + len2 * (cosa + sina); /* LR */ yy[0] = y0 + len2 * (sina - cosa); xx[1] = x0 + len2 * (cosa - sina); /* UR */ yy[1] = y0 + len2 * (sina + cosa); xx[2] = x0 + len2 * (-cosa - sina); /* UL */ yy[2] = y0 + len2 * (-sina + cosa); xx[3] = x0 + len2 * (-cosa + sina); /* LL */ yy[3] = y0 + len2 * (-sina - cosa); break; case F_RIGHT: angle += M_PI; case F_LEFT: /* Half square on the chosen side */ sincos (angle, &sina, &cosa); xx[0] = x0 + len2 * (cosa); /* LR */ yy[0] = y0 + len2 * (sina); xx[1] = x0 + len2 * (cosa - sina); /* UR */ yy[1] = y0 + len2 * (sina + cosa); xx[2] = x0 + len2 * (-cosa - sina); /* UL */ yy[2] = y0 + len2 * (-sina + cosa); xx[3] = x0 + len2 * (-cosa); /* LL */ yy[3] = y0 + len2 * (-sina); break; } if (project_info.three_D) GMT_2D_to_3D (xx, yy, 4); ps_patch (xx, yy, 4, g->rgb, outline); break; case F_SLIP: /* draw strike-slip arrows */ sincos (angle, &sina, &cosa); offx = gmtdefs.anot_offset * sina; /* get offsets from front line */ offy = gmtdefs.anot_offset * cosa; /* sense == F_LEFT == left-lateral, R_RIGHT = right lateral */ /* arrow "above" line */ sincos (angle + (f->f_sense * 30.0 * D2R), &sa, &ca); xp = x0 + f->f_sense * offx; yp = y0 - f->f_sense * offy; xx[0] = xp - len2 * cosa; yy[0] = yp - len2 * sina; xx[1] = xp + len2 * cosa; yy[1] = yp + len2 * sina; xx[2] = xx[1] - len2 * ca; yy[2] = yy[1] - len2 * sa; if (project_info.three_D) GMT_2D_to_3D (xx, yy, 3); ps_line (xx, yy, 3, 3, FALSE, TRUE); /* arrow "below" line */ sincos (angle - (f->f_sense * 150.0 *D2R), &sa, &ca); xp = x0 - f->f_sense * offx; yp = y0 + f->f_sense * offy; xx[0] = xp + len2 * cosa; yy[0] = yp + len2 * sina; xx[1] = xp - len2 * cosa; yy[1] = yp - len2 * sina; xx[2] = xx[1] - len2 * ca; yy[2] = yy[1] - len2 * sa; if (project_info.three_D) GMT_2D_to_3D (xx, yy, 3); ps_line (xx, yy, 3, 3, FALSE, TRUE); break; case F_FAULT: /* Normal fault ticks */ xx[0] = xx[1] = x0; yy[0] = yy[1] = y0; if (f->f_sense == F_CENTERED) { angle -= M_PI_2; sincos (angle, &sina, &cosa); xx[0] += len2 * cosa; yy[0] += len2 * sina; xx[1] -= len2 * cosa; yy[1] -= len2 * sina; } else { angle += (f->f_sense * M_PI_2); sincos (angle, &sina, &cosa); xx[1] += len2 * cosa; yy[1] += len2 * sina; } if (project_info.three_D) GMT_2D_to_3D (xx, yy, 2); ps_line (xx, yy, 2, 3, FALSE, TRUE); break; } dist += gap; } i++; } GMT_free ((void *)s); }