/*--------------------------------------------------------------------
 *	$Id: pswiggle.c,v 1.24 2004/08/17 01:35:45 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
 *--------------------------------------------------------------------*/
/*
 * pswiggle reads x,y,z from GMT_stdin and plots a wiggleplot using the
 * specified map projection. If the distance between 2 consecutive points
 * exceeds the value of dist_gap, a data gap is assumed.  The user may
 * select a preferred direction where s/he wants the positive anomalies
 * to be pointing towards.  Positive normal vectors to the trackline
 * will then always be within +- 90 degrees of this direction.
 * Separate colors may be specified for positive anomaly, outline of
 * anomaly, and trackline.  Plotting of the outline and track is optional.
 *
 * Author:	Paul Wessel
 * Date:	20-JUN-2000
 * Version:	4
 *
 */
 
#include "gmt.h"

#define MAX_POINTS 400

double *xx, *yy, *zz, *lon, *lat, *z;
double start_az, stop_az;

void GMT_draw_z_scale (double x0, double y0, double length, double zscale, BOOLEAN gave_xy, char *units);
void plot_wiggle (double *x, double *y, double *z, int np, double zscale, BOOLEAN fixed, double fix_az, struct GMT_FILL *fill, struct GMT_PEN *pen_o, struct GMT_PEN *pen_t, BOOLEAN paint_wiggle, BOOLEAN negative, BOOLEAN outline, BOOLEAN track);

