/*--------------------------------------------------------------------
 *	$Id: psrose.c,v 1.4 2001/04/11 23:07:29 pwessel Exp $
 *
 *	Copyright (c) 1991-2001 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
 *--------------------------------------------------------------------*/
/*
 * psrose reads a file [or standard input] with azimuth and length information
 * and draws a sector or rose diagram.  Several options for plot layout are available.
 * 2 diagrams are possible: Full circle (360) or half circle (180).  In the
 * latter case azimuths > 180 are reversed (-= 180).  psrose are fully
 * compatible with other gmtsystem v2.0 plot files and overlays.
 *
 * To be compatible with GMT, I assume radial distance to be "x"
 * and azimuth to be "y".  Hence, west = 0.0 and east = max_radius
 * south/north is -90,90 for halfcircle and 0,360 for full circle
 *
 * Author:	Paul Wessel
 * Date:	20-JUN-2000
 * Version:	3.4
 */
 
#include "gmt.h"

int n_alloc = GMT_CHUNK;

double *xx, *yy, *sum;
double *azimuth, *length;
double *mode_direction;
double *mode_length;

void plot3d(double x, double y, int pen);
void pie3d(double x, double y, double rad, double a1, double a2, int rgb[], BOOLEAN outline);
void arc3d(double x, double y, double r, double a1, double a2, int status);

