home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / psrose.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  21.9 KB  |  603 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)psrose.c    2.19  3/13/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * psrose reads a file [or standard input] with azimuth and length information
  9.  * and draws a sector or rose diagram.  Several options for plot layout are available.
  10.  * 2 diagrams are possible: Full circle (360) or half circle (180).  In the
  11.  * latter case azimuths > 180 are reversed (-= 180).  psrose are fully
  12.  * compatible with other gmtsystem v2.0 plot files and overlays.
  13.  *
  14.  * To be compatible with GMT, I assume radial distance to be "x"
  15.  * and azimuth to be "y".  Hence, west = 0.0 and east = max_radius
  16.  * south/north is -90,90 for halfcircle and 0,360 for full circle
  17.  *
  18.  * Author:    Paul Wessel
  19.  * Date:    20-JAN-1995
  20.  * Version:    3.0
  21.  */
  22.  
  23. #include "gmt.h"
  24.  
  25. int n_alloc = GMT_CHUNK;
  26.  
  27. double *xx, *yy, *sum;
  28. double *azimuth, *length;
  29. double mode_direction[50];    /* Max 50 modes allowed */
  30. double mode_length[50];
  31.  
  32. main (argc, argv)
  33. int argc;
  34. char **argv; {
  35.     int i, n_bins, bin, n_anot, n = 0, ar = 0, ag = 0, ab = 0;
  36.     int status, n_alpha, n_modes, ix, iy;
  37.     
  38.     BOOLEAN error = FALSE, inquire = FALSE, outline = FALSE, find_mean = FALSE, normalize = FALSE;
  39.     BOOLEAN half_only = FALSE, rose = FALSE, windrose = TRUE, sector_plot = FALSE;
  40.     BOOLEAN eq_area = FALSE, draw_modes = FALSE, automatic = FALSE;
  41.     
  42.     double scale = 3.0, max = 0.0, sector = 0.0, radius, az, x_origin, y_origin, tmp;
  43.     double angle1, angle2, angle, x, y, mean_theta, mean_radius, xr = 0.0, yr = 0.0;
  44.     double x1, x2, y1, y2, total = 0.0, total_arc, off, max_radius, az_offset;
  45.     double tailwidth = 0.02, headlength = 0.2, headwidth = 0.08, xy[2];
  46.     double west = 0.0, east = 0.0, south = 0.0, north = 0.0, asize, lsize, zscale = 1.0;
  47.     
  48.     char text[512], lw[10], le[10], ln[10], format[20], file[80];
  49.     
  50.     FILE *fp = NULL, *fpm = NULL;
  51.     
  52.     struct PEN pen;
  53.     struct FILL fill;
  54.     
  55.     argc = gmt_begin (argc, argv);
  56.     
  57.     if (gmtdefs.measure_unit) {
  58.         scale = 7.5;
  59.     }
  60.     gmt_init_fill (&fill, 0, 0, 0);
  61.     gmt_init_pen (&pen, 1);
  62.     asize = gmtdefs.anot_font_size / gmt_ppu[gmtdefs.measure_unit];
  63.     lsize = gmtdefs.label_font_size / gmt_ppu[gmtdefs.measure_unit];
  64.     strcpy (lw, "WEST");    strcpy (le, "EAST");    strcpy (ln, "NORTH");
  65.     strcpy (file, "(stdin)");
  66.     
  67.     /* Check and interpret the command line arguments */
  68.     
  69.     for (i = 1; i < argc; i++) {
  70.         if (argv[i][0] == '-') {
  71.             switch(argv[i][1]) {
  72.         
  73.                 /* Common parameters */
  74.             
  75.                 case 'B':
  76.                 case 'H':
  77.                 case 'K':
  78.                 case 'O':
  79.                 case 'P':
  80.                 case 'R':
  81.                 case 'U':
  82.                 case 'V':
  83.                 case 'X':
  84.                 case 'x':
  85.                 case 'Y':
  86.                 case 'y':
  87.                 case 'c':
  88.                 case ':':
  89.                 case '\0':
  90.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  91.                     west = 0.0;
  92.                     break;
  93.                     
  94.                 /* Supplemental options */
  95.                 
  96.                 case 'A':        /* Get Sector angle in degrees */
  97.                     sector = atof (&argv[i][2]);
  98.                     if (argv[i][strlen (argv[i])-1] == 'r') rose = TRUE;
  99.                     break;
  100.                 case 'C':        /* Read mode file and plot directions */
  101.                     draw_modes = TRUE;
  102.                     if ((fpm = fopen (&argv[i][2], "r")) == NULL) find_mean = TRUE;
  103.                     break;
  104.                 case 'E':
  105.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  106.                     break;
  107.                 case 'F':
  108.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  109.                         gmt_pen_syntax ('F');
  110.                         error++;
  111.                     }
  112.                     break;
  113.                 case 'G':        /* Set Gray shade */
  114.                     if (gmt_getfill (&argv[i][2], &fill)) {
  115.                         gmt_fill_syntax ('G');
  116.                         error++;
  117.                     }
  118.                     break;
  119.                 case 'I':        /* Compute statistics only - no plot */
  120.                     inquire = TRUE;
  121.                     break;
  122.                 case 'M':        /* Get arrow parameters */
  123.                     if (argv[i][2]) {
  124.                         n = sscanf (&argv[i][2], "%lf/%lf/%lf/%d/%d/%d", &tailwidth, &headlength, &headwidth, &ar, &ag, &ab);
  125.                         if (n != 6) {
  126.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -M option.  Correct syntax:\n", gmt_program);
  127.                             fprintf (stderr, "\t-M<tailwidth/headlength/headwidth/r/g/b>\n");
  128.                             error = TRUE;
  129.                         }
  130.                         else if (check_rgb (ar, ag, ab)) {
  131.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -M option.  RGB must all be in 0-255 range\n", gmt_program);
  132.                             error = TRUE;
  133.                         }
  134.                     }
  135.                     break;
  136.                 case 'N':        /* Make sectors area be proportional to frequency instead of radius */
  137.                     eq_area = TRUE;
  138.                     break;
  139.                 case 'S':        /* Get radius of unit circle in inches */
  140.                     scale = atof (&argv[i][2]);
  141.                     if (argv[i][strlen (argv[i])-1] == 'n') normalize = TRUE;
  142.                     break;
  143.                 case 'W':        /* Get pen width for outline */
  144.                     if (gmt_getpen (&argv[i][2], &pen)) {
  145.                         gmt_pen_syntax ('W');
  146.                         error++;
  147.                     }
  148.                     outline = TRUE;
  149.                     break;
  150.                 case 'Z':        /* Scale radii before using data */
  151.                     zscale = atof (&argv[i][2]);
  152.                     break;
  153.                 default:        /* Options not recognized */
  154.                     error = TRUE;
  155.                     gmt_default_error (argv[i][1]);
  156.                     break;
  157.             }
  158.         }
  159.         else {
  160.             strcpy (file, argv[i]);
  161.             if ((fp = fopen (file, "r")) == NULL) {
  162.                 fprintf (stderr, "psrose:  Cannot open file %s\n", file);
  163.                 exit (-1);
  164.             }
  165.         }
  166.     }
  167.     
  168.     if (argc == 1 || gmt_quick) {
  169.         fprintf (stderr,"psrose %s - Polar histogram (rose diagram) plotter\n\n", GMT_VERSION);
  170.         fprintf (stderr, "usage: psrose <infile> [-A<sector_angle>[r]] [-B<tickinfo>] \n");
  171.         fprintf (stderr, "    [-C[<modes>]] [-Eaz/el] [-F<r/g/b>] [-G<fill>] [-H] [-I] [-K] [-M<parameters>]\n");
  172.         fprintf (stderr, "    [-N] [-O] [-P] [-R<r0/r1/theta0/theta1>] [-Sscale[n]] [-U[label]] [-V] [-:]\n");
  173.         fprintf (stderr, "    [-Wpen] [-Xx_shift] [-Yy_shift] [-Zscale] [-cncopies]\n\n");
  174.         
  175.         if (gmt_quick) exit (-1);
  176.         
  177.         fprintf (stderr, "    <infile> has (length,azimuth) pairs.  If not given read standard input\n");
  178.         fprintf (stderr, "\n\tOPTIONS:\n");
  179.         fprintf (stderr, "    -A sector width in degrees for sector diagram [Default is windrose]\n");
  180.         fprintf (stderr, "       append r to get rose diagram\n");
  181.         explain_option ('B');
  182.         fprintf (stderr, "       (Remember: radial is x-direction, azimuthal is y-direction)\n");
  183.         fprintf (stderr, "    -C plot vectors listed in the <modes> file.  If no file, use mean direction\n");
  184.         fprintf (stderr, "    -D for inquire.  Only compute statistics - no plot is created\n");
  185.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  186.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  187.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  188.         fprintf (stderr, "    -G <fill> specify color for diagram.  Default is 0/0/0 (black)\n");
  189.         explain_option ('H');
  190.         explain_option ('K');
  191.         fprintf (stderr, "    -M Append tailwidth/headlength/headwidth/r/g/b to set arrow attributes [0.02/0.2/0.8/0/0/0\n");
  192.         fprintf (stderr, "    -N normalizes rose plots for area, i.e. takes sqrt(r) before plotting [FALSE]\n");
  193.         fprintf (stderr, "       Only applicable if normalization has been specified with -S<radius>n\n");
  194.         explain_option ('O');
  195.         explain_option ('P');
  196.         fprintf (stderr, "    -R specifies the region.  (r0 = 0, r1 = max_radius.  For azimuth:\n");
  197.         fprintf (stderr, "       Specify theta0/theta1 = -90/90 (half-circle) or 0/360 only)\n");
  198.         fprintf (stderr, "       If r0 = r1 = 0, psrose will compute reasonable a r1 value.\n");
  199.         fprintf (stderr, "    -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);
  200.         explain_option ('U');
  201.         explain_option ('V');
  202.         fprintf (stderr, "    -W sets pen attributes for outline of rose. [Default is no outline]\n");
  203.         explain_option ('X');
  204.         fprintf (stderr, "    -Z multiply the radii by scale before plotting\n");
  205.         explain_option ('c');
  206.         fprintf (stderr, "    -: Expect (azimuth,radius) input rather than (radius,azimuth) [%s]\n", gmt_choice[gmtdefs.xy_toggle]);
  207.         fprintf (stderr, "    (See man page for gmtdefaults to set these and other defaults to ON or OFF)\n");
  208.         exit (-1);
  209.     }
  210.  
  211.     /* Check that the options selected are mutually consistant */
  212.     
  213.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  214.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  215.         error++;
  216.     }
  217.     if (scale <= 0.0) {
  218.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  radius must be nonzero\n", gmt_program);
  219.         error++;
  220.     }
  221.     if (zscale == 0.0) {
  222.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  factor must be nonzero\n", gmt_program);
  223.         error++;
  224.     }
  225.     if (sector < 0.0) {
  226.         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option:  sector width must be positive\n", gmt_program);
  227.         error++;
  228.     }
  229.     if (!((south == -90.0 && north == 90.0) || (south == 0.0 && north == 360.0))) {
  230.         fprintf (stderr, "%s: GMT SYNTAX ERROR -R option:  theta0/theta1 must be either -90/90 or 0/360\n", gmt_program);
  231.         error++;
  232.     }
  233.     
  234.     if (error) exit (-1);
  235.  
  236.     if (fp == NULL) fp = stdin;
  237.     max_radius = east;
  238.     half_only = (south == -90.0);
  239.     if (rose) windrose = FALSE;
  240.     sector_plot = (sector > 0.0);
  241.     if (sector_plot) windrose = FALSE;    /* Draw rose diagram instead of sector diagram */
  242.     if (!normalize) eq_area = FALSE;    /* Only do this is data is normalized for length also */
  243.     if (!inquire && west == east) automatic = TRUE;
  244.     
  245.     if (half_only) {
  246.         total_arc = 180.0;
  247.         az_offset = 90.0;
  248.     }
  249.     else {
  250.         total_arc = 360.0;
  251.         az_offset = 0.0;
  252.     }
  253.     n_bins = (sector <= 0.0) ? 1 : total_arc / sector;
  254.     
  255.     sum = (double *) memory (CNULL, n_bins, sizeof (double), "psrose");
  256.     xx = (double *) memory (CNULL, n_bins+2, sizeof (double), "psrose");
  257.     yy = (double *) memory (CNULL, n_bins+2, sizeof (double), "psrose");
  258.     azimuth = (double *) memory (CNULL, n_alloc, sizeof (double), "psrose");
  259.     length = (double *) memory (CNULL, n_alloc, sizeof (double), "psrose");
  260.     
  261.     if (gmtdefs.io_header) fscanf (fp, "%[^\n]\n", text);
  262.         
  263.     /* Read data and do some stats */
  264.     
  265.     n = 0;
  266.     ix = gmtdefs.xy_toggle;    iy = 1 - ix;
  267.     
  268.     while (fscanf (fp, "%lf %lf", &xy[ix], &xy[iy]) != EOF) {
  269.         length[n] = xy[0];
  270.         azimuth[n] = xy[1];
  271.         
  272.         if (zscale != 1.0) length[n] *= zscale;
  273.         
  274.         /* Make sure azimuth is in 0 <= az < 360 range */
  275.         
  276.         while (azimuth[n] < 0.0) azimuth[n] += 360.0;
  277.         while (azimuth[n] >= 360.0) azimuth[n] -= 360.0;
  278.         
  279.         if (half_only) {    /* Flip azimuths about E-W line i.e. -90 < az <= 90 */
  280.             if (azimuth[n] > 90.0 && azimuth[n] <= 270.0) azimuth[n] -= 180.0;
  281.             if (azimuth[n] > 270.0) azimuth[n] -= 360.0;
  282.         }
  283.         
  284.         /* Double angle to find mean azimuth */
  285.         
  286.         xr += length[n] * cos (2.0 * azimuth[n] * D2R);
  287.         yr += length[n] * sin (2.0 * azimuth[n] * D2R);
  288.         
  289.         total += length[n];
  290.         
  291.         n++;
  292.         if (n == n_alloc) {    /* Get more memory */
  293.             n_alloc += GMT_CHUNK;
  294.             azimuth = (double *) memory ((char *)azimuth, n_alloc, sizeof (double), "psrose");
  295.             length = (double *) memory ((char *)length, n_alloc, sizeof (double), "psrose");
  296.         }
  297.     }
  298.     if (fp != stdin) fclose (fp);
  299.     
  300.     if (sector > 0.0) {    /* Sum up sector diagram info */
  301.         for (i = 0; i < n; i++) {
  302.             bin = (azimuth[i] + az_offset) / sector;
  303.             sum[bin] += length[i];
  304.         }
  305.     }
  306.         
  307.     mean_theta = 0.5 * R2D * d_atan2 (yr, xr);
  308.     if (mean_theta < 0.0) mean_theta += 360.0;
  309.     mean_radius = hypot (xr, yr) / total;
  310.     if (!normalize) mean_radius *= max_radius;
  311.     
  312.     if (sector > 0.0) {    /* Find max of the bins */
  313.         for (bin = 0; bin < n_bins; bin++) if (sum[bin] > max) max = sum[bin];
  314.         if (normalize) for (bin = 0; bin < n_bins; bin++) sum[bin] /= max;
  315.         if (eq_area) for (bin = 0; bin < n_bins; bin++) sum[bin] = sqrt (sum[bin]);
  316.     }
  317.     else {
  318.         for (i = 0; i < n; i++) if (length[i] > max) max = length[i];
  319.         if (normalize) {
  320.             max = 1.0 / max;
  321.             for (i = 0; i < n; i++) length[i] *= max;
  322.             max = 1.0 / max;
  323.         }
  324.     }
  325.     
  326.     if (inquire || gmtdefs.verbose) {
  327.         fprintf (stderr, "Info for %s: n = %d rmax = %lg mean r/az = (%lg/%lg) totlength = %lg\n",
  328.             file, n, max, mean_radius, mean_theta, total);
  329.         if (inquire) {
  330.             free ((char *)sum);
  331.             free ((char *)xx);
  332.             free ((char *)yy);
  333.             free ((char *)azimuth);
  334.             free ((char *)length);
  335.             exit (-1);
  336.         }
  337.     }
  338.     
  339.     if (automatic) {
  340.         if (frame_info.anot_int[0] == 0.0) {
  341.             tmp = pow (10.0, floor (d_log10 (max)));
  342.             if ((max / tmp) < 3.0) tmp *= 0.5;
  343.         }
  344.         else
  345.             tmp = frame_info.anot_int[0];
  346.         max_radius = ceil (max / tmp) * tmp;
  347.         if (frame_info.anot_int[0] == 0.0 || frame_info.grid_int[0]) {    /* Tickmarks not set */
  348.             frame_info.anot_int[0] = frame_info.grid_int[0] = tmp;
  349.             frame_info.plot = TRUE;
  350.         }
  351.     }
  352.     
  353.     if (frame_info.plot && frame_info.anot_int[1] == 0.0) frame_info.anot_int[1] = frame_info.grid_int[1] = 30.0;
  354.     
  355.     /* Ready to plot */
  356.     
  357.     project_info.projection = LINEAR;
  358.     project_info.region = TRUE;
  359.     project_info.pars[0] = project_info.pars[1] = (normalize) ? 1.0 : scale / max_radius;
  360.     map_setup (-scale, scale, -scale, scale);
  361.     
  362.     ps_plotinit (NULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  363.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  364.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  365.     echo_command (argc, argv);
  366.     if (gmtdefs.unix_time) timestamp (argc, argv);
  367.     
  368.     x_origin = scale;    y_origin = ((half_only) ? 0.0 : scale);
  369.     ps_transrotate (x_origin, y_origin, 0.0);
  370.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin - y_origin, 0.0);
  371.     if (!normalize) scale /= max_radius;
  372.     
  373.     ps_setline (pen.width);
  374.     ps_setpaint (pen.r, pen.g, pen.b);
  375.     if (windrose) {
  376.         for (i = 0; i < n; i++) {
  377.             angle = (90.0 - azimuth[i]) * D2R;
  378.             radius = length[i] * scale;
  379.             x = radius * cos (angle);
  380.             y = radius * sin (angle);
  381.             plot3d (0.0, 0.0, 3);
  382.             plot3d (x, y, 2);
  383.         }
  384.     }
  385.     
  386.     if (sector_plot && !rose) {    /* Draw pie slices for sector plot */
  387.             
  388.         for (bin = 0; bin < n_bins; bin++) {
  389.             az = bin * sector - az_offset;
  390.             angle1 = (90.0 - az - sector);
  391.             angle2 = angle1 + sector;
  392.             pie3d (0.0, 0.0, sum[bin] * scale, angle1, angle2, fill.r, fill.g, fill.b, FALSE);
  393.         }
  394.     }
  395.     else if (rose) {    /* Draw rose diagram */
  396.     
  397.         for (bin = i = 0; bin < n_bins; bin++, i++) {
  398.             az = bin * sector - az_offset + 0.5 * sector;
  399.             angle = 90.0 - az;
  400.             xx[i] = scale * sum[bin] * cos (D2R * angle);
  401.             yy[i] = scale * sum[bin] * sin (D2R * angle);
  402.         }
  403.         if (half_only) {
  404.             xx[i] = scale * 0.5 * (sum[0] + sum[n_bins-1]);
  405.             yy[i] = 0.0;
  406.             i++;
  407.             xx[i] = -xx[i-1];
  408.             yy[i] = 0.0;
  409.             i++;
  410.         }
  411.         if (project_info.three_D) two_d_to_three_d (xx, yy, i);
  412.         ps_polygon (xx, yy, i, fill.r, fill.g, fill.b, outline);
  413.     }
  414.     
  415.     if (sector_plot && outline && !rose) {    /* Draw a line outlining the pie slices */
  416.         ps_setpaint (pen.r, pen.g, pen.b);
  417.         angle1 = (half_only) ? 180.0 : 90.0;
  418.         angle2 = (half_only) ? 0.0 : 90.0;
  419.         angle1 *= D2R;
  420.         angle2 *= D2R;
  421.         x1 = (sum[0] * scale) * cos (angle1);
  422.         y1 = (sum[0] * scale) * sin (angle1);
  423.         x2 = (sum[n_bins-1] * scale) * cos (angle2);
  424.         y2 = (sum[n_bins-1] * scale) * sin (angle2);
  425.         plot3d (x1, y1, 3);
  426.         plot3d (x2, y2, 2);
  427.         for (bin = n_bins-1; bin >= 0; bin--) {
  428.             status = (bin == 0) ? 2 : 0;
  429.             az = bin * sector - az_offset;
  430.             angle1 = 90.0 - az - sector;
  431.             angle2 = angle1 + sector;
  432.             arc3d (0.0, 0.0, sum[bin] * scale, angle1, angle2, status);
  433.             angle1 *= D2R;
  434.         }
  435.     }
  436.     
  437.     if (draw_modes) {
  438.         
  439.         if (find_mean) {    /* Use the mean direction only */
  440.             n_modes = 1;
  441.             mode_direction[0] = mean_theta;
  442.             mode_length[0] = mean_radius;
  443.         }
  444.         else {    /* Get mode parameters from separate file */
  445.             n_modes = 0;
  446.             while (fgets (text, 512, fpm)) {
  447.                 sscanf (text, "%lf %lf", &mode_direction[n_modes], &mode_length[n_modes]);
  448.                 n_modes++;
  449.             }
  450.             fclose (fpm);
  451.         }
  452.         for (i = 0; i < n_modes; i++) {
  453.             if (eq_area) mode_length[i] = sqrt (mode_length[i]);
  454.             if (half_only && mode_direction[i] > 90.0 && mode_direction[i] <= 270.0) mode_direction[i] -= 180.0;
  455.             angle = D2R * (90.0 - mode_direction[i]);
  456.             xr = scale * mode_length[i] * cos (angle);
  457.             yr = scale * mode_length[i] * sin (angle);
  458.             gmt_vector3d (0.0, 0.0, xr, yr, project_info.z_level, tailwidth, headlength, headwidth, gmtdefs.vector_shape, ar, ag, ab, TRUE);
  459.         }
  460.         
  461.     }
  462.     
  463.     if (frame_info.plot) {    /* Draw grid lines etc */
  464.         ps_setline (gmtdefs.grid_pen);
  465.         off = max_radius * scale;
  466.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  467.         n_alpha = (frame_info.grid_int[1] > 0.0) ? total_arc / frame_info.grid_int[1] : -1;
  468.         for (i = 0; i <= n_alpha; i++) {
  469.             angle = D2R * i * frame_info.grid_int[1];
  470.             x = max_radius * scale * cos (angle);
  471.             y = max_radius * scale * sin (angle);
  472.             plot3d (0.0, 0.0, 3);
  473.             plot3d (x, y, 2);
  474.         }
  475.         
  476.         n_bins = (frame_info.grid_int[0] > 0.0) ? rint (max_radius / frame_info.grid_int[0]) : -1;
  477.         for (i = 1; i <= n_bins; i++)
  478.             arc3d (0.0, 0.0, i * frame_info.grid_int[0] * scale, 0.0, total_arc, 3);
  479.         y = lsize + 6.0 * gmtdefs.anot_offset;
  480.         gmt_text3d (0.0, off + y, project_info.z_level, gmtdefs.header_font_size, gmtdefs.header_font, frame_info.header, 0.0, 2, 0);
  481.         
  482.         get_format (frame_info.anot_int[0], format);
  483.         
  484.         if (half_only) {
  485.             strcpy (lw, "90\\312W");    strcpy (le, "90\\312E");    strcpy (ln, "0\\312");
  486.             y = -(3.0 * gmtdefs.anot_offset + font_height[gmtdefs.anot_font] * asize);
  487.             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);
  488.             y = -(5.0 * gmtdefs.anot_offset + font_height[gmtdefs.anot_font] * lsize + font_height[gmtdefs.label_font] * lsize);
  489.             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);
  490.             gmt_text3d (0.0, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, "0", 0.0, 10, 0);
  491.             n_anot = (frame_info.anot_int[0] > 0.0) ? rint (max_radius / frame_info.anot_int[0]) : -1;
  492.             for (i = 1; i <= n_anot; i++) {
  493.                 x = i * frame_info.anot_int[0];
  494.                 sprintf (text, format, x);
  495.                 x *= scale;
  496.                 gmt_text3d (x, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, text, 0.0, 10, 0);
  497.                 gmt_text3d (-x, -gmtdefs.anot_offset, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, text, 0.0, 10, 0);
  498.             }
  499.         }
  500.         else {
  501.             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);
  502.             plot3d (off, -off, 3);
  503.             plot3d ((max_radius - frame_info.grid_int[0]) * scale, -off, 2);
  504.             plot3d (off, -off, 3);
  505.             plot3d (off, gmtdefs.tick_length - off, 2);
  506.             plot3d ((max_radius - frame_info.grid_int[0]) * scale, -off, 3);
  507.             plot3d ((max_radius - frame_info.grid_int[0]) * scale, gmtdefs.tick_length - off, 2);
  508.             y = -(off + 5.0 * gmtdefs.anot_offset + font_height[gmtdefs.anot_font] * lsize + font_height[gmtdefs.label_font] * lsize);
  509.             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);
  510.             if (frame_info.label[0][0]) {
  511.                 strcat (format, " %s");
  512.                 sprintf (text, format, frame_info.grid_int[0], frame_info.label[0]);
  513.             }
  514.             else
  515.                 sprintf (text, format, frame_info.grid_int[0]);
  516.             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);
  517.         }
  518.         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);
  519.         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);
  520.         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);
  521.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  522.     }
  523.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin + y_origin, 0.0);
  524.     ps_rotatetrans (-x_origin, -y_origin, 0.0);
  525.  
  526.     ps_plotend (gmtdefs.last_page);
  527.     
  528.     gmt_end (argc, argv);
  529.     
  530.     free ((char *)sum);
  531.     free ((char *)xx);
  532.     free ((char *)yy);
  533.     free ((char *)azimuth);
  534.     free ((char *)length);
  535.     
  536.     exit (0);
  537. }
  538.  
  539. int arc3d (x, y, r, a1, a2, status)
  540. double x, y, r, a1, a2;
  541. int status; {
  542.     int i, n;
  543.     double da, a, *xx, *yy;
  544.     
  545.     if (project_info.three_D) {    /* Must project and draw line segments */
  546.         n = ceil (a2 - a1);
  547.         
  548.         xx = (double *) memory (CNULL, n + 1, sizeof (double), "arc3d");
  549.         yy = (double *) memory (CNULL, n + 1, sizeof (double), "arc3d");
  550.         
  551.         da = (a2 - a1) / n;
  552.         for (i = 0; i <= n; i++) {
  553.             a = (i == n) ? a2 : a1 + da * i;
  554.             xx[i] = x + r * cosd (a);
  555.             yy[i] = y + r * sind (a);
  556.         }
  557.         n++;    
  558.         two_d_to_three_d (xx, yy, n);
  559.         ps_line (xx, yy, n, status, FALSE, TRUE);
  560.         free ((char *)xx);
  561.         free ((char *)yy);
  562.     }
  563.     else
  564.         ps_arc (x, y, r, a1, a2, status);
  565. }
  566.  
  567. int plot3d (x, y, pen)
  568. double x, y;
  569. int pen; {
  570.     if (project_info.three_D) xyz_to_xy (x, y, project_info.z_level, &x, &y);
  571.     ps_plot (x, y, pen);
  572. }
  573.  
  574. int pie3d (x, y, rad, a1, a2, r, g, b, outline)
  575. double x, y, rad, a1, a2;
  576. int r, g, b;
  577. BOOLEAN outline; {
  578.     int i, n;
  579.     double da, a, *xx, *yy;
  580.     
  581.     if (project_info.three_D) {    /* Must project and draw line segments */
  582.         n = ceil (a2 - a1);
  583.         
  584.         xx = (double *) memory (CNULL, n + 2, sizeof (double), "arc3d");
  585.         yy = (double *) memory (CNULL, n + 2, sizeof (double), "arc3d");
  586.         
  587.         da = (a2 - a1) / n;
  588.         for (i = 0; i <= n; i++) {
  589.             a = (i == n) ? a2 : a1 + da * i;
  590.             xx[i] = x + rad * cosd (a);
  591.             yy[i] = y + rad * sind (a);
  592.         }
  593.         xx[i] = x;    yy[i] = y;
  594.         n += 2;    
  595.         two_d_to_three_d (xx, yy, n);
  596.         ps_polygon (xx, yy, n, r, g, b, outline);
  597.         free ((char *)xx);
  598.         free ((char *)yy);
  599.     }
  600.     else
  601.         ps_pie (x, y, rad, a1, a2, r, g, b, outline);
  602. }
  603.