main (int argc, char **argv)
{
	int i, j, n, k, n_files = 0, fno, n_alloc = GMT_CHUNK, n_args, n_expected_fields, n_fields, wantx, wanty;
	
	BOOLEAN error = FALSE, nofile = TRUE, outline = FALSE, track = FALSE, paint_wiggle, draw_scale = FALSE;
	BOOLEAN do_fill = FALSE, negative = FALSE, fix = FALSE, done, gave_xy = FALSE;
	
	double west = 0.0, east = 360.0, south = 0.0, north = 0.0, x_2, y_2, bias = 0.0, fix_az;
	double zscale = 0.0, dx, dy, dz, ds, x0, y0, length, dist_gap = DBL_MAX, *in, pref_az = 0.0;
	
	char line[BUFSIZ], *units, txt_a[32], txt_b[32], unit = 0;
	
	FILE *fp = NULL;
	
	struct GMT_FILL fill;
	struct GMT_PEN pen_o, pen_t;
	
	argc = GMT_begin (argc, argv);
	
	GMT_init_pen (&pen_o, GMT_PENWIDTH);
	GMT_init_pen (&pen_t, GMT_PENWIDTH);
	GMT_init_fill (&fill, gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
	
	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 'A':
					pref_az = atof (&argv[i][2]);
					break;
				case 'D':
					dist_gap = atof (&argv[i][2]);
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':
					if (GMT_getfill (&argv[i][2], &fill)) {
						GMT_fill_syntax ('G');
						error++;
					}
					do_fill = TRUE;
					break;
				case 'I':
					fix_az = atof (&argv[i][2]);
					fix = TRUE;
					break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':
					negative = TRUE;
					break;
				case 'C':
					bias = atof (&argv[i][2]);
					break;
				case 'T':
					if (GMT_getpen (&argv[i][2], &pen_t)) {
						GMT_pen_syntax ('T');
						error++;
					}
					track = TRUE;
					break;
				case 'W':
					if (GMT_getpen (&argv[i][2], &pen_o)) {
						GMT_pen_syntax ('W');
						error++;
					}
					outline = TRUE;
					break;
				case 'S':
					j = 0;
					if (argv[i][2] == 'x') gave_xy = TRUE, j = 1;
					k = sscanf (&argv[i][2+j], "%[^/]/%[^/]/%lf", txt_a, txt_b, &length);
					wantx = (gave_xy) ? GMT_IS_FLOAT : GMT_IS_LON;
					wanty = (gave_xy) ? GMT_IS_FLOAT : GMT_IS_LAT;
					error += GMT_verify_expectations (wantx, GMT_scanf_arg (txt_a, wantx, &x0), txt_a);
					error += GMT_verify_expectations (wanty, GMT_scanf_arg (txt_b, wanty, &y0), txt_b);
					units = strrchr (argv[i], '/');	units++;
					if (k != 3) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Correct syntax\n", GMT_program);
						fprintf (stderr, "	-S[x]<x0>/<y0>/<length>[/<units>]\n");
						error++;
					}
					draw_scale = TRUE;
					break;
					
				case 'Z':
					j = strlen (argv[i]) - 1;
					if (strchr ("cimpCIMP", (int)argv[i][j])) unit = argv[i][j];
					zscale = atof (&argv[i][2]);
					break;
					
				/* Illegal options */
			
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"pswiggle %s - Plot xyz-series along tracks\n\n", GMT_VERSION);
		fprintf (stderr, "usage: pswiggle <xyz-files> -J<params> -R<west/east/south/north> -Z<scale>\n");
		fprintf (stderr, "\t[-A<azimuth>] [-B<tickinfo>] [-C<center>] [-D<gap>] [-E<az/el>] [-G<fill>] [-H[<nrec>]] [-I<az>]\n");
		fprintf (stderr, "\t[-K] [-M[<flag>]] [-N] [-O] [-P] [-S[x]<lon0>/<lat0>/<length>/<units>] [-T<trackpen>] -U] [-V]\n");
		fprintf (stderr, "\t[[-W<outlinepen>] [-X<x_shift>] [-Y<y_shift>] [-c<ncopies>] [-:] [-bi[s][<n>] [-f[i|o]<colinfo>]\n");

		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<xyz_files> is one or more files.  If none, read standard input\n");
		GMT_explain_option ('j');
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A set azimuth for preferred positive wiggle orientation [0.0 (north)]\n");
		GMT_explain_option ('b');
		fprintf (stderr, "\t-C sets center value to be removed from z before plotting. [0]\n");
		fprintf (stderr, "\t-D means there is a datagap if 2 points are more than <gap> distance units apart\n");
		fprintf (stderr, "\t   gap is in km for map projections, user-units otherwise\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-G Specify color or pattern for positive areas. 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");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I set fixed projection azimuths for wiggles\n");
		GMT_explain_option ('K');
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N Fill negative wiggles instead [Default is positive]\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-S draws a simple vertical scale centered on <lon0>/<lat0>.  Use -Sx to specify cartesian coordinates instead.\n");
		fprintf (stderr, "\t   <length> is in z-units, append unit name for labeling\n");
		fprintf (stderr, "\t-T specifies track pen attributes. [Default is no track]\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W specifies outline pen attributes. [Default is no outline]\n");
		GMT_explain_option ('X');
		fprintf (stderr, "\t-Z gives the wiggle scale in data-units per %s\n", GMT_unit_names[gmtdefs.measure_unit]);
		GMT_explain_option ('c');
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 3 input columns\n");
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (!(outline || do_fill)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -G, -W\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 (zscale == 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  scale must be nonzero\n", GMT_program);
		error++;
	}
	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] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = 3;
        if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < 3) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 3 columns\n", GMT_program);
		error++;
	}
	      
	if (error) exit (EXIT_FAILURE);

	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]]);
	}

	GMT_map_setup (west, east, south, north);

	if (MAPPING) dist_gap *= 1000.0;    /* Distance in meters */
	
	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);
	if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);

	zscale = 1.0 / zscale;

	switch (unit) {	/* Adjust for possible unit selection */
		case 'C':
		case 'c':
			zscale *= GMT_u2u[GMT_CM][GMT_INCH];
			break;
		case 'I':
		case 'i':
			zscale *= GMT_u2u[GMT_INCH][GMT_INCH];
			break;
		case 'M':
		case 'm':
			zscale *= GMT_u2u[GMT_M][GMT_INCH];
			break;
		case 'P':
		case 'p':
			zscale *= GMT_u2u[GMT_PT][GMT_INCH];
			break;
		default:
			zscale *= GMT_u2u[gmtdefs.measure_unit][GMT_INCH];
			break;
	}
	
	/* Now convert angles to radians and keep it that way */
	start_az = (pref_az - 90.0) * D2R;
	stop_az  = (pref_az + 90.0) * D2R;
	fix_az *= D2R;
	
	GMT_map_clip_on (GMT_no_rgb, 3);
	
	/* Allocate memory */
		
	lon = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	lat = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	z   = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	xx  = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	yy  = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	zz  = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	
	if (n_files != 0) nofile = FALSE;
	done = FALSE;
	n_args = (argc > 1) ? argc : 2;
	
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over all the files */
		if (!nofile && argv[fno][0] == '-') continue;
		if (nofile) {
			fp = GMT_stdin;
			done = TRUE;
			if (gmtdefs.verbose) fprintf (stderr, "%s: Reading data from standard input\n", GMT_program);
#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;
		}
		
		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 (GMT_io.multi_segments) {
			GMT_fgets (line, BUFSIZ, fp);
			if (gmtdefs.verbose) ps_comment (line);
		}
		n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : 3;

		n_fields = GMT_input (fp, &n_expected_fields, &in);
		while (! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

			while (GMT_io.status & GMT_IO_SEGMENT_HEADER) {
				if (gmtdefs.verbose) ps_comment (GMT_io.segment_header);
				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			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_fields, n);
					exit (EXIT_FAILURE);
				}

				lon[n] = in[0];
				lat[n] = in[1];
				z[n] = in[2] - bias;
				n++;
				if (n >= n_alloc) {
					n_alloc += GMT_CHUNK;
					lon = (double *) GMT_memory ((void *)lon, (size_t)n_alloc, sizeof (double), GMT_program);
					lat = (double *) GMT_memory ((void *)lat, (size_t)n_alloc, sizeof (double), GMT_program);
					z = (double *) GMT_memory ((void *)z, (size_t)n_alloc, sizeof (double), GMT_program);
					xx = (double *) GMT_memory ((void *)xx, (size_t)n_alloc, sizeof (double), GMT_program);
					yy = (double *) GMT_memory ((void *)yy, (size_t)n_alloc, sizeof (double), GMT_program);
					zz = (double *) GMT_memory ((void *)zz, (size_t)n_alloc, sizeof (double), GMT_program);
				}

				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
					
			GMT_geo_to_xy (lon[0], lat[0], &xx[0], &yy[0]);
			zz[0] = z[0];
			paint_wiggle = (do_fill && ((negative && z[0] <= 0.0) || (!negative && z[0] >= 0.0)));
			j = 1;
			for (i = 1; i < n; i++) {	/* Convert to inches/cm and get distance increments */
				dx = lon[i] - lon[i-1];
				dy = lat[i] - lat[i-1];
				if (MAPPING) {	/* Distance in meters */
					dx *= cos (0.5 * (lat[i] + lat[i-1]) * D2R);
					ds = hypot (dx, dy) * project_info.M_PR_DEG;
				}
				else
					ds = hypot (dx, dy);
					
				GMT_geo_to_xy (lon[i], lat[i], &x_2, &y_2);

				if (j > 0 && (ds > dist_gap || GMT_is_dnan (z[i]))) {	/* Data gap, plot what we have */
					paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0)));
					plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track);
					j = 0;
				}
				else if (j >= MAX_POINTS) {
					paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0)));
					plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track);
					xx[0] = xx[j-1];
					yy[0] = yy[j-1];
					zz[0] = zz[j-1];
					j = 1;
				}
				else if (!GMT_is_dnan (z[i-1]) && (z[i]*z[i-1] < 0.0 || z[i] == 0.0)) {	/* Crossed 0, add new point and plot */
					dz = z[i] - z[i-1];
					xx[j] = (dz == 0.0) ? xx[j-1] : xx[j-1] + fabs (z[i-1] / dz) * (x_2 - xx[j-1]);
					yy[j] = (dz == 0.0) ? yy[j-1] : yy[j-1] + fabs (z[i-1] / dz) * (y_2 - yy[j-1]);
					zz[j] = 0.0;
					j++;
					paint_wiggle = (do_fill && ((negative && zz[j-2] <= 0.0) || (!negative && zz[j-2] >= 0.0)));
					plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track);
					xx[0] = xx[j-1];
					yy[0] = yy[j-1];
					zz[0] = zz[j-1];
					j = 1;
				}
				xx[j] = x_2;
				yy[j] = y_2;
				zz[j] = z[i];
				if (!GMT_is_dnan (z[i])) j++;
			}
			
			if (j > 1) {
				paint_wiggle = (do_fill && ((negative && zz[j-1] <= 0.0) || (!negative && zz[j-1] >= 0.0)));
				plot_wiggle (xx, yy, zz, j, zscale, fix, fix_az, &fill, &pen_o, &pen_t, paint_wiggle, negative, outline, track);
			}
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
			
	GMT_map_clip_off();
	GMT_map_basemap ();

	if (draw_scale) GMT_draw_z_scale (x0, y0, length, zscale, gave_xy, units);
	
	if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
	ps_plotend (gmtdefs.last_page);
	
	GMT_free ((void *)lon);
	GMT_free ((void *)lat);
	GMT_free ((void *)z);
	GMT_free ((void *)xx);
	GMT_free ((void *)yy);
	GMT_free ((void *)zz);
	
	GMT_end (argc, argv);
}
	