main (int argc, char **argv)
{
	int i, n_bins, bin, n_anot, n = 0, a_rgb[3];
	int status, n_alpha, n_modes, n_fields, n_expected_fields;
	
	BOOLEAN error = FALSE, inquire = FALSE, outline = FALSE, find_mean = FALSE, normalize = FALSE;
	BOOLEAN half_only = FALSE, rose = FALSE, windrose = TRUE, sector_plot = FALSE;
	BOOLEAN eq_area = FALSE, draw_modes = FALSE, automatic = FALSE, oriented = FALSE;
	
	double scale = 3.0, max = 0.0, sector = 0.0, radius, az, x_origin, y_origin, tmp;
	double angle1, angle2, angle, x, y, mean_theta, mean_radius, xr = 0.0, yr = 0.0;
	double x1, x2, y1, y2, total = 0.0, total_arc, off, max_radius, az_offset, *in;
	double tailwidth = 0.03, headlength = 0.12, headwidth = 0.1, one_or_two = 1.0, s, c;
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, asize, lsize, zscale = 1.0;
	
	char text[BUFSIZ], lw[10], le[10], ln[10], *file = CNULL, format[BUFSIZ];
	char *choice[2] = {"OFF", "ON"}, txt_a[32], txt_b[32], txt_c[32];
	
	FILE *fp = NULL, *fpm = NULL;
	
	struct GMT_PEN pen;
	struct GMT_FILL fill;
	
	argc = GMT_begin (argc, argv);
	
	if (gmtdefs.measure_unit == GMT_CM) {
		scale = 7.5 / 2.54;
		tailwidth = 0.075 / 2.54;
		headlength = 0.3 / 2.54;
		headwidth = 0.25 / 2.54;
	}

	a_rgb[0] = a_rgb[1] = a_rgb[2] = 0;
	GMT_init_fill (&fill, -1, -1, -1);
	GMT_init_pen (&pen, GMT_PENWIDTH);
	asize = gmtdefs.anot_font_size * GMT_u2u[GMT_PT][GMT_INCH];
	lsize = gmtdefs.label_font_size * GMT_u2u[GMT_PT][GMT_INCH];
	strcpy (lw, "WEST");	strcpy (le, "EAST");	strcpy (ln, "NORTH");
	mode_direction = (double *) GMT_memory (VNULL, (size_t)1, sizeof (double), "psrose");
	mode_length = (double *) GMT_memory (VNULL, (size_t)1, sizeof (double), "psrose");
	
	/* 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 '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);
					west = 0.0;
					break;
					
				/* Supplemental options */
				
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'A':		/* Get Sector angle in degrees */
					sector = atof (&argv[i][2]);
					if (argv[i][strlen (argv[i])-1] == 'r') rose = TRUE;
					break;
				case 'C':		/* Read mode file and plot directions */
					draw_modes = TRUE;
					if ((fpm = fopen (&argv[i][2], "r")) == NULL) find_mean = TRUE;
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':		/* Set Gray shade */
					if (GMT_getfill (&argv[i][2], &fill)) {
						GMT_fill_syntax ('G');
						error++;
					}
					break;
				case 'I':		/* Compute statistics only - no plot */
					inquire = TRUE;
					break;
				case 'M':		/* Get arrow parameters */
					if (argv[i][2]) {
						n = sscanf (&argv[i][2], "%[^/]/%[^/]/%[^/]/%d/%d/%d", txt_a, txt_b, txt_c, &a_rgb[0], &a_rgb[1], &a_rgb[2]);
						if (n != 6) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -M option.  Correct syntax:\n", GMT_program);
							fprintf (stderr, "\t-M<tailwidth/headlength/headwidth/r/g/b>\n");
							error = TRUE;
						}
						else if (GMT_check_rgb (a_rgb)) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -M option.  RGB must all be in 0-255 range\n", GMT_program);
							error = TRUE;
						}
						else {
							tailwidth  = GMT_convert_units (txt_a, GMT_INCH);
							headlength = GMT_convert_units (txt_b, GMT_INCH);
							headwidth  = GMT_convert_units (txt_c, GMT_INCH);
						}
					}
					break;
				case 'N':		/* Make sectors area be proportional to frequency instead of radius */
					eq_area = TRUE;
					break;
				case 'S':		/* Get radius of unit circle in inches */
					n = strlen (argv[i]) - 1;
					if (argv[i][n] == 'n') {
						normalize = TRUE;
						argv[i][n] = 0;
					}
					scale = GMT_convert_units (&argv[i][2], GMT_INCH);
					break;
				case 'T':		/* Oriented instead of directed data */
					oriented = TRUE;
					break;
				case 'W':		/* Get pen width for outline */
					if (GMT_getpen (&argv[i][2], &pen)) {
						GMT_pen_syntax ('W');
						error++;
					}
					outline = TRUE;
					break;
				case 'Z':		/* Scale radii before using data */
					zscale = atof (&argv[i][2]);
					break;
				default:		/* Options not recognized */
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			file = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"psrose %s - Polar histogram (rose diagram) plotter\n\n", GMT_VERSION);
		fprintf (stderr, "usage: psrose <infile> [-A<sector_angle>[r]] [-B<tickinfo>] \n");
		fprintf (stderr, "\t[-C[<modes>]] [-Eaz/el] [-G<fill>] [-H[<nrec>]] [-I] [-K] [-M<parameters>]\n");
		fprintf (stderr, "\t[-N] [-O] [-P] [-R<r0/r1/theta0/theta1>] [-Sscale[n]] [-T] [-U[label]] [-V] [-:]\n");
		fprintf (stderr, "\t[-Wpen] [-Xx_shift] [-Yy_shift] [-Zscale] [-cncopies] [-bi[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<infile> (in ASCII or binary) has (length,azimuth) pairs.  If not given read standard input\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A sector width in degrees for sector diagram [Default is windrose]\n");
		fprintf (stderr, "\t   append r to get rose diagram\n");
		GMT_explain_option ('B');
		fprintf (stderr, "\t   (Remember: radial is x-direction, azimuthal is y-direction)\n");
		fprintf (stderr, "\t-C plot vectors listed in the <modes> file.  If no file, use mean direction\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-G <fill> specify color for diagram.  [Default is no fill]\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I for inquire.  Only compute statistics - no plot is created\n");
		GMT_explain_option ('K');
		fprintf (stderr, "\t-M Append tailwidth/headlength/headwidth/r/g/b to set arrow attributes [0.03i/0.12i/0.1i/0/0/0\n");
		fprintf (stderr, "\t-N normalizes rose plots for area, i.e. takes sqrt(r) before plotting [FALSE]\n");
		fprintf (stderr, "\t   Only applicable if normalization has been specified with -S<radius>n\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-R specifies the region.  (r0 = 0, r1 = max_radius.  For azimuth:\n");
		fprintf (stderr, "\t   Specify theta0/theta1 = -90/90 (half-circle) or 0/360 only)\n");
		fprintf (stderr, "\t   If r0 = r1 = 0, psrose will compute reasonable a r1 value.\n");
		fprintf (stderr, "\t-S specifies the radius of the unit circle in %s [%lg]. Normalize r if n is appended\n", GMT_unit_names[gmtdefs.measure_unit], scale);
		fprintf (stderr, "\t-T indicates that the vectors are oriented (two-headed), not directed [Default]\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets pen attributes for outline of rose. [Default is no outline]\n");
		GMT_explain_option ('X');
		fprintf (stderr, "\t-Z multiply the radii by scale before plotting\n");
		GMT_explain_option ('c');
		fprintf (stderr, "\t-: Expect (azimuth,radius) input rather than (radius,azimuth) [%s]\n", choice[gmtdefs.xy_toggle]);
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns.\n");
		fprintf (stderr, "\t(See man page for gmtdefaults to set these and other defaults to ON or OFF)\n");
		exit (EXIT_FAILURE);
	}

	/* Check that the options selected are mutually consistant */
	
	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 (scale <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  radius must be nonzero\n", GMT_program);
		error++;
	}
	if (fabs (zscale) < GMT_CONV_LIMIT) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  factor must be nonzero\n", GMT_program);
		error++;
	}
	if (sector < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -A option:  sector width must be positive\n", GMT_program);
		error++;
	}
	if (!((south == -90.0 && north == 90.0) || (south == 0.0 && north == 360.0))) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -R option:  theta0/theta1 must be either -90/90 or 0/360\n", GMT_program);
		error++;
	}
	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] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < 2) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 2 columns\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

	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 (file && (fp = GMT_fopen (file, GMT_io.r_mode)) == NULL) {
		fprintf (stderr, "%s:  Cannot open file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	if (fp == NULL) {
		fp = GMT_stdin;
#ifdef SET_IO_MODE
		GMT_setmode (0);
#endif
	}
	
	max_radius = east;
	half_only = (fabs (south + 90.0) < GMT_CONV_LIMIT);
	if (rose) windrose = FALSE;
	sector_plot = (sector > 0.0);
	if (sector_plot) windrose = FALSE;	/* Draw rose diagram instead of sector diagram */
	if (!normalize) eq_area = FALSE;	/* Only do this is data is normalized for length also */
	if (!inquire && !project_info.region_supplied) automatic = TRUE;
	if (oriented) one_or_two = 2.0;
	
	if (half_only) {
		total_arc = 180.0;
		az_offset = 90.0;
	}
	else {
		total_arc = 360.0;
		az_offset = 0.0;
	}
	n_bins = (sector <= 0.0) ? 1 : irint (total_arc / sector);
	
	sum = (double *) GMT_memory (VNULL, (size_t)n_bins, sizeof (double), GMT_program);
	xx = (double *) GMT_memory (VNULL, (size_t)(n_bins+2), sizeof (double), GMT_program);
	yy = (double *) GMT_memory (VNULL, (size_t)(n_bins+2), sizeof (double), GMT_program);
	azimuth = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	length = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	
	if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (text, BUFSIZ, fp);
		
	/* Read data and do some stats */
	
	n = 0;
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2;
	
	n_fields = GMT_input (fp,  &n_expected_fields, &in);

	while (! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

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

		length[n] = in[0];
		azimuth[n] = in[1];
		
		if (zscale != 1.0) length[n] *= zscale;
		
		/* Make sure azimuth is in 0 <= az < 360 range */
		
		while (azimuth[n] < 0.0) azimuth[n] += 360.0;
		while (azimuth[n] >= 360.0) azimuth[n] -= 360.0;
		
		if (half_only) {	/* Flip azimuths about E-W line i.e. -90 < az <= 90 */
			if (azimuth[n] > 90.0 && azimuth[n] <= 270.0) azimuth[n] -= 180.0;
			if (azimuth[n] > 270.0) azimuth[n] -= 360.0;
		}
		
		/* Double angle to find mean azimuth */
		
		sincos (one_or_two * azimuth[n] * D2R, &s, &c);
		xr += length[n] * c;
		yr += length[n] * s;
		
		total += length[n];
		
		n++;
		if (n == n_alloc) {	/* Get more memory */
			n_alloc += GMT_CHUNK;
			azimuth = (double *) GMT_memory ((void *)azimuth, (size_t)n_alloc, sizeof (double), GMT_program);
			length = (double *) GMT_memory ((void *)length, (size_t)n_alloc, sizeof (double), GMT_program);
		}
		n_fields = GMT_input (fp, &n_expected_fields, &in);
	}
	if (fp != GMT_stdin) GMT_fclose (fp);
	
	if (sector > 0.0) {	/* Sum up sector diagram info */
		for (i = 0; i < n; i++) {
			bin = (int) ((azimuth[i] + az_offset) / sector);
			sum[bin] += length[i];
		}
	}
		
	mean_theta = R2D * d_atan2 (yr, xr) / one_or_two;
	if (mean_theta < 0.0) mean_theta += 360.0;
	mean_radius = hypot (xr, yr) / total;
	if (!normalize) mean_radius *= max_radius;
	
	if (sector > 0.0) {	/* Find max of the bins */
		for (bin = 0; bin < n_bins; bin++) if (sum[bin] > max) max = sum[bin];
		if (normalize) for (bin = 0; bin < n_bins; bin++) sum[bin] /= max;
		if (eq_area) for (bin = 0; bin < n_bins; bin++) sum[bin] = sqrt (sum[bin]);
	}
	else {
		for (i = 0; i < n; i++) if (length[i] > max) max = length[i];
		if (normalize) {
			max = 1.0 / max;
			for (i = 0; i < n; i++) length[i] *= max;
			max = 1.0 / max;
		}
	}
	
	if (inquire || gmtdefs.verbose) {
		if (file) strcpy (text, file); else strcpy (text, "<stdin>");
		sprintf (format, "%%s: Info for %%s: n = %%d rmax = %s mean r/az = (%s/%s) totlength = %s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, text, n, max, mean_radius, mean_theta, total);
		if (inquire) {
			GMT_free ((void *)sum);
			GMT_free ((void *)xx);
			GMT_free ((void *)yy);
			GMT_free ((void *)azimuth);
			GMT_free ((void *)length);
			exit (EXIT_FAILURE);
		}
	}
	
	if (automatic) {
		if (fabs (frame_info.anot_int[0]) < GMT_CONV_LIMIT) {
			tmp = pow (10.0, floor (d_log10 (max)));
			if ((max / tmp) < 3.0) tmp *= 0.5;
		}
		else
			tmp = frame_info.anot_int[0];
		max_radius = ceil (max / tmp) * tmp;
		if (fabs (frame_info.anot_int[0]) < GMT_CONV_LIMIT || fabs (frame_info.grid_int[0]) < GMT_CONV_LIMIT) {	/* Tickmarks not set */
			frame_info.anot_int[0] = frame_info.grid_int[0] = tmp;
			frame_info.plot = TRUE;
		}
	}
	
	if (frame_info.plot && fabs (frame_info.anot_int[1]) < GMT_CONV_LIMIT) frame_info.anot_int[1] = frame_info.grid_int[1] = 30.0;
	
	/* Ready to plot */
	
	project_info.projection = LINEAR;
	project_info.region = TRUE;
	project_info.pars[0] = project_info.pars[1] = (normalize) ? 1.0 : scale / max_radius;
	GMT_map_setup (-scale, scale, -scale, scale);
	
	ps_plotinit (NULL, 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);
	
	x_origin = scale;	y_origin = ((half_only) ? 0.0 : scale);
	ps_transrotate (x_origin, y_origin, 0.0);
	if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin - y_origin, 0.0);
	if (!normalize) scale /= max_radius;
	
	GMT_setpen (&pen);
	if (windrose) {
		for (i = 0; i < n; i++) {
			angle = (90.0 - azimuth[i]) * D2R;
			radius = length[i] * scale;
			x = radius * cos (angle);
			y = radius * sin (angle);
			plot3d (0.0, 0.0, 3);
			plot3d (x, y, 2);
		}
	}
	
	if (sector_plot && !rose) {	/* Draw pie slices for sector plot */
			
		for (bin = 0; bin < n_bins; bin++) {
			az = bin * sector - az_offset;
			angle1 = (90.0 - az - sector);
			angle2 = angle1 + sector;
			pie3d (0.0, 0.0, sum[bin] * scale, angle1, angle2, fill.rgb, FALSE);
		}
	}
	else if (rose) {	/* Draw rose diagram */
	
		for (bin = i = 0; bin < n_bins; bin++, i++) {
			az = bin * sector - az_offset + 0.5 * sector;
			angle = 90.0 - az;
			xx[i] = scale * sum[bin] * cos (D2R * angle);
			yy[i] = scale * sum[bin] * sin (D2R * angle);
		}
		if (half_only) {
			xx[i] = scale * 0.5 * (sum[0] + sum[n_bins-1]);
			yy[i] = 0.0;
			i++;
			xx[i] = -xx[i-1];
			yy[i] = 0.0;
			i++;
		}
		if (project_info.three_D) GMT_2D_to_3D (xx, yy, i);
		ps_polygon (xx, yy, i, fill.rgb, outline);
	}
	
	if (sector_plot && outline && !rose) {	/* Draw a line outlining the pie slices */
		angle1 = (half_only) ? 180.0 : 90.0;
		angle2 = (half_only) ? 0.0 : 90.0;
		angle1 *= D2R;
		angle2 *= D2R;
		sincos (angle1, &s, &c);
		x1 = (sum[0] * scale) * c;
		y1 = (sum[0] * scale) * s;
		sincos (angle2, &s, &c);
		x2 = (sum[n_bins-1] * scale) * c;
		y2 = (sum[n_bins-1] * scale) * s;
		plot3d (x1, y1, 3);
		plot3d (x2, y2, 2);
		for (bin = n_bins-1; bin >= 0; bin--) {
			status = (bin == 0) ? 2 : 0;
			az = bin * sector - az_offset;
			angle1 = 90.0 - az - sector;
			angle2 = angle1 + sector;
			arc3d (0.0, 0.0, sum[bin] * scale, angle1, angle2, status);
			angle1 *= D2R;
		}
	}
	
	if (draw_modes) {
		
		if (find_mean) {	/* Use the mean direction only */
			n_modes = 1;
			mode_direction[0] = mean_theta;
			mode_length[0] = mean_radius;
		}
		else {	/* Get mode parameters from separate file */
			n_modes = 0;
			n_alloc = 1;
			while (fgets (text, BUFSIZ, fpm)) {
				sscanf (text, "%lf %lf", &mode_direction[n_modes], &mode_length[n_modes]);
				n_modes++;
				if (n_modes == n_alloc) {
					n_alloc += GMT_SMALL_CHUNK;
					mode_direction = (double *) GMT_memory ((void *)mode_direction, (size_t)n_alloc, sizeof (double), "psrose");
					mode_length = (double *) GMT_memory ((void *)mode_length, (size_t)n_alloc, sizeof (double), "psrose");
				}
			}
			fclose (fpm);
			mode_direction = (double *) GMT_memory ((void *)mode_direction, (size_t)n_modes, sizeof (double), "psrose");
			mode_length = (double *) GMT_memory ((void *)mode_length, (size_t)n_modes, sizeof (double), "psrose");
		}
		for (i = 0; i < n_modes; i++) {
			if (eq_area) mode_length[i] = sqrt (mode_length[i]);
			if (half_only && mode_direction[i] > 90.0 && mode_direction[i] <= 270.0) mode_direction[i] -= 180.0;
			angle = D2R * (90.0 - mode_direction[i]);
			sincos (angle, &s, &c);
			xr = scale * mode_length[i] * c;
			yr = scale * mode_length[i] * s;
			GMT_vector3d (0.0, 0.0, xr, yr, project_info.z_level, tailwidth, headlength, headwidth, gmtdefs.vector_shape, a_rgb, TRUE);
		}
		
	}
	
	if (frame_info.plot) {	/* Draw grid lines etc */
		GMT_setpen (&gmtdefs.grid_pen);
		off = max_radius * scale;
		n_alpha = (frame_info.grid_int[1] > 0.0) ? irint (total_arc / frame_info.grid_int[1]) : -1;
		for (i = 0; i <= n_alpha; i++) {
			angle = D2R * i * frame_info.grid_int[1];
			sincos (angle, &s, &c);
			x = max_radius * scale * c;
			y = max_radius * scale * s;
			plot3d (0.0, 0.0, 3);
			plot3d (x, y, 2);
		}
		
		n_bins = (frame_info.grid_int[0] > 0.0) ? irint (max_radius / frame_info.grid_int[0]) : -1;
		for (i = 1; i <= n_bins; i++)
			arc3d (0.0, 0.0, i * frame_info.grid_int[0] * scale, 0.0, total_arc, 3);
		ps_setpaint (gmtdefs.basemap_frame_rgb);
		y = lsize + 6.0 * gmtdefs.anot_offset;
		GMT_text3d (0.0, off + y, project_info.z_level, gmtdefs.header_font_size, gmtdefs.header_font, frame_info.header, 0.0, 2, 0);
		
		GMT_get_format (frame_info.anot_int[0], frame_info.unit[0], format);
		
		if (half_only) {
			int which;
			which = (gmtdefs.degree_format >= 10);	/* 0 is small [Default], 1 is large */
			sprintf (lw, "90%sW\0", GMT_degree_symbol[which]);
			sprintf (le, "90%sE\0", GMT_degree_symbol[which]);
			sprintf (ln, "0%s\0", GMT_degree_symbol[which]);
			y = -(3.0 * gmtdefs.anot_offset + GMT_font_height[gmtdefs.anot_font] * asize);
			if (frame_info.label[0][0]) GMT_text3d (0.0, y, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, frame_info.label[0], 0.0, 10, 0);
			y = -(5.0 * gmtdefs.anot_offset + GMT_font_height[gmtdefs.anot_font] * lsize + GMT_font_height[gmtdefs.label_font] * lsize);
			if (frame_info.label[1][0]) GMT_text3d (0.0, y, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, frame_info.label[1], 0.0, 10, 0);
			GMT_text3d (0.0, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, "0", 0.0, 10, 0);
			n_anot = (frame_info.anot_int[0] > 0.0) ? irint (max_radius / frame_info.anot_int[0]) : -1;
			for (i = 1; i <= n_anot; i++) {
				x = i * frame_info.anot_int[0];
				sprintf (text, format, x);
				x *= scale;
				GMT_text3d (x, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, text, 0.0, 10, 0);
				GMT_text3d (-x, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, text, 0.0, 10, 0);
			}
		}
		else {
			GMT_text3d (0.0, -off - 2.0 * gmtdefs.anot_offset, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, "SOUTH", 0.0, 10, 0);
			plot3d (off, -off, 3);
			plot3d ((max_radius - frame_info.grid_int[0]) * scale, -off, 2);
			plot3d (off, -off, 3);
			plot3d (off, gmtdefs.tick_length - off, 2);
			plot3d ((max_radius - frame_info.grid_int[0]) * scale, -off, 3);
			plot3d ((max_radius - frame_info.grid_int[0]) * scale, gmtdefs.tick_length - off, 2);
			y = -(off + 5.0 * gmtdefs.anot_offset + GMT_font_height[gmtdefs.anot_font] * lsize + GMT_font_height[gmtdefs.label_font] * lsize);
			if (frame_info.label[1][0]) GMT_text3d (0.0, y, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, frame_info.label[1], 0.0, 10, 0);
			if (frame_info.label[0][0]) {
				strcat (format, " %s");
				sprintf (text, format, frame_info.grid_int[0], frame_info.label[0]);
			}
			else
				sprintf (text, format, frame_info.grid_int[0]);
			GMT_text3d ((max_radius - 0.5 * frame_info.grid_int[0]) * scale, -(off + gmtdefs.anot_offset), project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, text, 0.0, 10, 0);
		}
		GMT_text3d (off + 2.0 * gmtdefs.anot_offset, 0.0, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, le, 0.0, 5, 0);
		GMT_text3d (-off - 2.0 * gmtdefs.anot_offset, 0.0, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, lw, 0.0, 7, 0);
		GMT_text3d (0.0, off + 2.0 * gmtdefs.anot_offset, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, ln, 0.0, 2, 0);
		ps_setpaint (gmtdefs.background_rgb);
	}
	if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin + y_origin, 0.0);
	ps_rotatetrans (-x_origin, -y_origin, 0.0);

	ps_plotend (gmtdefs.last_page);
	
	GMT_end (argc, argv);
	
	GMT_free ((void *)sum);
	GMT_free ((void *)xx);
	GMT_free ((void *)yy);
	GMT_free ((void *)azimuth);
	GMT_free ((void *)length);
	
	exit (EXIT_SUCCESS);
}

