home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / psxyz.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-11  |  25.9 KB  |  773 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)psxyz.c    2.35  4/11/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.  * psxyz will read <x,y,z> triplets from stdin and plot
  9.  * the positions in 3-D using symbols selected by the user. If no
  10.  * symbol size is specified, psxyz expects <x,y,z,size> as input, where
  11.  * size is the symbol size.  Several map projections are supported.
  12.  * For linear projections a 3-D basemap is provided.  All symbols are
  13.  * projected onto the xy plane (so that circles becomes ellipses) except
  14.  * BAR and CUBE which is fully 3-D.
  15.  * PostScript code is then written to stdout.
  16.  *
  17.  * Author:    Paul Wessel
  18.  * Date:      15-JAN-1991-1995
  19.  * Version:   2.0, based on old v1.1
  20.  *
  21.  */
  22.  
  23. #include "gmt.h"
  24.  
  25. #define LINE 0
  26. #define BAR 1
  27. #define CROSS 2
  28. #define POINT 3
  29. #define CIRCLE 4
  30. #define SQUARE 5
  31. #define TRIANGLE 6
  32. #define DIAMOND 7
  33. #define CUBE 8
  34. #define COLUMN 9
  35. #define VECTOR 10
  36. #define VECTOR2 11
  37.  
  38. #define POINTSIZE 0.005
  39.  
  40. double *xp, *yp;
  41.  
  42. int compare();
  43.  
  44. struct DATA1 {
  45.     double x, y, z, value, dist;
  46.     float x_size, y_size;
  47. } *data1;
  48.  
  49. struct DATA2 {
  50.     double x, y, z;
  51. } *data2;
  52.  
  53. main (argc, argv)
  54. int argc;
  55. char **argv; {
  56.     int     i, j, symbol = 0, n, ix = 0, iy = 1, n_files = 0, fno, red[3], green[3], blue[3];
  57.     int    n_alloc = GMT_CHUNK, n_read, n_col, n_args, three, four, five, bset = -1, n_expected = 3;
  58.     
  59.     BOOLEAN    error = FALSE, nofile = TRUE, polygon = FALSE, done, read_vector = FALSE;
  60.     BOOLEAN    read_size = FALSE, outline = FALSE, multi_segments = FALSE, get_rgb = FALSE;
  61.     BOOLEAN convert_angles = FALSE, skip_if_outside = TRUE;
  62.     
  63.     double xy[2], west = 0.0, east = 0.0, south = 0.0, north = 0.0, new_z_level = 0.0;
  64.     double symbol_size_x = 0.0, symbol_size_y = 0.0, lux[3], x, y, z, tmp, base;
  65.     double x2, y2, v_width = 0.03, h_length = 0.12, h_width = 0.1, v_w, h_l, h_w, v_shrink, v_norm;
  66.     
  67.     char line[512], symbol_type, col[7][100], *cpt, *more, EOL_flag = '>';
  68.     
  69.     FILE *fp = NULL;
  70.     
  71.     struct PEN pen;
  72.     struct FILL fill;
  73.     
  74.     argc = gmt_begin (argc, argv);
  75.     
  76.     gmt_init_pen (&pen, 1);
  77.     gmt_init_fill (&fill, -1, -1, -1);      /* Default is no fill */
  78.  
  79.     for (i = 0; i < 3; i++) red[i] = green[i] = blue[i] = -1;
  80.     base = gmt_NaN;
  81.  
  82.     /* Check and interpret the command line arguments */
  83.     
  84.     if (gmtdefs.measure_unit) {
  85.         v_width = 0.075;    h_length = 0.3;    h_width = 0.25;
  86.     }
  87.  
  88.     for (i = 1; i < argc; i++) {
  89.         if (argv[i][0] == '-') {
  90.             switch(argv[i][1]) {
  91.                 /* Common parameters */
  92.             
  93.                 case 'B':
  94.                 case 'H':
  95.                 case 'J':
  96.                 case 'K':
  97.                 case 'O':
  98.                 case 'P':
  99.                 case 'R':
  100.                 case 'U':
  101.                 case 'V':
  102.                 case 'X':
  103.                 case 'x':
  104.                 case 'Y':
  105.                 case 'y':
  106.                 case 'c':
  107.                 case ':':
  108.                 case '\0':
  109.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  110.                     break;
  111.                 
  112.                 /* Supplemental parameters */
  113.             
  114.                 case 'C':
  115.                     cpt = &argv[i][2];
  116.                     get_rgb = TRUE;
  117.                     break;
  118.                 case 'E':
  119.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  120.                     break;
  121.                 case 'F':
  122.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  123.                         gmt_pen_syntax ('F');
  124.                         error++;
  125.                     }
  126.                     break;
  127.                 case 'G':        /* Set color for polygon */
  128.                     if (gmt_getfill (&argv[i][2], &fill)) {
  129.                         gmt_fill_syntax ('G');
  130.                         error++;
  131.                     }
  132.                     polygon = TRUE;
  133.                     break;
  134.                 case 'L':        /* Draw the outline */
  135.                     outline = TRUE;
  136.                     break;
  137.                 case 'M':        /* Multiple line segments */
  138.                     multi_segments = TRUE;
  139.                     if (argv[i][2]) EOL_flag = argv[i][2];
  140.                     break;
  141.                                 case 'N':               /* Do not clip to map */
  142.                                         skip_if_outside = FALSE;
  143.                                         break;
  144.                 case 'S':        /* Get symbol [and size] */
  145.                     sscanf (&argv[i][2], "%c%lf/%lf", &symbol_type, &symbol_size_x, &symbol_size_y);
  146.                     if (symbol_size_y == 0.0) symbol_size_y = symbol_size_x;
  147.                     if (symbol_size_x == 0.0 && symbol != POINT) read_size = TRUE;
  148.                     symbol_size_x *= 0.5;
  149.                     symbol_size_y *= 0.5;
  150.                     switch (symbol_type) {
  151.                         case 'b':
  152.                             symbol = BAR;
  153.                             for (j = 3; argv[i][j] && bset < 0; j++) if (argv[i][j] == 'b') bset = j;
  154.                             if (bset > 0) base = atof (&argv[i][j]);
  155.                             break;
  156.                         case 'c':
  157.                             symbol = CIRCLE;
  158.                             symbol_size_x *= 2.0;
  159.                             break;
  160.                         case 'd':
  161.                             symbol = DIAMOND;
  162.                             break;
  163.                         case 'o':
  164.                             symbol = COLUMN;
  165.                             for (j = 3, bset = -1; argv[i][j] && bset < 0; j++) if (argv[i][j] == 'b') bset = j;
  166.                             if (bset > 0) base = atof (&argv[i][j]);
  167.                             break;
  168.                         case 'p':
  169.                             symbol = POINT;
  170.                             break;
  171.                         case 's':
  172.                             symbol = SQUARE;
  173.                             break;
  174.                         case 't':
  175.                             symbol = TRIANGLE;
  176.                             break;
  177.                         case 'u':
  178.                             symbol = CUBE;
  179.                             break;
  180.                         case 'V':
  181.                             convert_angles = TRUE;
  182.                         case 'v':
  183.                             symbol = VECTOR;
  184.                             read_size = FALSE;
  185.                             n_expected += 2;
  186.                             if (argv[i][3]) sscanf (&argv[i][3], "%lf/%lf/%lf", &v_width, &h_length, &h_width);
  187.                             for (j = 3; argv[i][j] && argv[i][j] != 'n'; j++);
  188.                             if (argv[i][j]) {       /* Normalize option used */
  189.                                 v_norm = atof (&argv[i][j+1]);
  190.                                 if (v_norm > 0.0) {
  191.                                     v_shrink = 1.0 / v_norm;
  192.                                     symbol = VECTOR2;
  193.                                 }
  194.                             }
  195.                             read_vector = TRUE;
  196.                             break;
  197.                         case 'x':
  198.                             symbol = CROSS;
  199.                             break;
  200.                         default:
  201.                             error = TRUE;
  202.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Unrecognized symbol type %c\n", gmt_program, symbol_type);
  203.                             break;
  204.                     }
  205.                     if (read_size) n_expected++;
  206.                     break;
  207.                 case 'W':        /* Set line attributes */
  208.                     if (gmt_getpen (&argv[i][2], &pen)) {
  209.                         gmt_pen_syntax ('W');
  210.                         error++;
  211.                     }
  212.                     break;
  213.                 case 'Z':
  214.                     new_z_level = atof (&argv[i][2]);
  215.                     break;
  216.  
  217.                 default:        /* Options not recognized */
  218.                     error = TRUE;
  219.                     gmt_default_error (argv[i][1]);
  220.                     break;
  221.             }
  222.         }
  223.         else
  224.             n_files++;
  225.     }
  226.     
  227.     if (argc == 1 || gmt_quick) {    /* Display usage */
  228.         fprintf (stderr,"psxyz %s - Plot lines, polygons, and symbols in 3-D\n\n", GMT_VERSION);
  229.         fprintf(stderr,"usage: psxyz <xyzfiles> -R<west/east/south/north/zmin/zmax> -J<params>\n");
  230.         fprintf(stderr, "    -Jz<params> [-B<tickinfo>] [-C<cpt>] [-E<az/el>] [-F<r/g/b>] [-G<fill>] [-H] [-K] [-L]\n");
  231.         fprintf(stderr, "    [-M[<flag>]] [-N] [-O] [-P] [-S<symbol><size>[/size_y]] [-U] [-V] [-W<pen>] [-X<x_shift>]\n");
  232.         fprintf(stderr, "    [-Y<y_shift>] [-c<ncopies>]\n");
  233.         
  234.         if (gmt_quick) exit(-1);
  235.         
  236.         fprintf (stderr, "    <xyzfiles> is one or more files.  If none, read standard input\n");
  237.         explain_option ('j');
  238.         explain_option ('Z');
  239.         explain_option ('r');
  240.         fprintf (stderr, "\n\tOPTIONS:\n");
  241.         explain_option ('b');
  242.         fprintf (stderr, "    -C Use cpt-file to assign colors based on z-value in 3rd column\n");
  243.         fprintf (stderr, "       Must be used with -S\n");
  244.         fprintf (stderr, "    -E set user viewpoint azimuth and elevation [180/90].\n");
  245.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d].\n",
  246.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  247.         fprintf (stderr, "    -G Specify color (for symbols/polygons) or pattern (for polygons). fill can be either\n");
  248.         fprintf (stderr, "       1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
  249.         fprintf (stderr, "       2) p[or P]<iconsize>/<pattern> for predefined patterns (0-31).\n");
  250.         fprintf (stderr, "       If -G is specified but not -S, psxyz draws a filled polygon.\n");
  251.         fprintf (stderr, "       Default is no fill (transparent symbols or polygons)\n");
  252.         explain_option ('H');
  253.         explain_option ('K');
  254.         fprintf (stderr, "    -L close polygon OR draw symbol outline with current pen (see -W).\n");
  255.         fprintf (stderr, "    -M Input files each consist of multiple segments separated by one record\n");
  256.         fprintf (stderr, "       whose first character is <flag> [>].\n");
  257.         fprintf (stderr, "    -N Do Not skip symbols that fall outside map border [Default will ignore those outside]\n");
  258.         explain_option ('O');
  259.         explain_option ('P');
  260.         fprintf (stderr, "    -S to select symbol type and symbol size (in %s).  Choose between\n", gmt_unit_names[gmtdefs.measure_unit]);
  261.         fprintf (stderr, "       (b)ar, (c)ircle, (d)iamond, c(o)lumn, (p)oint, (s)quare, (t)riangle, c(u)be, (v)ector, (x)cross\n");
  262.         fprintf (stderr, "       If no size is specified, psxyz expects the 4th column to have sizes.\n");
  263.         fprintf (stderr, "       [Note: if -C is selected then 4th means 5th column, etc.]\n");
  264.         fprintf (stderr, "       column and cube are true 3-D objects (give size as xsize/ysize), the rest is 2-D perspective only.\n");
  265.         fprintf (stderr, "       Bar (or Column): Append b<base> to give the y- (or z-) value of the base [Default = 0 (1 for log-scales)]\n");
  266.         fprintf (stderr, "         Vectors: the direction and length must be in input columns 4 and 5.\n");
  267.         fprintf (stderr, "         Furthermore, <size> means arrowwidth/headlength/headwith [Default is %lg/%lg/%lg]\n", v_width, h_length, h_width);
  268.         fprintf (stderr, "           If -SV rather than -Sv is selected, psxy will expect azimuth and length\n");
  269.         fprintf (stderr, "           and convert azimuths based on the chosen map projection\n");        explain_option ('U');
  270.         explain_option ('V');
  271.         fprintf (stderr, "    -W sets pen attributes [width = %d, color = (%d/%d/%d), solid line].\n", 
  272.             pen.width, pen.r, pen.g, pen.b);
  273.         explain_option ('X');
  274.         explain_option ('c');
  275.         explain_option (':');
  276.         explain_option ('.');
  277.         exit(-1);
  278.     }
  279.  
  280.     /* Check that the options selected are mutually consistant */
  281.     
  282.     if (!project_info.region_supplied) {
  283.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  284.         error++;
  285.     }
  286.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  287.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  288.         error++;
  289.     }
  290.     if (get_rgb && symbol == 0) {
  291.         error++;
  292.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Must also specify a symbol (see -S)\n", gmt_program);
  293.     }
  294.     
  295.     if (error) exit (-1);
  296.     
  297.     if (get_rgb) read_cpt (cpt);
  298.     
  299.     if (symbol == 0 && outline) polygon = TRUE;
  300.  
  301.     if (bset < 0) base = (project_info.xyz_projection[2] == LOG10) ? 1.0 : 0.0;
  302.     
  303.     if (n_files > 0)
  304.         nofile = FALSE;
  305.     else
  306.         n_files = 1;
  307.     n_args = (argc > 1) ? argc : 2;
  308.     
  309.     map_setup (west, east, south, north);
  310.  
  311.     if (fill.r >= 0 && (symbol == COLUMN || symbol == CUBE)) {    /* Modify the color for each facet */
  312.         lux[0] = fabs (z_project.sin_az * z_project.cos_el);
  313.         lux[1] = fabs (z_project.cos_az * z_project.cos_el);
  314.         lux[2] = fabs (z_project.sin_el);
  315.         tmp = MAX (lux[0], MAX (lux[1], lux[2]));
  316.         for (i = 0; i < 3; i++) {
  317.             lux[i] = (lux[i] / tmp) - 0.5;
  318.             red[i] = fill.r;
  319.             green[i] = fill.g;
  320.             blue[i] = fill.b;
  321.             illuminate (lux[i], &red[i], &green[i], &blue[i]);
  322.         }
  323.     }
  324.     
  325.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  326.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  327.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  328.     echo_command (argc, argv);
  329.     if (gmtdefs.unix_time) timestamp (argc, argv);
  330.  
  331.     ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  332.     
  333.     if (new_z_level != 0.0) project_info.z_level = new_z_level;
  334.     z_to_zz (base, &tmp);    base = tmp;
  335.  
  336.     if (frame_info.plot) {    /* First plot background axes */
  337.         frame_info.header[0] = 0;    /* No header for grdview and psxyz */
  338.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  339.         map_basemap ();
  340.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  341.     }
  342.     ps_setline (pen.width);
  343.     if (pen.texture) ps_setdash (pen.texture, pen.offset);
  344.     ps_setpaint (pen.r, pen.g, pen.b);
  345.     
  346.     three = (get_rgb) ? 4 : 3;
  347.     four = (get_rgb) ? 5 : 4;
  348.     five = (get_rgb) ? 6 : 5;
  349.     
  350.     ix = (gmtdefs.xy_toggle);    iy = 1 - ix;
  351.     done = FALSE;
  352.     for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over all input files */
  353.         if (!nofile && argv[fno][0] == '-') continue;
  354.         if (nofile) {
  355.             fp = stdin;
  356.             done = TRUE;
  357.         }
  358.         else if ((fp = fopen (argv[fno], "r")) == NULL) {
  359.             fprintf (stderr, "psxy: Cannot open file %s\n", argv[fno]);
  360.             continue;
  361.         }
  362.  
  363.         if (!nofile && gmtdefs.verbose) {
  364.             fprintf (stderr, "psxyz: Working on file %s\n", argv[fno]);
  365.             sprintf (line, "File: %s\0", argv[fno]);
  366.             ps_comment (line);
  367.         }
  368.         if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
  369.         
  370.         if (multi_segments) {
  371.             fgets (line, 512, fp);
  372.             if (gmtdefs.verbose) ps_comment (line);
  373.         }
  374.  
  375.         if (symbol > 0) {    /* symbol part */
  376.         
  377.             gmt_world_map = TRUE;
  378.         
  379.             data1 = (struct DATA1 *) memory (CNULL, n_alloc, sizeof (struct DATA1), "psxyz");
  380.             n = 0;
  381.             while (fgets (line, 512, fp)) {
  382.                 n_col = sscanf (line, "%s %s %s %s %s %s", col[0], col[1], col[2], col[3], col[4], col[5]);
  383.             
  384.                 if (multi_segments && col[0][0] == EOL_flag) continue;
  385.             
  386.                 if (n_col < n_expected) {
  387.                     fprintf (stderr, "psxyz: Mismatch between expected (%d) and actual (%d) columns!\n", n_expected, n_col);
  388.                     fprintf (stderr, "--> %s\n", line);
  389.                     exit (-1);
  390.                 }
  391.                 x = atof (col[ix]);
  392.                 y = atof (col[iy]);
  393.                 z = atof (col[2]);
  394.             
  395.                 if (skip_if_outside) {
  396.                     map_outside (x, y);
  397.                     if ( abs (gmt_x_status_new) > 1 || abs (gmt_y_status_new) > 1) continue;
  398.                 }
  399.  
  400.                 project3D (x, y, z, &data1[n].x, &data1[n].y, &data1[n].z);
  401.                 if (get_rgb) data1[n].value = atof (col[3]);
  402.                 if (read_size) {
  403.                     data1[n].x_size = atof (col[three]);
  404.                     data1[n].y_size = (n_col == five) ? atof (col[four]) : data1[n].x_size;
  405.                 }
  406.                 else if (read_vector) {
  407.                     data1[n].x_size = atof (col[three]);    /* direction */
  408.                     data1[n].y_size = atof (col[four]);    /* length */
  409.                 }
  410.                 else {
  411.                     data1[n].x_size = symbol_size_x;
  412.                     data1[n].y_size = symbol_size_y;
  413.                 }
  414.                 n++;
  415.                 if (n == n_alloc) {
  416.                     n_alloc += GMT_CHUNK;
  417.                     data1 = (struct DATA1 *) memory ((char *)data1, n_alloc, sizeof (struct DATA1), "psxyz");
  418.                 }
  419.             }
  420.             data1 = (struct DATA1 *) memory ((char *)data1, n, sizeof (struct DATA1), "psxyz");
  421.             
  422.             /* Sort according to distance  from viewer */
  423.             
  424.             sort_on_distance (data1, n);
  425.             
  426.             for (i = 0; i < n; i++) {
  427.                 if (get_rgb) {
  428.                     get_rgb24 (data1[i].value, &fill.r, &fill.g, &fill.b);
  429.                     if (symbol == COLUMN || symbol == CUBE) {
  430.                         for (j = 0; j < 3; j++) {
  431.                             red[j] = fill.r;
  432.                             green[j] = fill.g;
  433.                             blue[j] = fill.b;
  434.                             illuminate (lux[j], &red[j], &green[j], &blue[j]);
  435.                         }
  436.                     }
  437.                 }
  438.  
  439.                 switch (symbol) {
  440.                     case BAR:
  441.                         bar3D (data1[i].x, data1[i].y, data1[i].z, base, (double)data1[i].x_size, fill.r, fill.g, fill.b, outline);
  442.                         break;
  443.                     case COLUMN:
  444.                         column3D (data1[i].x, data1[i].y, data1[i].z, base, (double)data1[i].x_size, (double)data1[i].y_size, red, green, blue, outline);
  445.                         break;
  446.                     case CROSS:
  447.                         cross3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size);
  448.                         break;
  449.                     case POINT:
  450.                         cross3D (data1[i].x, data1[i].y, data1[i].z, POINTSIZE);
  451.                         break;
  452.                     case CIRCLE:
  453.                         circle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.r, fill.g, fill.b, outline);
  454.                         break;
  455.                     case SQUARE:
  456.                         square3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.r, fill.g, fill.b, outline);
  457.                         break;
  458.                     case TRIANGLE:
  459.                         triangle3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.r, fill.g, fill.b, outline);
  460.                         break;
  461.                     case DIAMOND:
  462.                         diamond3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, fill.r, fill.g, fill.b, outline);
  463.                         break;
  464.                     case CUBE:
  465.                         cube3D (data1[i].x, data1[i].y, data1[i].z, (double)data1[i].x_size, (double)data1[i].y_size, red, green, blue, outline);
  466.                         break;
  467.                     case VECTOR:
  468.                         if (data1[i].y_size <= 0.0) continue;
  469.                         if (convert_angles) {
  470.                             azim_2_angle (data1[i].x, data1[i].y, 0.1, data1[i].x_size, &tmp);
  471.                             data1[i].x_size = tmp;
  472.                         }
  473.                         x2 = data1[i].x + data1[i].y_size * cosd (data1[i].x_size);
  474.                         y2 = data1[i].y + data1[i].y_size * sind (data1[i].x_size);
  475.                         gmt_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, v_width, h_length, h_width, gmtdefs.vector_shape, fill.r, fill.g, fill.b, outline);
  476.                         break;
  477.                     case VECTOR2:
  478.                         if (data1[i].y_size <= 0.0) continue;
  479.                         if (convert_angles) {
  480.                             azim_2_angle (data1[i].x, data1[i].y, 0.1, data1[i].x_size, &tmp);
  481.                             data1[i].x_size = tmp;
  482.                         }
  483.                         x2 = data1[i].x + data1[i].y_size * cosd (data1[i].x_size);
  484.                         y2 = data1[i].y + data1[i].y_size * sind (data1[i].x_size);
  485.                         if (data1[i].y_size < v_norm) {    /* Scale arrow attributes down with length */
  486.                             v_w = v_width * data1[i].y_size * v_shrink;
  487.                             h_l = h_length * data1[i].y_size * v_shrink;
  488.                             h_w = h_width * data1[i].y_size * v_shrink;
  489.                             gmt_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, v_w, h_l, h_w, gmtdefs.vector_shape, fill.r, fill.g, fill.b, outline);
  490.                         }
  491.                         else    /* Leave as specified */
  492.                             gmt_vector3d (data1[i].x, data1[i].y, x2, y2, data1[i].z, v_width, h_length, h_width, gmtdefs.vector_shape, fill.r, fill.g, fill.b, outline);
  493.                         break;
  494.                 }
  495.             }
  496.             free ((char *)data1);
  497.         }
  498.         else {    /* Line/polygon part */
  499.             data2 = (struct DATA2 *) memory (CNULL, n_alloc, sizeof (struct DATA2), "psxyz");
  500.             more = fgets (line, 512, fp);
  501.             n = 0;
  502.             while (more) {
  503.                 while ((more && !multi_segments) || (more && multi_segments && line[0] != EOL_flag)) {
  504.                     n_read = sscanf (line, "%lf %lf %lf", &xy[ix], &xy[iy], &z);
  505.                     data2[n].x = xy[0];    data2[n].y = xy[1];
  506.                     data2[n].z = (n_read == 2) ? 0.0 : z;
  507.                     n++;
  508.                     if (n == n_alloc) {
  509.                         n_alloc += GMT_CHUNK;
  510.                         data2 = (struct DATA2 *) memory ((char *)data2, n_alloc, sizeof (struct DATA2), "psxyz");
  511.                     }
  512.                     more = fgets (line, 512, fp);
  513.                 }
  514.         
  515.                 if (polygon) {  /* Explicitly close polygon */
  516.                     data2[n].x = data2[0].x;
  517.                     data2[n].y = data2[0].y;
  518.                     data2[n].z = data2[0].z;
  519.                     n++;
  520.                 }
  521.                 data2 = (struct DATA2 *) memory ((char *)data2, n, sizeof (struct DATA2), "psxyz");
  522.                 n_alloc = n;
  523.                 
  524.                 xp = (double *) memory (CNULL, n, sizeof (double), "psxyz");
  525.                 yp = (double *) memory (CNULL, n, sizeof (double), "psxyz");
  526.                 for (i = 0; i < n; i++) geoz_to_xy (data2[i].x, data2[i].y, data2[i].z, &xp[i], &yp[i]);
  527.                 
  528.                 if (polygon) {
  529.                     if (fill.use_pattern)
  530.                         ps_imagefill (xp, yp, n, fill.pattern_no, fill.pattern, fill.inverse, fill.icon_size, FALSE);
  531.                     else
  532.                         ps_polygon (xp, yp, n, fill.r, fill.g, fill.b, outline);
  533.                 }
  534.                 else
  535.                     ps_line (xp, yp, n, 3, FALSE, TRUE);
  536.         
  537.                 free ((char *)xp);
  538.                 free ((char *)yp);
  539.                 n = 0;
  540.             
  541.                 if (more) more = fgets (line, 512, fp);
  542.             }
  543.             free ((char *)data2);
  544.         }
  545.         if (fp != stdin) fclose (fp);
  546.         
  547.     }
  548.     
  549.     if (pen.texture) ps_setdash (CNULL, 0);
  550.  
  551.     if (project_info.three_D) {    /* Redraw front axes */
  552.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  553.         vertical_axis (2);
  554.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  555.     }
  556.     
  557.     ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  558.     ps_plotend (gmtdefs.last_page);
  559.     
  560.     gmt_end (argc, argv);
  561. }
  562.  
  563. int column3D (x, y, z, base, x_size, y_size, r, g, b, outline)
  564. double x, y, z, base, x_size, y_size;
  565. int r[3], g[3], b[3], outline; {
  566.     int i, j, k;
  567.     double zz, xp[4], yp[4], zp[4], plot_x[4], plot_y[4], top, sign;
  568.     
  569.     top = z;
  570.     if (top < base) d_swap (top, base);
  571.     
  572.     for (i = 0; i < 3; i++) {
  573.         sign = -1.0;
  574.         zz = base;
  575.         switch (z_project.face[i]) {
  576.             case 0:    /* yz plane positive side */
  577.                 sign = 1.0;
  578.             case 1:    /* negative side */
  579.                 xp[0] = xp[1] = xp[2] = xp[3] = x + sign * x_size;
  580.                 yp[0] = yp[3] = y - y_size;    yp[1] = yp[2] = y + y_size;
  581.                 zp[0] = zp[1] = base;    zp[2] = zp[3] = top;
  582.                 break;
  583.             case 2:    /* xz plane positive side */
  584.                 sign = 1.0;
  585.             case 3:    /* negative side */
  586.                 xp[0] = xp[3] = x - x_size;    xp[1] = xp[2] = x + x_size;
  587.                 yp[0] = yp[1] = yp[2] = yp[3] = y + sign * y_size;
  588.                 zp[0] = zp[1] = base;    zp[2] = zp[3] = top;
  589.                 break;
  590.             case 4:    /* xy plane positive side */
  591.                 zz = top;
  592.             case 5:    /* negative side */
  593.                 xp[0] = xp[3] = x - x_size;    yp[0] = yp[1] = y - y_size;
  594.                 xp[1] = xp[2] = x + x_size;    yp[2] = yp[3] = y + y_size;
  595.                 zp[0] = zp[1] = zp[2] = zp[3] = zz;
  596.                 break;
  597.         }
  598.         k = z_project.face[i] / 2;
  599.         for (j = 0; j < 4; j++) xyz_to_xy (xp[j], yp[j], zp[j], &plot_x[j], &plot_y[j]);
  600.         ps_patch (plot_x, plot_y, 4, r[k], g[k], b[k], outline);
  601.     }
  602. }
  603.  
  604. int cube3D (x, y, z, x_size, y_size, r, g, b, outline)
  605. double x, y, z, x_size, y_size;
  606. int r[3], g[3], b[3], outline; {
  607.     int i, j, k;
  608.     double xp[4], yp[4], zp[4], plot_x[4], plot_y[4], sign;
  609.     
  610.     for (i = 0; i < 3; i++) {
  611.         sign = -1.0;
  612.         switch (z_project.face[i]) {
  613.             case 4:    /* xy plane positive side */
  614.                 sign = 1.0;
  615.             case 5:    /* negative side */
  616.                 xp[0] = xp[3] = x - x_size;    yp[0] = yp[1] = y - y_size;
  617.                 xp[1] = xp[2] = x + x_size;    yp[2] = yp[3] = y + y_size;
  618.                 zp[0] = zp[1] = zp[2] = zp[3] = z + sign * x_size;
  619.                 break;
  620.             case 2:    /* xz plane positive side */
  621.                 sign = 1.0;
  622.             case 3:    /* negative side */
  623.                 xp[0] = xp[3] = x - x_size;    xp[1] = xp[2] = x + x_size;
  624.                 yp[0] = yp[1] = yp[2] = yp[3] = y + sign * y_size;
  625.                 zp[0] = zp[1] = z - x_size;    zp[2] = zp[3] = z + x_size;
  626.                 break;
  627.             case 0:    /* yz plane positive side */
  628.                 sign = 1.0;
  629.             case 1:    /* negative side */
  630.                 xp[0] = xp[1] = xp[2] = xp[3] = x + sign * x_size;
  631.                 yp[0] = yp[3] = y - y_size;    yp[1] = yp[2] = y + y_size;
  632.                 zp[0] = zp[1] = z - x_size;    zp[2] = zp[3] = z + x_size;
  633.                 break;
  634.         }
  635.         k = z_project.face[i] / 2;
  636.         for (j = 0; j < 4; j++) xyz_to_xy (xp[j], yp[j], zp[j], &plot_x[j], &plot_y[j]);
  637.         ps_patch (plot_x, plot_y, 4, r[k], g[k], b[k], outline);
  638.     }
  639. }
  640.  
  641. int cross3D (x, y, z, size)
  642. double x, y, z, size; {
  643.     double xp[2], yp[2], plot_x, plot_y;
  644.     
  645.     xp[0] = x - size;    xp[1] = x + size;
  646.     yp[0] = y - size;    yp[1] = y + size;
  647.     xyz_to_xy (xp[0], yp[0], z, &plot_x, &plot_y);
  648.     ps_plot (plot_x, plot_y, 3);
  649.     xyz_to_xy (xp[1], yp[1], z, &plot_x, &plot_y);
  650.     ps_plot (plot_x, plot_y, 2);
  651.     xyz_to_xy (xp[1], yp[0], z, &plot_x, &plot_y);
  652.     ps_plot (plot_x, plot_y, 3);
  653.     xyz_to_xy (xp[0], yp[1], z, &plot_x, &plot_y);
  654.     ps_plot (plot_x, plot_y, 2);
  655. }
  656.     
  657. int bar3D (x, y, z, base, size, r, g, b, outline)
  658. double x, y, z, base, size;
  659. int r, g, b, outline; {
  660.     int i;
  661.     double xp[4], yp[4], plot_x[4], plot_y[4];
  662.     
  663.     xp[0] = xp[3] = x - size;    xp[1] = xp[2] = x + size;
  664.     yp[0] = yp[1] = base;    yp[2] = yp[3] = y;
  665.     for (i = 0; i < 4; i++) xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]);
  666.     ps_patch (plot_x, plot_y, 4, r, g, b, outline);
  667. }
  668.  
  669. int square3D (x, y, z, size, r, g, b, outline)
  670. double x, y, z, size;
  671. int r, g, b, outline; {
  672.     int i;
  673.     double xp[4], yp[4], plot_x[4], plot_y[4];
  674.     
  675.     xp[0] = xp[3] = x - size;    xp[1] = xp[2] = x + size;
  676.     yp[0] = yp[1] = y - size;    yp[2] = yp[3] = y + size;
  677.     for (i = 0; i < 4; i++) xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]);
  678.     ps_patch (plot_x, plot_y, 4, r, g, b, outline);
  679. }
  680.  
  681. int circle3D (x, y, z, size, r, g, b, outline)
  682. double x, y, z, size;
  683. int r, g, b, outline; {
  684.     /* Must plot a squashed circle */
  685.     int i;
  686.     double xx, yy, a, da, plot_x[51], plot_y[51];
  687.     
  688.     da = 2.0 * M_PI / 50.0;
  689.     for (i = 0; i <= 50; i++) {
  690.         a = i * da;
  691.         xx = x + size * cos (a);
  692.         yy = y + size * sin (a);
  693.         xyz_to_xy (xx, yy, z, &plot_x[i], &plot_y[i]);
  694.     }
  695.     ps_polygon (plot_x, plot_y, 51, r, g, b, outline);
  696. }
  697.  
  698. int triangle3D (x, y, z, size, r, g, b, outline)
  699. double x, y, z, size;
  700. int r, g, b, outline; {
  701.     int i;
  702.     double xp[3], yp[3], plot_x[3], plot_y[3];
  703.     
  704.     xp[0] = x - size;    yp[0] = yp[1] = y - 0.5773502  * size;
  705.     xp[1] = x + size;    xp[2] = x;     yp[2] = y + 1.1547004 * size;
  706.     for (i = 0; i < 3; i++) xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]);
  707.     ps_patch (plot_x, plot_y, 3, r, g, b, outline);
  708. }
  709.  
  710. int diamond3D (x, y, z, size, r, g, b, outline)
  711. double x, y, z, size;
  712. int r, g, b, outline; {
  713.     int i;
  714.     double xp[4], yp[4], plot_x[4], plot_y[4];
  715.     
  716.     xp[0] = xp[2] = x;    xp[1] = x - size;    xp[3] = x + size;
  717.     yp[0] = y - size;    yp[1] = yp[3] = y;    yp[2] = y + size;
  718.     for (i = 0; i < 4; i++) xyz_to_xy (xp[i], yp[i], z, &plot_x[i], &plot_y[i]);
  719.     ps_patch (plot_x, plot_y, 4, r, g, b, outline);
  720. }
  721.  
  722. int sort_on_distance (data, n)
  723. struct DATA1 *data;
  724. int n; {
  725.     /* This function sorts the data array such that points farthest away are plotted first */
  726.     int i;
  727.     double dx, dy, x0, y0, x, y, dr, a, b, c;
  728.  
  729.     x0 = 0.5 * (project_info.xmin + project_info.xmax);
  730.     y0 = 0.5 * (project_info.ymin + project_info.ymax);
  731.     
  732.     dx = 0.5 * (project_info.xmax - project_info.xmin);
  733.     dy = 0.5 * (project_info.ymax - project_info.ymin);
  734.     dr = hypot (dx, dy);
  735.     
  736.     x = x0 - dr * z_project.sin_az;
  737.     y = y0 - dr * z_project.cos_az;
  738.     
  739.     if (z_project.cos_az == 0.0) {
  740.         a = 1.0;
  741.         b = 0.0;
  742.         c = x;
  743.     }
  744.     else {
  745.         a = -tan (z_project.view_azimuth);
  746.         b = -1.0;
  747.         c = y - x * a;
  748.     }
  749.     
  750.     for (i = 0; i < n; i++) data[i].dist = fabs (a * data[i].x + b * data[i].y + c);
  751.     
  752.     qsort ((char *)data, n, sizeof (struct DATA1), compare);
  753. }
  754.  
  755. int compare (point_1, point_2)
  756. struct DATA1 *point_1, *point_2; {
  757.     int first;
  758.     
  759.     if (point_1->dist > point_2->dist)
  760.         return (-1);
  761.     else if (point_1->dist < point_2->dist)
  762.         return (1);
  763.     else {
  764.         first = (point_1->z < point_2->z);
  765.         if (first && z_project.view_elevation >= 0.0)
  766.             return (-1);
  767.         else if (first && z_project.view_elevation < 0.0)
  768.             return (1);
  769.         else
  770.             return (0);
  771.     }
  772. }
  773.