void plot_wiggle (double *x, double *y, double *z, int np, double zscale, BOOLEAN fixed, double fix_az, struct GMT_FILL *fill, struct GMT_PEN *pen_o, struct GMT_PEN *pen_t, BOOLEAN paint_wiggle, BOOLEAN negative, BOOLEAN outline, BOOLEAN track)
{
	double dx, dy, len, az = 0.0, s, c, x_inc, y_inc;
	int n, i;
	static char *name[2] = {"Positive", "Negative"}, comment[32];

	if (fixed) {
		az = fix_az;
		sincos (az, &s, &c);
	}

	if (paint_wiggle || outline) {
		if (!GMT_n_alloc) GMT_get_plot_array ();
		GMT_x_plot[0] = x[0];
		GMT_y_plot[0] = y[0];
		n = 1;
		for (i = 0; i < np; i++) {
			if (!fixed && i < np-1) {
				dx = x[i+1] - x[i];
				dy = y[i+1] - y[i];
				if (!(dx == 0.0 && dy == 0.0)) az = -d_atan2 (dy, dx) - TWO_PI;	/* Azimuth of normal to track */
				while (az < start_az) az += TWO_PI;
				if (az > stop_az) az -= M_PI;
			}
			if (fabs (z[i]) > 0.0) {
				if (!fixed) sincos (az, &s, &c);
				len = zscale * z[i];
				x_inc = len * s;
				y_inc = len * c;
			}
			else
				x_inc = y_inc = 0.0;
				
			GMT_x_plot[n] = x[i] + x_inc;
			GMT_y_plot[n] = y[i] + y_inc;
			n++;
			if (n == GMT_n_alloc) GMT_get_plot_array ();
		}
		GMT_x_plot[n] = x[np-1];
		GMT_y_plot[n] = y[np-1];
		n++;
		
		if (paint_wiggle) {
			for (i = np - 2; i >= 0; i--, n++) {	/* Go back to 1st point along track */
				if (n == GMT_n_alloc) GMT_get_plot_array ();
				GMT_x_plot[n] = x[i];
				GMT_y_plot[n] = y[i];
			}
		}
		if (project_info.three_D) GMT_2D_to_3D (GMT_x_plot, GMT_y_plot, project_info.z_level, n);
	}
	
	
	if (paint_wiggle) { /* First shade wiggles */
		sprintf (comment, "%s wiggle", name[negative]);
		ps_comment (comment);
		GMT_fill (GMT_x_plot, GMT_y_plot, n, fill, FALSE);
	}
	
	if (outline) { /* Then draw wiggle outline */
		sprintf (comment, "Wiggle line");
		ps_comment (comment);
		GMT_setpen (pen_o);
		ps_line (&GMT_x_plot[1], &GMT_y_plot[1], np, 3, FALSE, TRUE);
	}
	
	if (track) {	/* Finally draw track line */
		if (project_info.three_D) GMT_2D_to_3D (x, y, project_info.z_level, np);
		sprintf (comment, "Track line");
		ps_comment (comment);
		GMT_setpen (pen_t);
		ps_line (x, y, np, 3, FALSE, TRUE);
	}
}

