psvelomeca.c ----c语言源代码

/*--------------------------------------------------------------------
 *    The GMT-system:    @(#)psvelomeca.c    2.6  09/20/93
 *
 *    Copyright (c) 1991 by P. Wessel and W. H. F. Smith
 *    See README file for copying and redistribution conditions.
 *--------------------------------------------------------------------*/
/*

psvelomeca will read pairs (or ) from inputfile and
plot symbols on a map. Focal mechanisms, velocity ellipses, strain
crosses, or strain wedges, may be specified, some of which require
additional columns of data.  Only one symbol may be plotted at a time.
PostScript code is written to stdout.


 Authors:     Genevieve Patau and Kurt Feigl
 Date:          20 Sept. 1993
 Version:     2.6
 Roots:       based on psxy.c   

 */

#include "gmt.h"        /* to have gmt environment */

#include "velomeca.h"

#define MECA 1
#define CINE 2
#define ANISO 3
#define WEDGE 4
#define CROSS 5

#define POINTSIZE 0.015
#define TIRET_WIDTH 3
#define veclen(x, y) sqrt((x) * (x) + (y) * (y))

     /* parameters for writing text */
#define ANGLE 0.0
#define FORM 0

int *plot_pen;

main (argc, argv)
int argc;
char **argv; {

    int     i, j, symbol = 0, n, ix = 0, iy = 1, n_files = 0, fno;
    int    n_args;
        int     n_rec = 0;
   
    BOOLEAN    error = FALSE, nofile = TRUE, polygon = FALSE;
        BOOLEAN shade_uncert = FALSE, clip = TRUE;
    BOOLEAN    outline = FALSE;
        BOOLEAN rescale_sigma = FALSE;
    BOOLEAN done, no_size_needed, greenwich;
        BOOLEAN read_cmt = FALSE, read_aki = FALSE, read_planes = FALSE;
        BOOLEAN read_change_position = FALSE, write_change_position = FALSE;
    BOOLEAN read_vector = FALSE, read_ellipse = FALSE, read_rotated_ellipse = FALSE;
        BOOLEAN des_ellipse = TRUE, des_arrow = TRUE;
        BOOLEAN transparence = FALSE;
        BOOLEAN read_anisotropy = FALSE;
        BOOLEAN read_wedge = FALSE;
        BOOLEAN read_cross = FALSE;
   
    double xy[2], xynew[2], west = 0.0, east = 0.0, south = 0.0, north = 0.0;
    double plot_x, plot_y, scale = 0.0, y0;
        double plot_xnew, plot_ynew;
        double vxy[2], vxytip[2], plot_vx, plot_vy;
        double eps1,eps2,spin,spinsig,theta;
    double direction, small_axis, great_axis;
        double sigma_x, sigma_y, corr_xy;
        double confidence = 0., sigma_scale = 1.0, conrad = 1.0;
    double v_width = 0.01, h_length = 0.12, h_width = 0.03, vector_shape = 0.4;
        double wedge_amp = 1.e7;
    double t11 = 1.0, t12 = 0.0, t21 = 0.0, t22 = 1.0;
        double delaz;
        double hl,hw,vw;

        char station_name[20];
        char event_title[512];

    char line[512], symbol_type, col[12][40], EOL_flag = '>';
   
    FILE *fp = NULL;
   
    struct PEN pen, epen;
    struct FILL fill, efill;     /* efill is for uncertainty wedge */
        struct nodal_plane NP1, NP2;
        st_me meca;

        double fault;
        int justify, fontsize;
        double zero_360();
        double ps_mechanism();
        double ps_meca();
        double size;
        double computed_dip2();
        double computed_rake2();

    argc = gmt_begin (argc, argv);
   
    gmt_init_pen (&pen, 1);
    gmt_init_pen (&epen, 1);
    gmt_init_fill (&fill, gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
    gmt_init_fill (&efill, 255-gmtdefs.basemap_frame_rgb[0], 255-gmtdefs.basemap_frame_rgb[1], 255-gmtdefs.basemap_frame_rgb[2]);
   
   
    /* Check and interpret the command line arguments */
   
    for (i = 1; !error && i < argc; i++) {
        if (argv[i][0] == '-') {
            switch(argv[i][1]) {
       
                /* Common parameters */
           
                case 'B':
                case 'H':
                case 'J':
                case 'K':
                case 'O':
                case 'P':
                case 'R':
                case 'U':
                case 'V':
                case 'X':
                case 'x':
                case 'Y':
                case 'y':
                case '#':
                case ':':
                case '\0':
                    error += get_common_args (argv[i], &west, &east, &south, &north);
                    break;
               
                /* Supplemental parameters */
           
                                case 'A':       /* Change size of arrow head */
                                        sscanf(&argv[i][3], "%lf/%lf/%lf", &v_width, &h_length, &h_width);
                                        break;
                                case 'C':       /* Change position */
                                        read_change_position = TRUE;
                                        break;
                                case 'D':       /* Rescale Sigmas */
                                        rescale_sigma = TRUE;
                    sscanf (&argv[i][2], "%lf",&sigma_scale);
                                        break;
                case 'E':        /* Set Gray shade for uncertainty polygon - currently only works for Sw */
                    gmt_getfill (&argv[i][2], &efill);
                    shade_uncert = TRUE;
                    break;
                case 'F':
                    sscanf (&argv[i][2], "%d/%d/%d",&gmtdefs.basemap_frame_rgb[0],
                        &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2]);
                    break;
                case 'G':        /* Set Gray shade for polygon */
                    gmt_getfill (&argv[i][2], &fill);
                    polygon = TRUE;
                    break;
                case 'L':        /* Draw the outline */
                    outline = TRUE;
                    break;
                case 'N':        /* Do not clip at boundary*/
                    clip = FALSE;
                    break;
                case 'S':        /* Get symbol [and size] */
                                        symbol_type = argv[i][2];
                                        if(symbol_type == 'c' || symbol_type == 'a' || symbol_type == 'p')
                                            sscanf(&argv[i][3], "%lf", &scale);
                                        if(symbol_type == 'e' || symbol_type == 'r') {
                                            sscanf(&argv[i][3], "%lf/%lf/%d", &scale, &confidence, &fontsize);
                        /* confidence scaling */
                        conrad = sqrt( -2.0 * log(1.0 - confidence));
                      }
                                        if(symbol_type == 'n' || symbol_type == 'x' )
                                            sscanf(&argv[i][3], "%lf", &scale);
                                        if(symbol_type == 'w')
                                            sscanf(&argv[i][3], "%lf/%lf", &scale, &wedge_amp);

                    switch (symbol_type) {
                                                case 'c':
                                                        symbol = MECA;
                                                        read_cmt = TRUE;
                                                        break;
                                                case 'a':
                                                        symbol = MECA;
                                                        read_aki = TRUE;
                                                        break;
                                                case 'p':
                                                        symbol = MECA;
                                                        read_planes = TRUE;
                                                        break;
                                                case 'e':
                                                        symbol = CINE;
                                                        read_ellipse = TRUE;
                                                        break;
                                                case 'r':
                                                        symbol = CINE;
                                                        read_rotated_ellipse = TRUE;
                                                        break;
                                                case 'n':
                                                        symbol = ANISO;
                                                        read_anisotropy = TRUE;
                                                        break;
                                                case 'w':
                                                        symbol = WEDGE;
                                                        read_wedge = TRUE;
                                                        break;
                                                case 'x':
                                                        symbol = CROSS;
                                                        read_cross = TRUE;
                                                        break;
                        default:
                            error = TRUE;
                            break;
                    }
                    break;
                                case 'T':
                                        transparence = TRUE;
                    break;
                case 'W':        /* Set line attributes */
                    gmt_getpen (&argv[i][2], &pen);
                    break;
                   
                /* Illegal options */
           
                default:        /* Options not recognized */
                    error = TRUE;
                    break;
            }
        }
        else
            n_files++;
    }
   
    /* Check that the options selected are mutually consistant */
   
    no_size_needed = (read_ellipse || read_rotated_ellipse || read_cmt || read_planes || read_aki || read_anisotropy || read_cross || read_wedge );
    error += check_region (west, east, south, north);
    error += check_rgb (pen.r, pen.g, pen.b)
          + check_rgb (fill.r, fill.g, fill.b)
          + check_rgb (efill.r, efill.g, efill.b)
          + check_rgb (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
        /* Only one allowed */
    if ((read_ellipse + read_rotated_ellipse + read_cmt +  read_aki + read_planes + read_anisotropy + read_cross + read_wedge ) > 1) error = TRUE;   
    if (!no_size_needed && (symbol > 1 && scale <= 0.0)) error = TRUE;
        if (rescale_sigma && ! (read_ellipse || read_wedge)) error = TRUE;

    if (argc == 1 || error) {    /* Display usage */
        fprintf (stderr,"psvelomeca plots symbols on maps\n\n");
        fprintf (stderr,"usage: psvelomeca -J -R [-B]\n");
        fprintf (stderr, "    [-F] [-G] [-H] [-K] [-L] [-M] [-N] [-O]\n");
        fprintf (stderr, "    [-P] [-S] [-U[/dx/dy/][]] [-V] [-W] [-c]\n");
        fprintf (stderr, "    [-X] [-Y] [-:]\n\n");
       
        if (gmt_quick) exit(-1);
       
        fprintf (stderr, "    is one or more files.  If no, read standard input\n");
        explain_option ('j');
        explain_option ('R');
        explain_option ('b');
                fprintf (stderr, "      -A Change the size of arrow head; specify arrow_width, head_length, head_width;");
                fprintf (stderr, "         default is 0.01/0.12/0.03;");
                fprintf (stderr, "      -D Multiply uncertainties by sigscale. (Se and Sw only) \n");
        fprintf (stderr, "      -E Set color used for uncertainty wedges in -Sw option.\n");
        fprintf (stderr, "      -F Set color used for Frame and anotation [%d/%d/%d]\n",
            gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
        fprintf (stderr, "    -G Specify color (for symbols/polygons) or pattern (for polygons). fill can be either\n");
        fprintf (stderr, "       1) (each 0-255) for color or (0-255) for gray-shade [0].\n");
        fprintf (stderr, "       2) p[or P]/ for predefined patterns (0-31).\n");
        explain_option ('H');
        explain_option ('K');
        fprintf (stderr, "    -L draw line or symbol outline using the current pen (see -W).\n");
        fprintf (stderr, "    -N Do Not skip/clip symbols that fall outside map border [Default will ignore those outside]\n");
        explain_option ('O');
        explain_option ('P');
        fprintf (stderr, "    -S to select symbol type and scale. Choose between\n");
        fprintf (stderr, "       focal mechanism from CMT or with AKI & RICHARD's convention\n");
        fprintf (stderr, "       or velocity arrows with error ellipses\n");
                fprintf (stderr, "         (c) Focal mechanisms in CMT convention\n");
                fprintf (stderr, "             X, Y, strike1, dip1, rake1, strike2, dip2, rake2, moment, newX, newY, event_title\n");
                fprintf (stderr, "             with moment in 2 columns : mantiss and exponent corresponding to seismic moment in dynes-cm\n");
                fprintf (stderr, "         (a) Focal mechanism in AKI & RICHARD's convention:\n");
                fprintf (stderr, "             X, Y, strike, dip, rake, mag, newX, newY, event_title\n");
                fprintf (stderr, "         (p) Focal mechanism defined with\n");
                fprintf (stderr, "             X, Y, strike1, dip1, strike2, fault, mag, newX, newY, event_title\n");
                fprintf (stderr, "             fault = -1/+1 for a normal/inverse fault\n");
                fprintf (stderr, "         (e) Velocity ellipses: in X,Y,Vx,Vy,SigX,SigY,CorXY,name format \n");
                fprintf (stderr, "         (r) Velocity ellipses: in X,Y,Vx,Vy,a,b,theta,name format \n");
                fprintf (stderr, "         (n) Anisotropy : in X,Y,Vx,Vy \n");
                fprintf (stderr, "         (w) Rotational wedges: in X,Y,Spin,Spinsig \n");
                fprintf (stderr, "         (x) Strain crosses : in X,Y,Eps1,Eps2,Theta \n");
                fprintf (stderr, "      -T draw nodal planes and circumference only to provide a transparent beach ball. \n");
        explain_option ('U');
        explain_option ('V');
        fprintf (stderr, "    -W sets pen attributes [width = %d, color = (%d/%d/%d), texture = solid line].\n",
            pen.width, pen.r, pen.g, pen.b);
        explain_option ('X');
        explain_option ('c');
        explain_option (':');
        explain_option ('.');
        exit(-1);
    }

    if (n_files > 0)
        nofile = FALSE;
    else
        n_files = 1;
    n_args = (argc > 1) ? argc : 2;
   
    greenwich = (west < 0.0 || east <= 0.0);
   
    map_setup (west, east, south, north);

    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, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
    echo_command (argc, argv);
    if (gmtdefs.unix_time && !gmtdefs.overlay) timestamp (argc, argv);
   
    ps_setline (pen.width);
    if (pen.texture) ps_setdash (pen.texture, pen.offset);
    ps_setpaint (pen.r, pen.g, pen.b);
   
    if (clip) map_clip_on (-1, -1, -1, 3);
   
    geo_to_xy (west, 0.0, &plot_x, &y0);    /* Zero level for bars */
   
    ix = (gmtdefs.xy_toggle);    iy = 1 - ix;
    done = FALSE;
    for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over all input files */
        if (!nofile && argv[fno][0] == '-') continue;
        if (nofile) {
            fp = stdin;
            done = TRUE;
        }
        else if ((fp = fopen (argv[fno], "r")) == NULL) {
            fprintf (stderr, "psvelomeca: Cannot open file %s\n", argv[fno]);
            continue;
        }

        if (!nofile && gmtdefs.verbose) {
            fprintf (stderr, "psvelomeca: Working on file %s\n", argv[fno]);
            sprintf (line, "File: %s\0", argv[fno]);
            ps_comment (line);
        }
        if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
       
        if (gmtdefs.verbose && (read_ellipse || read_rotated_ellipse))   
          fprintf (stderr, "psvelomeca: 2-D confidence interval and scaling factor %f %f\n",confidence, conrad);

        while (fgets (line, 512, fp)) {
                        n_rec++;
                        if (read_ellipse || read_rotated_ellipse) {
              sscanf (line, "%s %s %s %s %s %s %s %[^\n]\n",
                  col[0], col[1], col[2], col[3], col[4], col[5], col[6], station_name);
              if (strlen(station_name) <= 0) sprintf(station_name,"\n");
            }
                        else if (read_cmt) {
                          sscanf (line, "%s %s %s %s %s %s %s %s %s %s %s %s %[^\n]\n",
                                  col[0], col[1], col[2], col[3], col[4], col[5], col[6], col[7], col[8], col[9], col[10], col[11], event_title);
              if (strlen(event_title) <= 0) sprintf(event_title,"\n");
            }
                        else if (read_aki) {
                          sscanf (line, "%s %s %s %s %s %s %s %s %[^\n]\n",
                                  col[0], col[1], col[2], col[3], col[4], col[5], col[6], col[7], event_title);
              if (strlen(event_title) <= 0) sprintf(event_title,"\n");
            }
                        else {
                          sscanf (line, "%s %s %s %s %s %s %s %s %s %[^\n]\n",
                  col[0], col[1], col[2], col[3], col[4], col[5], col[6], col[7], col[8], event_title);
              if (strlen(event_title) <= 0) sprintf(event_title,"\n");
            }

            xy[ix] = atof (col[0]);
            xy[iy] = atof (col[1]);
           
                        if (read_cmt) {
                                meca.NP1.str = atof(col[2]);
                                meca.NP1.dip = atof(col[3]);
                                meca.NP1.rake = atof(col[4]);
                                meca.NP2.str = atof(col[5]);
                                meca.NP2.dip = atof(col[6]);
                                meca.NP2.rake = atof(col[7]);
                                meca.moment.mant = atof(col[8]);
                                meca.moment.exponent = atoi(col[9]);
                                if(meca.moment.exponent == 0)
                                    meca.magms = atof(col[8]);

                                if(read_change_position) {
                                    xynew[ix] = atof(col[10]);
                                    xynew[iy] = atof(col[11]);
                                    if(fabs(xynew[ix]) > EPSIL || fabs(xynew[iy]) > EPSIL)
                                        write_change_position = TRUE;
                                    else
                                        write_change_position = FALSE;
                                }
                        }
            else if(read_aki) {
                                NP1.str = atof(col[2]);
                                NP1.dip = atof(col[3]);
                                NP1.rake = atof(col[4]);
                                meca.magms = atof(col[5]);   

                                meca.moment.exponent = 0;
                                define_second_plane(NP1, &NP2);

                                meca.NP1.str = NP1.str;
                                meca.NP1.dip = NP1.dip;
                                meca.NP1.rake = NP1.rake;
                                meca.NP2.str = NP2.str;
                                meca.NP2.dip = NP2.dip;
                                meca.NP2.rake = NP2.rake;

                                if(read_change_position) {
                                    xynew[ix] = atof(col[6]);
                                    xynew[iy] = atof(col[7]);
                                    if(fabs(xynew[ix]) > EPSIL || fabs(xynew[iy]) > EPSIL)
                                        write_change_position = TRUE;
                                    else
                                        write_change_position = FALSE;
                                }
                        }
                        else if(read_planes) {
                                meca.NP1.str = atof(col[2]);
                                meca.NP1.dip = atof(col[3]);
                                meca.NP2.str = atof(col[4]);
                                fault = atof(col[5]);
                                meca.magms = atof(col[6]);

                                meca.moment.exponent = 0;
                                meca.NP2.dip = computed_dip2(meca.NP1.str, meca.NP1.dip, meca.NP2.str);
                                meca.NP2.rake = computed_rake2(meca.NP1.str, meca.NP1.dip, meca.NP2.str, meca.NP2.dip, fault);
                                meca.NP2.rake = computed_rake2(meca.NP1.str, meca.NP1.dip, meca.NP2.str, meca.NP2.dip, fault);
                                meca.NP1.rake = computed_rake2(meca.NP2.str, meca.NP2.dip, meca.NP1.str, meca.NP1.dip, fault);

                                if(read_change_position) {
                                    xynew[ix] = atof(col[7]);
                                    xynew[iy] = atof(col[8]);
                                    if(fabs(xynew[ix]) > EPSIL || fabs(xynew[iy]) > EPSIL)
                                        write_change_position = TRUE;
                                    else
                                        write_change_position = FALSE;
                                }
                        }

                        else if(read_ellipse) {
                                    vxy[ix] = atof(col[2]);
                                    vxy[iy] = atof(col[3]);
                                    sigma_x = atof(col[4]);
                                    sigma_y = atof(col[5]);
                                    corr_xy = atof(col[6]);
                                    /* rescale uncertainties if necessary */
                                    if (rescale_sigma) {
                      sigma_x = sigma_scale * sigma_x;
                      sigma_y = sigma_scale * sigma_y;
                         }
                                    if(fabs(sigma_x) < EPSIL && fabs(sigma_y) < EPSIL)
                                        des_ellipse = FALSE;
                                    else {
                                        des_ellipse = TRUE;
                                        ellipse_convert(sigma_x, sigma_y, corr_xy, conrad, &small_axis, &great_axis, &direction);

                                    /* convert to degrees */
                                       direction = direction / D2R;
                                    }
                        }
                        else if(read_rotated_ellipse) {
                                    vxy[ix] = atof(col[2]);
                                    vxy[iy] = atof(col[3]);
                                    great_axis = conrad*atof(col[4]);
                                    small_axis = conrad*atof(col[5]);
                                    direction = atof(col[6]);
                                    if(fabs(great_axis) < EPSIL && fabs(small_axis) < EPSIL)
                                        des_ellipse = FALSE;
                                    else
                                        des_ellipse = TRUE;
                        }
                        else if(read_anisotropy) {
                                    vxy[ix] = atof(col[2]);
                                    vxy[iy] = atof(col[3]);
                        }
                        else if(read_cross) {
                                    eps1 = atof(col[2]);
                                    eps2 = atof(col[3]);
                                    theta = atof(col[4]);
                  }
                        else if(read_wedge) {
                                    spin = atof(col[2]);
                                    spinsig = atof(col[3]);
                                    if (rescale_sigma) spinsig = spinsig * sigma_scale;
                        }

            map_outside (xy[0], xy[1]);
            if ( abs (gmt_x_status_new) > 1 || abs (gmt_y_status_new) > 1) continue;

            geo_to_xy (xy[0], xy[1], &plot_x, &plot_y);
           
               
            switch (symbol) {
                                case MECA:
                                        if(write_change_position) {
                                            geo_to_xy(xynew[0], xynew[1], &plot_xnew, &plot_ynew);
                                            ps_circle(plot_x, plot_y, POINTSIZE, fill.r, fill.g, fill.b, 1);
                                            ps_setline(TIRET_WIDTH);
                                            ps_plot(plot_x, plot_y, 3);
                                            ps_plot(plot_xnew, plot_ynew, 2);
                                            ps_setline(1);
                                            plot_x = plot_xnew;
                                            plot_y = plot_ynew;
                                        }
                                        get_trans(xy[0], xy[1], &t11, &t12, &t21, &t22);
                                        delaz = acos(t11 / sqrt(squared(t11) + squared(t21))) / D2R;
                                        meca.NP1.str = zero_360(meca.NP1.str + delaz);
                                        meca.NP2.str = zero_360(meca.NP2.str + delaz);
                                        size = transparence == 1 ? ps_meca(plot_x, plot_y, meca, scale, outline) : ps_mechanism(plot_x, plot_y, meca, scale, fill.r, fill.g, fill.b, outline);
                                        justify = 2;
                                        fontsize = 9;
                                        ps_text(plot_x, plot_y + 1.2 * size, fontsize, event_title, ANGLE, justify, FORM);
                                        ps_plot(plot_x, plot_y + 1.2 * size, 3);
                                        break;
                                case CINE:
                                        des_arrow = veclen((vxy[0]), (vxy[1])) < 1.e-8 ? FALSE : TRUE;
                                        trace_arrow(xy[0], xy[1], vxy[0], vxy[1], scale, &plot_x, &plot_y, &plot_vx, &plot_vy);
                                        get_trans (xy[0], xy[1], &t11, &t12, &t21, &t22);
                                        if(des_ellipse)
                                            paint_ellipse (plot_vx, plot_vy, direction, great_axis, small_axis, scale,
                                            t11,t12,t21,t22, shade_uncert, efill.r, efill.g, efill.b, outline);
                                        if(des_arrow) {
                      /* verify that arrow is not ridiculously small */
                      if (veclen (plot_x-plot_vx,plot_y-plot_vy) <= 1.5 * h_length) {
                        hl = veclen (plot_x-plot_vx,plot_y-plot_vy) * 0.6;
                        hw = hl * h_width/h_length;
                        vw = hl * v_width/h_length;
                        if (vw < 2./(double)gmtdefs.dpi) vw = 2./(double)gmtdefs.dpi;
                      }
                      else {
                        hw = h_width;
                        hl = h_length;
                        vw = v_width;
                      }
                      ps_vector(plot_x, plot_y, plot_vx, plot_vy, vw, hl, hw, vector_shape,
                            fill.r, fill.g, fill.b, outline);
                      justify = plot_vx - plot_x > 0. ? 7 : 5;
                                          if(fontsize > 0) {
                                ps_text(plot_x + (6 - justify) / 25.4 , plot_y, fontsize,
                                    station_name, ANGLE, justify,FORM); /* 1 inch = 2.54 cm */
                                          }
                                        }
                                        else {
                                            ps_plot(plot_x, plot_y, 3);
                                            ps_circle(plot_x, plot_y, POINTSIZE, fill.r, fill.g, fill.b, TRUE);
                                            justify = 10;
                                            if(fontsize > 0) {
                                                 ps_text(plot_x, plot_y - 1. / 25.4, fontsize, station_name, ANGLE, justify, FORM);
                                            }
                                            /*  1 inch = 2.54 cm */
                                        }
                                        i = 0;
                                        while(col[7][i] != '\0') {
                                            col[7][i] = ' ';
                                            i++;
                                        }
                                        i = 0;
                                        while(station_name[i] != '\0') {
                                            station_name[i] = ' ';
                                            i++;
                                        }
                                        break;
                                case ANISO:
                                        trace_arrow(xy[0], xy[1], vxy[0], vxy[1], scale, &plot_x, &plot_y, &plot_vx, &plot_vy);
                                        ps_plot(plot_x, plot_y, 3);
                                        ps_plot(plot_vx, plot_vy,2);
                                        break;
                                case CROSS:
                                        vector_shape = 0.1; /* triangular arrowheads */
                                        trace_cross(xy[0],xy[1],eps1,eps2,theta,scale,v_width,h_length,h_width,vector_shape,outline,pen);
                                        break;
                                case WEDGE:
                                        sprintf (line, "begin wedge number %i",n_rec);
                    ps_comment (line);
                                        geo_to_xy(xy[0], xy[1], &plot_x, &plot_y);
                                        get_trans (xy[0], xy[1], &t11, &t12, &t21, &t22);
                                        paint_wedge (plot_x, plot_y, spin, spinsig, scale, wedge_amp, t11,t12,t21,t22,
                             polygon, fill.r, fill.g, fill.b,
                             shade_uncert, efill.r, efill.g, efill.b,outline);
                                        break;
            }
        }
    if (fp != stdin) fclose (fp);
    }
   
    if (gmtdefs.verbose)
            fprintf (stderr, "psvelomeca: Number of records read: %i\n", n_rec);

    if (gmtdefs.verbose && rescale_sigma)
            fprintf (stderr, "psvelomeca: Rescaling uncertainties by a factor of %f\n", sigma_scale);


    if (clip) map_clip_off ();
   
    if (pen.texture) ps_setdash (CNULL, 0);
    if (frame_info.plot) {
        ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
        map_basemap ();
        ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
    }
    ps_plotend (gmtdefs.last_page);
   
    gmt_end (argc, argv);
}