/*-------------------------------------------------------------------- * $Id: psrose.c,v 1.4.4.3 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * 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.1 */ #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\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 [-A[r]] [-B] \n"); fprintf (stderr, "\t[-C[]] [-Eaz/el] [-G] [-H[]] [-I] [-K] [-M]\n"); fprintf (stderr, "\t[-N] [-O] [-P] [-R] [-Sscale[n]] [-T] [-U[label]] [-V] [-:]\n"); fprintf (stderr, "\t[-Wpen] [-Xx_shift] [-Yy_shift] [-Zscale] [-cncopies] [-bi[s][]]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\t (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 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 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 -Sn\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, ""); sprintf (format, "%%s: Info for %%s: n = %%d rmax = %s mean r/az = (%s/%s) totlength = %s\n", 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 >= 100); /* 0 is small [Default], 1 is large */ sprintf (lw, "90%sW", GMT_degree_symbol[which]); sprintf (le, "90%sE", GMT_degree_symbol[which]); sprintf (ln, "0%s", 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); }