void GMT_draw_z_scale (double x0, double y0, double length, double zscale, BOOLEAN gave_xy, char *units)
{
	double dy, off, xx[4], yy[4];
	char txt[128];
	
	GMT_setpen (&gmtdefs.tick_pen);

	if (!gave_xy) {
		GMT_geo_to_xy (x0, y0, &xx[0], &yy[0]);
		x0 = xx[0];	y0 = yy[0];
	}
		
	if (units)
		sprintf (txt, "%g %s", length, units);
	else
		sprintf (txt, "%g", length);
	dy = 0.5 * length * zscale;	
	GMT_xyz_to_xy (x0 + gmtdefs.map_scale_height, y0 - dy, 0.0, &xx[0], &yy[0]);
	GMT_xyz_to_xy (x0, y0 - dy, 0.0, &xx[1], &yy[1]);
	GMT_xyz_to_xy (x0, y0 + dy, 0.0, &xx[2], &yy[2]);
	GMT_xyz_to_xy (x0 + gmtdefs.map_scale_height, y0 + dy, 0.0, &xx[3], &yy[3]);
	ps_line (xx, yy, 4, 3, FALSE, TRUE);
	off = ((gmtdefs.map_scale_height > 0.0) ? gmtdefs.tick_length : 0.0) + gmtdefs.annot_offset[0];
	GMT_text3D (x0 + off, y0, project_info.z_level, gmtdefs.annot_font_size[0], gmtdefs.annot_font[0], txt, 0.0, 5, 0);
}