void arc3d (double x, double y, double r, double a1, double a2, int status)
{
	int i, n;
	double da, a, *xx, *yy, s, c;
	
	if (project_info.three_D) {	/* Must project and draw line segments */
		n = (int)ceil (a2 - a1);
		
		xx = (double *) GMT_memory (VNULL, (size_t)(n + 1), sizeof (double), GMT_program);
		yy = (double *) GMT_memory (VNULL, (size_t)(n + 1), sizeof (double), GMT_program);
		
		a1 *= D2R;
		a2 *= D2R;
		da = (a2 - a1) / n;
		for (i = 0; i <= n; i++) {
			a = (i == n) ? a2 : a1 + da * i;
			sincos (a, &s, &c);
			xx[i] = x + r * c;
			yy[i] = y + r * s;
		}
		n++;	
		GMT_2D_to_3D (xx, yy, n);
		ps_line (xx, yy, n, status, FALSE, TRUE);
		GMT_free ((void *)xx);
		GMT_free ((void *)yy);
	}
	else
		ps_arc (x, y, r, a1, a2, status);
}

void plot3d (double x, double y, int pen)
{
	if (project_info.three_D) GMT_xyz_to_xy (x, y, project_info.z_level, &x, &y);
	ps_plot (x, y, pen);
}

void pie3d (double x, double y, double rad, double a1, double a2, int rgb[], BOOLEAN outline)
{
	int i, n;
	double da, a, *xx, *yy, s, c;
	
	if (project_info.three_D) {	/* Must project and draw line segments */
		n = (int)ceil (a2 - a1);
		
		xx = (double *) GMT_memory (VNULL, (size_t)(n + 2), sizeof (double), GMT_program);
		yy = (double *) GMT_memory (VNULL, (size_t)(n + 2), sizeof (double), GMT_program);
		
		a1 *= D2R;
		a2 *= D2R;
		da = (a2 - a1) / n;
		for (i = 0; i <= n; i++) {
			a = (i == n) ? a2 : a1 + da * i;
			sincos (a, &s, &c);
			xx[i] = x + rad * c;
			yy[i] = y + rad * s;
		}
		xx[i] = x;	yy[i] = y;
		n += 2;	
		GMT_2D_to_3D (xx, yy, n);
		ps_polygon (xx, yy, n, rgb, outline);
		GMT_free ((void *)xx);
		GMT_free ((void *)yy);
	}
	else
		ps_pie (x, y, rad, a1, a2, rgb, outline);
}