/*--------------------------------------------------------------------
 *	$Id: psxyz.c,v 1.63 2004/07/08 23:23:08 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 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 <x,y,z> 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 <x,y,z,size> 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:   4
 *
 */

#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 RECT		17
#define XDASH		18
#define YDASH		19
#define ZDASH		20
#define TEXT		21
#define PENTAGON	22
#define OCTAGON		23
#define CUSTOM		24
#define FRONT		-100
#define QUOTED_LINE	-200

/* 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 */

	int convert_angles;	/* If 2, convert azimuth to angle on map, 1 special case for -JX, 0 plain case */
	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 */
	int v_just;		/* How to justify vector: head point given (3), head (2), center(1), tail (0 - Default) */
	BOOLEAN v_double_heads;		/* If TRUE, Add 8 (|= 8) to outline to specify double-headed vector (FALSE is single-headed) */

	char string[64];	/* Character code to plot (could be octal) */

	struct FRONTLINE f;	/* parameters needed for a front */
	struct CUSTOM_SYMBOL *custom;	/* pointer to a custom symbol */

	struct GMT_CONTOUR G;	/* For labelled lines */
} S;


struct DATA1 {
	double lon, lat;
	double x, y, z, value, dist;
	float x_size, y_size, z_size;
	int symbol, flag;
	char *string;
	struct GMT_FILL f;
	struct GMT_PEN p;
} *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 bar3D(double x, double y, double z, double base, double size, 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);
void extract_label (char *line, char *label);

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, n_expected = 0, this_outline;
	
	BOOLEAN	error = FALSE, nofile = TRUE, polygon = FALSE, done, close_polygon = FALSE, penset_OK = TRUE;
	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, font_size;
	double x2, y2, v_w, h_l, h_w, dx, dy, xx[2], yy[2], s, c;
	
	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 */
	GMT_contlabel_init (&S.G);

	for (i = 0; i < 9; i++) rgb[i] = -1;

	/* Initialize PSXYZ_SYMBOL structure */

	S.symbol = S.n_required = S.convert_angles  = 0;
	S.size_x = S.size_y = S.v_norm = S.v_shrink = 0.0;
	S.read_size = S.user_unit = 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 'b':
				case 'c':
				case 'f':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
			
				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 <xyzfiles> -R<west/east/south/north/zmin/zmax> -J<params>\n");
		fprintf(stderr, "\t[-Jz<params>] [-B<tickinfo>] [-C<cpt>] [-E<az/el>] [-G<fill>] [-H[<nrec>]] [-K] [-L]\n");
		fprintf(stderr, "\t[-M[<flag>]] [-N] [-O] [-P] [-S<symbol><size>[/size_y]] [-U] [-V] [-W[<pen>]] [-X<x_shift>]\n");
		fprintf(stderr, "\t[-Y<y_shift>] [-:] [-c<ncopies>] [-bi[s][<n>]] [-f[i|o]<colinfo>]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<xyzfiles> is one or more files.  If none, read standard input\n");
		GMT_explain_option ('j');
		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 lines/polygons and looks for -Z<tval> 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) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
		fprintf (stderr, "\t   2) p[or P]<dpi>/<pattern> 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 ('Z');
		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   -(xdash), st(a)r, (b)ar, (c)ircle, (d)iamond, (e)llipse, (f)ront, octa(g)on, (h)exagon\n");
		fprintf (stderr, "\t   (i)nvtriangle, (k)ustom, (%c)etter, pe(n)tagon, c(o)lumn, (p)oint, (q)uoted line, (r)ect, (s)quare\n", 'l');
		fprintf (stderr, "\t   (t)riangle, c(u)be, (v)ector, (w)edge, (x)cross, (y)dash, (z)dash\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, F, H, I, N, S, T are adjusted to have same area of a circle of given size\n");
		fprintf (stderr, "\t   Bar (or Column): Append b<base> 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.  If projection is linear then\n");
		fprintf (stderr, "\t     we scale the minor/major axes by the map scale.\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 <dir> 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   Kustom: append <symbolname> immediately after 'k'; this will look for <symbolname>.def in\n");
		fprintf (stderr, "\t     the current directory or in $GMTHOME/share/custom.\n");
		GMT_list_custom_symbols ();
		fprintf (stderr, "\t   Letter: append /<string> after symbol size, and optionally %%<font>\n");
		fprintf (stderr, "\t   Quoted line: Give [d|f|n|l|x]<info>[:<labelinfo>]. z must be constant for each segment\n");
		fprintf (stderr, "\t     <code><info> controls the placement of labels along the lines.  Choose among five algorithms:\n");
		GMT_cont_syntax (7, 1);
		fprintf (stderr, "\t   <labelinfo> controls the specifics of the labels.  Append what you need:\n");
		GMT_label_syntax (7, 1);
		fprintf (stderr, "\t   Rectangles: the x- and y-dimension must be in input columns 4 and 5.\n");
		fprintf (stderr, "\t   Vectors: the direction and length must be in input columns 4 and 5.\n");
		fprintf (stderr, "\t     Furthermore, <size> means arrowwidth/headlength/headwith [Default is %g/%g/%g]\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");
		fprintf (stderr, "\t     Insert h(head), +(center), or t(ail) after -Sv|V to justify vector w.r.t. input (x,y).\n");
		fprintf (stderr, "\t     Insert s(egment) if (x,y) is tail and columns 3 and 4 holds head (x,y).\n");
		fprintf (stderr, "\t     Upper case H, B, T, S will give double-headed vector [Default is t]\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 = %gp, color = black, solid line].\n", pen.width);
		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 ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	/* Check that the options selected are mutually consistent */
	
	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 || S.symbol == CUSTOM)) {	/* 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[GMT_IN] && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && 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 || S.symbol == QUOTED_LINE) 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[GMT_IN]) ? GMT_io.ncol[GMT_IN] : n_required;
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN]=  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 .gmtcommands4 */
	
	if (GMT_io.binary[GMT_IN] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[GMT_IN], type[GMT_io.single_precision[GMT_IN]]);
	}

	if (cpt_given) GMT_read_cpt (cpt);
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	if (S.symbol == QUOTED_LINE) {
		if (GMT_contlabel_prep (&S.G, NULL, 0)) exit (EXIT_FAILURE);
		penset_OK = FALSE;	/* Since it is set in PSL */
	}
	
	GMT_map_setup (west, east, south, north);

	if ((cpt_given || 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, gmtdefs.encoding.name, 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;

	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;

	if (penset_OK) GMT_setpen (&pen);

	if (S.symbol == TEXT && fill_set && !outline) ps_setpaint (fill.rgb);
	if (S.symbol == TEXT) {
		font_size = S.size_x * 72.0;	/* To get points */
		ps_setfont (S.font_no);		/* Set the required font */
	}
	if (S.convert_angles && MAPPING) S.convert_angles = 2;
	
	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 (GMT_IN);
#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;
		}

		GMT_io.file_no++;
		GMT_io.seg_no = -1;
		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[GMT_IN]) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (line, BUFSIZ, fp);
				if (i == 0 && S.G.label_type == 2) extract_label (&line[1], S.G.label);	/* Set initial header as potential label */
			}
		}		
		
		if (S.symbol > 0) {	/* symbol part (not counting FRONT and QUOTED_LINE) */
		
			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")) || (p = strstr (GMT_io.segment_header, "	-G"))) {
						strcpy (line, &p[3]);
						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")) || (p = strstr (GMT_io.segment_header, "	-Z")))) {	/* Set symbol r/g/b via cpt-lookup */
						z = (p[3] == 'N') ? GMT_d_NaN : atof (&p[3]);
						GMT_get_rgb24 (z, fill.rgb);
					}
					else
						fill = default_fill;

					if ((p = strstr (GMT_io.segment_header, " -W")) || (p = strstr (GMT_io.segment_header, "	-W"))) {
						strcpy (line, &p[3]);
						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 {
						pen = 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;
					GMT_io.current_record[i--] = 0;	/* Chop off the \n */
					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) */
					data1[n].flag = S.convert_angles;	/* Azim vs Dir. */
					if (data1[i].y_size < S.v_norm) data1[n].flag |= 8;	/* Flag to shrink vector attributes */
				}
				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 */
					data1[n].flag = S.convert_angles;	/* Azim vs Dir. */
					if (S.convert_angles == 1) {	/* Got axes in user units, change to inches */
						data1[n].y_size *= (float)project_info.x_scale;
						data1[n].z_size *= (float)project_info.x_scale;
					}
					else if (S.convert_angles == 0){	/* Got axes in current plot 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 if (S.symbol == RECT) {
					data1[n].x_size = (float)(0.5 * in[three+S.read_size]);	/* x-dim */
					data1[n].y_size = (float)(0.5 * in[four+S.read_size]);	/* y-dim */
				}
				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 && S.v_just < 3 && 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) {	/* Lookup z to get rgb */
					data1[n].value = in[3];
					GMT_get_rgb24 (data1[i].value, fill.rgb);
				}
				if (S.symbol == COLUMN) {
					if (GMT_is_dnan (S.base))
						tmp = 0.0;
					else
						GMT_z_to_zz (S.base, &tmp);
					data1[n].z_size = (float)tmp;
				}
				if (S.symbol == BAR) {
					if (GMT_is_dnan (S.base))
						tmp = 0.0;
					else
						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;
				data1[n].f = fill;
				data1[n].p = pen;

				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 (data1[i].symbol == COLUMN || data1[i].symbol == CUBE) {
					for (j = 0; j < 3; j++) {
						rgb[3*j] = data1[i].f.rgb[0];
						rgb[3*j+1] = data1[i].f.rgb[1];
						rgb[3*j+2] = data1[i].f.rgb[2];
						GMT_illuminate (lux[j], &rgb[3*j]);
					}
				}
				
				if (GMT_cpt_skip) continue;	/* Skip this because cpt says so */

				if (outline) GMT_setpen (&data1[i].p);
				
				switch (data1[i].symbol) {
					case NONE:
						break;
					case STAR:
						GMT_star3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.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, data1[i].f.rgb, outline);
						}
						else
							bar3D (data1[i].x, data1[i].y, data1[i].z, data1[i].z_size, (double)data1[i].x_size, data1[i].f.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:
						GMT_cross3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size);
						break;
					case POINT:
						GMT_cross3D (data1[i].x, data1[i].y, data1[i].z, POINTSIZE);
						break;
					case CIRCLE:
						GMT_circle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case SQUARE:
						GMT_square3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case HEXAGON:
						GMT_hexagon3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case PENTAGON:
						GMT_pentagon3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case OCTAGON:
						GMT_octagon3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case TRIANGLE:
						GMT_triangle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case ITRIANGLE:
						GMT_itriangle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.rgb, outline);
						break;
					case DIAMOND:
						GMT_diamond3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, data1[i].f.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 == 2)
							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
							GMT_ellipse3D (data1[i].x, data1[i].y, data1[i].z, data1[i].x_size, data1[i].y_size, data1[i].z_size, data1[i].f.rgb, outline);
						break;
					case TEXT:
						font_size = 2.0 * data1[i].x_size * 72.0;	/* To get points */
						if (outline && fill_set) {
							ps_setpaint (data1[i].f.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 & 2) {
							GMT_azim_to_angle (data1[i].lon, data1[i].lat, 0.1, data1[i].x_size, &tmp);
							data1[i].x_size = (float)tmp;
						}
						if (S.v_just == 3) {
							GMT_geo_to_xy (in[3], in[4], &x2, &y2);
							if (GMT_is_dnan (x2) || GMT_is_dnan (y2)) {
								fprintf (stderr, "%s: Warning: Vector head coordinates contain NaNs near line %d. Skipped\n", GMT_program, n_total_read);
								continue;
							}
						}
						else {
							sincos (data1[i].x_size * D2R, &s, &c);
							x2 = data1[i].x + data1[i].y_size * c;
							y2 = data1[i].y + data1[i].y_size * s;
							if (S.v_just) {
								dx = S.v_just * 0.5 * (x2 - data1[i].x);	dy = S.v_just * 0.5 * (y2 - data1[i].y);
								data1[i].x -= dx;	data1[i].y -= dy;
								x2 -= dx;		y2 -= dy;
							}
						}
						this_outline = (S.v_double_heads) ? outline + 8 : outline;
						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, data1[i].f.rgb, this_outline);
						break;
					case VECTOR2:
						if (data1[i].y_size <= 0.0) continue;
						if (data1[i].flag & 2) {
							GMT_azim_to_angle (data1[i].lon, data1[i].lat, 0.1, data1[i].x_size, &tmp);
							data1[i].x_size = (float)tmp;
						}
						if (S.v_just == 3) {
							GMT_geo_to_xy (in[3], in[4], &x2, &y2);
							if (GMT_is_dnan (x2) || GMT_is_dnan (y2)) {
								fprintf (stderr, "%s: Warning: Vector head coordinates contain NaNs near line %d. Skipped\n", GMT_program, n_total_read);
								continue;
							}
						}
						else {
							sincos (data1[i].x_size * D2R, &s, &c);
							x2 = data1[i].x + data1[i].y_size * c;
							y2 = data1[i].y + data1[i].y_size * s;
							if (S.v_just) {
								dx = S.v_just * 0.5 * (x2 - data1[i].x);	dy = S.v_just * 0.5 * (y2 - data1[i].y);
								data1[i].x -= dx;	data1[i].y -= dy;
								x2 -= dx;		y2 -= dy;
							}
						}
						if (data1[i].flag & 8) {	/* 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;
							this_outline = (S.v_double_heads) ? outline + 8 : outline;
							GMT_vector3D (data1[i].x, data1[i].y, x2, y2, data1[i].z, v_w, h_l, h_w, gmtdefs.vector_shape, data1[i].f.rgb, this_outline);
						}
						else {	/* Leave as specified */
							this_outline = (S.v_double_heads) ? outline + 8 : outline;
							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, data1[i].f.rgb, this_outline);
						}
						break;
					case PIE:
						GMT_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, data1[i].f.rgb, outline);
						break;
					case XDASH:
						GMT_xyz_to_xy (data1[i].x-(double)data1[i].x_size, data1[i].y, data1[i].z, &xx[0], &yy[0]);
						GMT_xyz_to_xy (data1[i].x+(double)data1[i].x_size, data1[i].y, data1[i].z, &xx[1], &yy[1]);
						ps_segment (xx[0], yy[0], xx[1], yy[1]);
						break;
					case YDASH:
						GMT_xyz_to_xy (data1[i].x, data1[i].y-(double)data1[i].x_size, data1[i].z, &xx[0], &yy[0]);
						GMT_xyz_to_xy (data1[i].x, data1[i].y+(double)data1[i].x_size, data1[i].z, &xx[1], &yy[1]);
						ps_segment (xx[0], yy[0], xx[1], yy[1]);
						break;
					case ZDASH:
						GMT_xyz_to_xy (data1[i].x, data1[i].y, data1[i].z-(double)data1[i].x_size, &xx[0], &yy[0]);
						GMT_xyz_to_xy (data1[i].x, data1[i].y, data1[i].z+(double)data1[i].x_size, &xx[1], &yy[1]);
						ps_segment (xx[0], yy[0], xx[1], yy[1]);
						break;
					case RECT:
						GMT_rect3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, (double)data1[i].y_size, data1[i].f.rgb, outline);
						break;
					case CUSTOM:
						GMT_draw_custom_symbol (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, S.custom, &pen, &fill, 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, " -W")) || (p = strstr (GMT_io.segment_header, "	-W"))) {
						strcpy (line, &p[3]);
						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 {
							if (penset_OK) GMT_setpen (&pen);
							outline = TRUE;
						}
					}
					else {
						if (penset_OK) GMT_setpen (&default_pen);
						outline = default_outline;
						pen = default_pen;
					}
					if ((p = strstr (GMT_io.segment_header, " -G")) || (p = strstr (GMT_io.segment_header, "	-G"))) {
						strcpy (line, &p[3]);
						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")) || (p = strstr (GMT_io.segment_header, "	-Z")))) {	/* Set polygon r/g/b via cpt-lookup */
						if(!strncmp (&p[3], "NaN", 3)) {	/* Got -ZNaN */
							GMT_get_rgb24 (GMT_d_NaN, fill.rgb);
							fill.use_pattern = FALSE;
							polygon = TRUE;
						}
						else {
							if (sscanf (&p[3], "%lg", &z) == 1) {
								GMT_get_rgb24 (z, fill.rgb);
								fill.use_pattern = FALSE;
								polygon = TRUE;
							}
						}
						if (!close_polygon) {	/* Use -Zz -Ccpt to set pen color instead */
							polygon = FALSE;
							ps_setpaint (fill.rgb);
						}
					}
					else {
						polygon = default_polygon;
						fill = default_fill;
					}
					if (S.G.label_type == 2) {	/* Update segment header */
						extract_label (&GMT_io.segment_header[1], S.G.label);
					}

					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);
				
					if (polygon) {
						for (i = 0; i < n; i++) GMT_geoz_to_xy (data2[i].x, data2[i].y, data2[i].z, &xp[i], &yp[i]);
						GMT_fill (xp, yp, n, &fill, outline);
					}
					else if (S.symbol == QUOTED_LINE) {
						for (i = 0; i < n; i++) {
							GMT_geo_to_xy (data2[i].x, data2[i].y, &xp[i], &yp[i]);
						}
						S.G.z_level = data2[0].z;	/* Any z will do since we require line to be in x-y plate */
						S.G.line_pen = pen;
						GMT_hold_contour (&xp, &yp, n, "N/A", 'A', S.G.label_angle, close_polygon, &S.G);
					}
					else {
						for (i = 0; i < n; i++) GMT_geoz_to_xy (data2[i].x, data2[i].y, data2[i].z, &xp[i], &yp[i]);
						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 (S.symbol == QUOTED_LINE) GMT_contlabel_plot (&S.G);
	
	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 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 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, one, colon;
	BOOLEAN check, old_style;
	char symbol_type, txt_a[32], txt_b[32], txt_c[32], text_cp[128], *c;

	p->n_required = p->convert_angles = 0;
	p->user_unit = p->shrink = p->read_vector = FALSE;

	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_font_lookup (c, GMT_font, N_FONTS);
			if (p->font_no >= N_FONTS) {
				fprintf (stderr, "%s: -Sl contains bad font (set to %s (0))\n", GMT_program, GMT_font[gmtdefs.annot_font[0]].name);
				p->font_no = gmtdefs.annot_font[0];
			}
		}
		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 (text[0] == 'k') {	/* Custom symbol spec */
		for (j = strlen (text); j > 0 && text[j] != '/'; j--);;
		if (j == 0) {	/* No slash, i.e., no symbol size given */
			if (p->size_x == 0.0) p->size_x = p->given_size_x;
			n = sscanf (text, "%c%s", &symbol_type, text_cp);
		}
		else {
			text[j] = ' ';
			n = sscanf (text, "%c%s %s", &symbol_type, text_cp, txt_a);
			text[j] = '/';
			p->given_size_x = p->size_x = GMT_convert_units (txt_a, GMT_INCH);
		}
	}
	else if (strchr ("-aAbCcDdeEfGgHhIiNnopqrSsTtuVvwxy", (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 ("-aAbCcDdeEfGgHhIiNnopqrSsTtuVvwxy", (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<size_x>[c|i|m|p|u][b<base>]				*/
		/* Column:	-So<size_x>[c|i|m|p][/<ysize>[c|i|m|p]][u][b<base>]	*/
		/* Cube:	-Su<size_x>[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<base> 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 '-':
			p->symbol = XDASH;
			p->size_x *= 0.5;		/* We will use +- to get full width */
			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 = 1;
		case 'e':
			p->symbol = ELLIPSE;
			p->n_required = 3;
			check = FALSE;
			break;

		case 'f':	/* Fronts:   -Sf<spacing>/<size>[dir][type][:<offset>]	*/
			p->symbol = FRONT;
			p->f.f_off = 0.0;
			strcpy (text_cp, text);
			if ((c = strchr (text_cp, ':'))) {	/* Gave :<offset>, 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 (GMTMANSECTION)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 'G':
			p->size_x *= 0.94805944897;	/* To equal area of circle with same size */
		case 'g':
			p->symbol = OCTAGON;
			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 'N':
			p->size_x *= 1.55512030155;	/* To equal area of circle with same size */
		case 'n':
			p->symbol = PENTAGON;
			break;
		case 'o':
			p->symbol = COLUMN;
			if (bset) p->base = atof (&text[bset+1]);
			break;
		case 'p':
			p->symbol = POINT;
			check = FALSE;
			break;
		case 'q':	/* Quoted lines: -Sq[d|n|l|x]<info>[:<labelinfo>] */
			p->symbol = QUOTED_LINE;
			for (j = 1, colon = 0; text[j]; j++) if (text[j] == ':') colon = j;
			if (colon) {	/* Gave :<labelinfo> */
				text[colon] = 0;
				decode_error += GMT_contlabel_info ('S', &text[1], &p->G);
				decode_error += GMT_contlabel_specs (&text[colon+1], &p->G);
			}
			else
				decode_error += GMT_contlabel_info ('S', &text[1], &p->G);
			break;
			check = FALSE;
		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 = 1;
		case 'v':
			p->symbol = VECTOR;
			switch (text[1]) {	/* Check if s(egment), h(ead), b(alance center), or t(ail) have been specified */
				case 'S':	/* Input (x,y) refers to vector head (the tip), double heads */
					p->v_double_heads = TRUE;
				case 's':	/* Input (x,y) refers to vector head (the tip), single head  */
					p->v_just = 3;
					one = 2;
					break;
				case 'H':	/* Input (x,y) refers to vector head (the tip), double heads */
					p->v_double_heads = TRUE;
				case 'h':	/* Input (x,y) refers to vector head (the tip), single head */
					p->v_just = 2;
					one = 2;
					break;
				case 'B':	/* Input (x,y) refers to balance point of vector, double heads */
					p->v_double_heads = TRUE;
				case 'b':	/* Input (x,y) refers to balance point of vector, single head */
					p->v_just = 1;
					one = 2;
					break;
				case 'T':	/* Input (x,y) refers to tail of vector, double heads */
					p->v_double_heads = TRUE;
				case 't':	/* Input (x,y) refers to tail of vector [Default], single head */
					p->v_just = 0;
					one = 2;
					break;
				default:	/* No modifier given, default to tail, single head */
					p->v_just = 0;
					one = 1;
					break;
			}
			for (j = one; text[j] && text[j] != 'n'; j++);
			len = strlen(text) - 1;
			if (text[j] == 'n') {	/* 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[one]) {
				/* It is possible that the user have appended a unit modifier after
				 * the <size> 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 column 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[one], "%[^/]/%[^/]/%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<shrink> part */
			p->read_vector = TRUE;
			p->n_required = 2;
			check = FALSE;
			break;
		case 'w':
			p->symbol = PIE;
			p->n_required = 2;
			break;
		case 'r':
			p->symbol = RECT;
			p->n_required = 2;
			check = FALSE;
			break;
		case 'x':
			p->symbol = CROSS;
			break;
		case 'y':
			p->symbol = YDASH;
			p->size_x *= 0.5;		/* We will use +- to get full width */
			break;
		case 'z':
			p->symbol = ZDASH;
			p->size_x *= 0.5;		/* We will use +- to get full width */
			break;
		case 'k':
			p->symbol = CUSTOM;
			p->custom = GMT_get_custom_symbol (text_cp);
			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; 
	double *px, *py, *xp, *yp;

	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_2Dz_to_3D (GMT_x_plot, GMT_y_plot, z, 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_2Dz_to_3D (GMT_x_plot, GMT_y_plot, z, GMT_n_plot);
		GMT_fill (GMT_x_plot, GMT_y_plot, n_use, &fill, outline);
						
	}
	else {
		GMT_2Dz_to_3D (xp, yp, z, plot_n);
		GMT_fill (xp, yp, plot_n, &fill, outline);
	}

	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]));	/* Don't 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, project_info.z_level, 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, project_info.z_level, 3);
							ps_patch (xx, yy, 3, g->rgb, outline);
							break;
					}
					break;
				
				case F_CIRCLE:	/* Circles */	
					switch (f->f_sense) {
						case F_CENTERED:
							GMT_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;
							GMT_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, project_info.z_level, 4);
					ps_patch (xx, yy, 4, g->rgb, outline);
					break;
				
				case F_SLIP: /* draw strike-slip arrows */
					sincos (angle, &sina, &cosa);
				        offx = gmtdefs.annot_offset[0] * sina; /* get offsets from front line */
					offy = gmtdefs.annot_offset[0] * 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, project_info.z_level, 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, project_info.z_level, 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, project_info.z_level, 2);
					ps_line (xx, yy, 2, 3, FALSE, TRUE);
					break;
			}
			dist += gap;
		}
		i++;
	}
	GMT_free ((void *)s);
}

void extract_label (char *line, char *label)
{
	int i = 0, j, j0;
	char *p;
	
	if ((p = strstr (line, " -L")) || (p = strstr (line, "	-L"))) 	/* Get label specified wih -L option */
		i = p - line;
	else {								/* Bypass whitespace and pick first word */
		while (line[i] && (line[i] == ' ' || line[i] == '\t')) i++;
		
	}
	if ((p = strchr (&line[i], '\"'))) {	/* Gave several words as label */
		for (j0 = j = (int)(p - line) + 1; line[j] != '\"'; j++);
		if (line[j] == '\"')	/* Found the matching quote */
			strncpy (label, &line[j0], j-j0);
		else {			/* Missing the matching quote */
			sscanf (&line[i], "%s", label);
			fprintf (stderr, "%s: Warning: Label (%s) not terminated by matching quote\n", GMT_program, label);
		}
	}
	else
		sscanf (&line[i], "%s", label);
}
