home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_plot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-09  |  94.4 KB  |  3,086 lines
  1.  
  2. /*--------------------------------------------------------------------
  3.  *    The GMT-system:    @(#)gmt_plot.c    2.58  09 Aug 1995
  4.  *
  5.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  6.  *    See README file for copying and redistribution conditions.
  7.  *--------------------------------------------------------------------*/
  8. /*
  9.  *
  10.  *            G M T _ P L O T . C
  11.  *
  12.  *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  13.  * gmt_plot.c contains code related to plotting maps
  14.  *
  15.  * Author:    Paul Wessel
  16.  * Date:    27-JUL-1992
  17.  * Version:    2.1
  18.  *
  19.  *
  20.  * Functions include:
  21.  *
  22.  *    basemap_3D :        Plots 3-D basemap
  23.  *    color_image :        Calls the desired colorimage operator
  24.  *    echo_command :        Puts the command line into the PostScript file as comments
  25.  *    geoplot :        As ps_plot, but using lon/lat directly
  26.  *    gmt_plot_line :        Same for polygons
  27.  *    get_angle :        Sub function to get annotation angles
  28.  *    get_anot_label :    Construct degree/minute label
  29.  *    lambert_map_boundary :    Plot basemap for Lambert conformal conic projection
  30.  *    linear_map_boundary :    Plot basemap for Linear projections
  31.  *    logx_grid :        Draw log10 x grid lines
  32.  *    logy_grid :        Draw log10 y grid lines
  33.  *    prepare_label :        Gets angle and justification for frame anotations
  34.  *    map_anotate :        Annotate basemaps
  35.  *    map_basemap :        Generic basemap function
  36.  *    map_boundary :        Draw the maps boundary
  37.  *    map_clip_off :        Deactivate map region clip path
  38.  *    map_clip_on :        Activate map region clip path
  39.  *    map_gridcross :        Draw grid crosses on maps
  40.  *    map_gridlines :        Draw gridlines on maps
  41.  *    map_latline :        Draw a latitude line
  42.  *    map_lattick :        Draw a latitude tick marck
  43.  *    map_lonline :        Draw a longitude line
  44.  *    map_lontick :        Draw a longitude tick mark
  45.  *    map_symbol :        Plot map annotation
  46.  *    map_symbol_ew :          for east/west sides
  47.  *    map_symbol_ns :          for south/north sides
  48.  *    map_tick :        Draw the ticks
  49.  *    map_tickmarks :        Plot tickmarks on maps
  50.  *    merc_map_boundary :    Plot basemap for Mercator projection
  51.  *    ellipse_map_boundary :    Plot basemap for Mollweide and Hammer-Aitoff projections
  52.  *    oblmrc_map_boundary :    Plot basemap for Oblique Mercator projection
  53.  *    polar_map_boundary :    Plot basemap for Polar stereographic projection
  54.  *    rect_map_boundary :    Plot plain basemap for projections with rectangular boundaries
  55.  *    basic_map_boundary :    Plot plain basemap for most projections
  56.  *    timestamp :        Plot UNIZ time stamp with optional string
  57.  *    wesn_map_boundary :    Plot plain basemap for projections with geographical boundaries
  58.  *    x_axis :        Draw x axis
  59.  *    xyz_axis3D :        Draw 3-D axis
  60.  *    y_axis :        Draw y axis
  61.  *    polar_adjust :        Adjust label justification for polar projection
  62.  *    
  63.  
  64.  */
  65.  
  66. #include "gmt.h"
  67.  
  68. struct XINGS {
  69.     double xx[2], yy[2];    /* Cartesian coordinates of intersection with map boundary */
  70.     double angle[2];    /* Angles of intersection */
  71.     int sides[2];        /* Side id of intersection */
  72.     int nx;            /* Number of intersections (1 or 2) */
  73. };
  74.  
  75.  
  76. /*    LINEAR PROJECTION MAP FUNCTIONS    */
  77.  
  78. int linear_map_boundary (w, e, s, n)
  79. double w, e, s, n; {
  80.     double x1, x2, y1, y2, x_length, y_length;
  81.     
  82.     geo_to_xy (w, s, &x1, &y1);
  83.     geo_to_xy (e, n, &x2, &y2);
  84.     if (!project_info.xyz_pos[0]) d_swap (x1, x2);
  85.     if (!project_info.xyz_pos[1]) d_swap (y1, y2);
  86.     x_length = fabs (x2 - x1);
  87.     y_length = fabs (y2 - y1);
  88.     
  89.     if (frame_info.side[3])    /* West or left y-axis */
  90.         y_axis (x1, y1, y_length, s, n, frame_info.anot_int[1], frame_info.frame_int[1],
  91.                 frame_info.label[1], project_info.xyz_projection[1], frame_info.anot_type[1], TRUE, frame_info.side[3]-1);
  92.     
  93.     if (frame_info.side[1])    /* East or right y-axis */
  94.         y_axis (x2, y1, y_length, s, n, frame_info.anot_int[1], frame_info.frame_int[1],
  95.             frame_info.label[1], project_info.xyz_projection[1], frame_info.anot_type[1], FALSE, frame_info.side[1]-1);
  96.     
  97.     if (frame_info.side[0])    /* South or lower x-axis */
  98.         x_axis (x1, y1, x_length, w, e, frame_info.anot_int[0], frame_info.frame_int[0],
  99.             frame_info.label[0], project_info.xyz_projection[0], frame_info.anot_type[0], TRUE, frame_info.side[0]-1);
  100.     
  101.     if (frame_info.side[2])    /* North or upper x-axis */
  102.         x_axis (x1, y2, x_length, w, e, frame_info.anot_int[0], frame_info.frame_int[0],
  103.             frame_info.label[0], project_info.xyz_projection[0], frame_info.anot_type[0], FALSE, frame_info.side[2]-1);
  104. }
  105.  
  106. int x_axis (x0, y0, length, val0, val1, anotation_int, tickmark_int, label, axistype, anottype, below, anotate)
  107. double x0, y0, length, val0, val1, anotation_int, tickmark_int;
  108. char *label;
  109. int axistype, anottype, below, anotate; {
  110.     int i, i_a, i_f, justify, n_anotations = 0, n_tickmarks = 0, test;
  111.     
  112.     BOOLEAN left = FALSE, do_anot, do_tick;
  113.     
  114.     double axis_scale, val, v0, v1, anot_off, label_off, start_val_a, start_val_f, end_val;
  115.     double tvals_a[9], tvals_f[9], sign, dy, tmp, xx, len, start_log_a, start_log_f;
  116.     
  117.     char annotation[80], format[20];
  118.     
  119.     ps_setfont (gmtdefs.anot_font);
  120.     justify = (below) ? 10 : 2;
  121.     sign = (below) ? -1.0 : 1.0;
  122.     dy = sign * gmtdefs.tick_length;
  123.     len = (gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0;
  124.     left = (anotation_int < 0.0);
  125.     anotation_int = fabs (anotation_int);
  126.     tickmark_int = fabs (tickmark_int);
  127.     
  128.     do_anot = (anotation_int > 0.0);
  129.     do_tick = (tickmark_int > 0.0);
  130.         
  131.     /* Find number of decimals needed, if any */
  132.     
  133.     get_format (anotation_int, format);
  134.  
  135.     anot_off = sign * (len + gmtdefs.anot_offset);
  136.     label_off = sign * (len + 2.5 * gmtdefs.anot_offset + (gmtdefs.anot_font_size / gmt_ppu[gmtdefs.measure_unit]) * font_height[gmtdefs.anot_font]);
  137.     
  138.     /* Ready to draw axis */
  139.     
  140.     if (below) ps_comment ("Start of lower x-axis"); else ps_comment ("Start of upper x-axis");
  141.     ps_transrotate (x0, y0, 0.0);
  142.     ps_setline (gmtdefs.frame_pen);
  143.     ps_plot (0.0, 0.0, 3);
  144.     ps_plot (length, 0.0, -2);
  145.     ps_setline (gmtdefs.tick_pen);
  146.     
  147.     i_a = i_f = 0;
  148.     switch (axistype) {
  149.         case POW:    /* Anotate in pow(x) */
  150.             (*x_forward) (val0, &v0);
  151.             (*x_forward) (val1, &v1);
  152.             if (anottype == 2) {
  153.                 val = (anotation_int == 0.0) ? 0.0 : floor (v0 / anotation_int) * anotation_int;
  154.                 if (fabs (val - v0) > SMALL) val += anotation_int;
  155.                 start_val_a = val;    end_val = v1;
  156.                 val = (tickmark_int == 0.0) ? 0.0 : floor (v0 / tickmark_int) * tickmark_int;
  157.                 if (fabs (val - v0) > SMALL) val += tickmark_int;
  158.                 start_val_f = val;
  159.             }
  160.             else {
  161.                 val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  162.                 if (fabs (val - val0) > SMALL) val += anotation_int;
  163.                 start_val_a = val;    end_val = val1;
  164.                 val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  165.                 if (fabs (val - val0) > SMALL) val += tickmark_int;
  166.                 start_val_f = val;
  167.             }
  168.             break;
  169.         case LOG10:    /* Anotate in d_log10 (x) */
  170.             v0 = d_log10 (val0);
  171.             v1 = d_log10 (val1);
  172.             val = pow (10.0, floor (v0));
  173.             test = rint (anotation_int) - 1;
  174.             if (test < 0 || test > 2) test = 0;
  175.             if (test == 0) {
  176.                 n_anotations = 1;
  177.                 tvals_a[0] = 10.0;
  178.             }
  179.             else if (test == 1) {
  180.                 tvals_a[0] = 1.0;
  181.                 tvals_a[1] = 2.0;
  182.                 tvals_a[2] = 5.0;
  183.                 n_anotations = 3;
  184.             }
  185.             else if (test == 2) {
  186.                 n_anotations = 9;
  187.                 for (i = 0; i < n_anotations; i++) tvals_a[i] = i + 1;
  188.             }
  189.             test = rint (tickmark_int) - 1;
  190.             if (test < 0 || test > 2) test = 0;
  191.             if (test == 0) {
  192.                 n_tickmarks = 1;
  193.                 tvals_f[0] = 10.0;
  194.             }
  195.             if (test == 1) {
  196.                 tvals_f[0] = 1.0;
  197.                 tvals_f[1] = 2.0;
  198.                 tvals_f[2] = 5.0;
  199.                 n_tickmarks = 3;
  200.             }
  201.             else if (test == 2) {
  202.                 n_tickmarks = 9;
  203.                 for (i = 0; i < n_tickmarks; i++) tvals_f[i] = i + 1;
  204.             }
  205.             i_a = 0;
  206.             start_log_a = val = pow (10.0, floor (v0));
  207.             while ((val0 - val) > SMALL) {
  208.                 if (i_a < n_anotations)
  209.                     val = start_log_a * tvals_a[i_a];
  210.                 else {
  211.                     val = (start_log_a *= 10.0);
  212.                     i_a = 0;
  213.                 }
  214.                 i_a++;
  215.             }
  216.             i_a--;
  217.             start_val_a = val;
  218.             i_f = 0;
  219.             start_log_f = val = pow (10.0, floor (v0));
  220.             while ((val0 - val) > SMALL) {
  221.                 if (i_f < n_tickmarks)
  222.                     val = start_log_f * tvals_f[i_f];
  223.                 else {
  224.                     val = (start_log_f *= 10.0);
  225.                     i_f = 0;
  226.                 }
  227.                 i_f++;
  228.             }
  229.             i_f--;
  230.             start_val_f = val;
  231.             end_val = val1;
  232.             break;
  233.         case LINEAR:
  234.             val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  235.             if (fabs (val - val0) > SMALL) val += anotation_int;
  236.             start_val_a = val;    end_val = val1;
  237.             val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  238.             if (fabs (val - val0) > SMALL) val += tickmark_int;
  239.             start_val_f = val;
  240.             (*x_forward) (val0, &v0);
  241.             (*x_forward) (val1, &v1);
  242.             break;
  243.     }
  244.     
  245.     axis_scale = length / (v1 - v0);
  246.     
  247.     /* Do anotations with tickmarks */
  248.     
  249.     val = (anotation_int == 0.0) ? end_val + 1.0 : start_val_a;
  250.     
  251.     while (do_anot && val <= end_val) {
  252.     
  253.         i_a++;
  254.  
  255.         switch (anottype) {
  256.             case 0:
  257.                 (*x_forward) (val, &tmp);
  258.                 xx = (tmp - v0) * axis_scale;
  259.                 sprintf (annotation, format, val);
  260.                 break;
  261.             case 1:
  262.                 sprintf (annotation, "%d\0", (int) d_log10 (val));
  263.                 xx = (d_log10 (val) - v0) * axis_scale;
  264.                 break;
  265.             case 2:
  266.                 if (axistype == POW) {
  267.                     xx = (val - v0) * axis_scale;
  268.                     (*x_inverse) (&tmp, val);
  269.                     sprintf (annotation, format, tmp);
  270.                 }
  271.                 else {
  272.                     xx = (d_log10 (val) - v0) * axis_scale;
  273.                     sprintf (annotation, "10@+%d@+\0", (int) d_log10 (val));
  274.                 }
  275.                 break;
  276.         }
  277.                 
  278.         if (left) xx = length - xx;
  279.         
  280.         ps_plot (xx, 0.0, 3);
  281.         ps_plot (xx, dy, -2);
  282.         if (anotate) ps_text (xx, anot_off, gmtdefs.anot_font_size, annotation, 0.0, justify, 0);
  283.         
  284.         if (axistype == LOG10) {
  285.             if (i_a < n_anotations)
  286.                 val = start_log_a * tvals_a[i_a];
  287.             else {
  288.                 val = (start_log_a *= 10.0);
  289.                 i_a = 0;
  290.             }
  291.         }
  292.         else
  293.             val = start_val_a + i_a * anotation_int;
  294.             
  295.     }
  296.  
  297.     /* Now do frame tickmarks */
  298.     
  299.     dy *= 0.5;
  300.     
  301.     val = (tickmark_int == 0.0) ? end_val + 1.0 : start_val_f;
  302.  
  303.     while (do_tick && val <= end_val) {
  304.     
  305.         i_f++;
  306.         
  307.         switch (anottype) {
  308.             case 0:
  309.                 (*x_forward) (val, &tmp);
  310.                 xx = (tmp - v0) * axis_scale;
  311.                 break;
  312.             case 1:
  313.                 xx = (d_log10 (val) - v0) * axis_scale;
  314.                 break;
  315.             case 2:
  316.                 xx =  ((axistype == POW) ? (val - v0) : (d_log10 (val) - v0)) * axis_scale;
  317.                 break;
  318.         }
  319.                 
  320.         if (left) xx = length - xx;
  321.         
  322.         ps_plot (xx, 0.0, 3);
  323.         ps_plot (xx, dy, -2);
  324.         
  325.         if (axistype == LOG10) {
  326.             if (i_f < n_tickmarks)
  327.                 val = start_log_f * tvals_f[i_f];
  328.             else {
  329.                 val = start_log_f *= 10.0;
  330.                 i_f = 0;
  331.             }
  332.         }
  333.         else
  334.             val = start_val_f + i_f * tickmark_int;
  335.     }
  336.  
  337.     /* Finally do label */
  338.     
  339.     ps_setfont (gmtdefs.label_font);
  340.     if (label && anotate) ps_text (0.5 * length, label_off, gmtdefs.label_font_size, label, 0.0, justify, 0);
  341.     ps_rotatetrans  (-x0, -y0, 0.0);
  342.     ps_comment ("End of x-axis");
  343. }
  344.  
  345. int y_axis (x0, y0, length, val0, val1, anotation_int, tickmark_int, label, axistype, anottype, left_side, anotate)
  346. double x0, y0, length, val0, val1, anotation_int, tickmark_int;
  347. char *label;
  348. int axistype, anottype, left_side, anotate; {
  349.     int i, i_a, i_f, anot_justify, label_justify, n_anotations = 0, n_tickmarks = 0, ndec, test;
  350.     
  351.     BOOLEAN left = FALSE, do_anot, do_tick;
  352.     
  353.     double axis_scale, val, v0, v1, anot_off, label_off, start_val_a, start_val_f, end_val;
  354.     double tvals_a[9], tvals_f[9], sign, len, dy, tmp, xx, off, angle, start_log_a, start_log_f;
  355.     
  356.     char annotation[80], text_l[80], text_u[80], format[20];
  357.     
  358.     ps_setfont (gmtdefs.anot_font);
  359.     sign = (left_side) ? 1.0 : -1.0;
  360.     len = (gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0;
  361.     dy = sign * gmtdefs.tick_length;
  362.     left = (anotation_int < 0.0);
  363.     anotation_int = fabs (anotation_int);
  364.     tickmark_int = fabs (tickmark_int);
  365.     
  366.     do_anot = (anotation_int > 0.0);
  367.     do_tick = (tickmark_int > 0.0);
  368.     
  369.     ndec = get_format (anotation_int, format);
  370.     
  371.     /* Ready to draw axis */
  372.     
  373.     if (left_side) ps_comment ("Start of left y-axis"); else ps_comment ("Start of right y-axis");
  374.     ps_transrotate (x0, y0, 90.0);
  375.     ps_setline (gmtdefs.frame_pen);
  376.     ps_plot (0.0, 0.0, 3);
  377.     ps_plot (length, 0.0, -2);
  378.     ps_setline (gmtdefs.tick_pen);
  379.     
  380.     i_a = i_f = 0;
  381.     switch (axistype) {
  382.         case POW:
  383.             (*y_forward) (val0, &v0);
  384.             (*y_forward) (val1, &v1);
  385.             if (anottype == 2) {
  386.                 val = (anotation_int == 0.0) ? 0.0 : floor (v0 / anotation_int) * anotation_int;
  387.                 if (fabs (val - v0) > SMALL) val += anotation_int;
  388.                 start_val_a = val;    end_val = v1;
  389.                 val = (tickmark_int == 0.0) ? 0.0 : floor (v0 / tickmark_int) * tickmark_int;
  390.                 if (fabs (val - v0) > SMALL) val += tickmark_int;
  391.                 start_val_f = val;
  392.                 sprintf (text_l, "%d\0", (int)floor (v0));
  393.                 sprintf (text_u, "%d\0", (int)ceil (v1));
  394.             }
  395.             else {    /* Anotate in x */
  396.                 val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  397.                 if (fabs (val - val0) > SMALL) val += anotation_int;
  398.                 start_val_a = val;    end_val = val1;
  399.                 val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  400.                 if (fabs (val - val0) > SMALL) val += tickmark_int;
  401.                 start_val_f = val;
  402.                 sprintf (text_l, "%d\0", (int)floor (val0));
  403.                 sprintf (text_u, "%d\0", (int)ceil (val1));
  404.             }
  405.             break;
  406.         case LOG10:
  407.             v0 = d_log10 (val0);
  408.             v1 = d_log10 (val1);
  409.             val = pow (10.0, floor (v0));
  410.             test = rint (anotation_int) - 1;
  411.             if (test < 0 || test > 2) test = 0;
  412.             if (test == 0) {
  413.                 n_anotations = 1;
  414.                 tvals_a[0] = 10.0;
  415.                 if (fabs (val - val0) > SMALL) val *= 10.0;
  416.             }
  417.             else if (test == 1) {
  418.                 tvals_a[0] = 1.0;
  419.                 tvals_a[1] = 2.0;
  420.                 tvals_a[2] = 5.0;
  421.                 n_anotations = 3;
  422.             }
  423.             else if (test == 2) {
  424.                 n_anotations = 9;
  425.                 for (i = 0; i < n_anotations; i++) tvals_a[i] = i + 1;
  426.             }
  427.             test = rint (tickmark_int) - 1;
  428.             if (test < 0 || test > 2) test = 0;
  429.             if (test == 0) {
  430.                 n_tickmarks = 1;
  431.                 tvals_f[0] = 10.0;
  432.             }
  433.             else if (test == 1) {
  434.                 tvals_f[0] = 1.0;
  435.                 tvals_f[1] = 2.0;
  436.                 tvals_f[2] = 5.0;
  437.                 n_tickmarks = 3;
  438.             }
  439.             else if (test == 2) {
  440.                 n_tickmarks = 9;
  441.                 for (i = 0; i < n_tickmarks; i++) tvals_f[i] = i + 1;
  442.             }
  443.             i_a = 0;
  444.             start_log_a = val = pow (10.0, floor (v0));
  445.             while ((val0 - val) > SMALL) {
  446.                 if (i_a < n_anotations)
  447.                     val = start_log_a * tvals_a[i_a];
  448.                 else {
  449.                     val = (start_log_a *= 10.0);
  450.                     i_a = 0;
  451.                 }
  452.                 i_a++;
  453.             }
  454.             i_a--;
  455.             start_val_a = val;
  456.             i_f = 0;
  457.             start_log_f = val = pow (10.0, floor (v0));
  458.             while ((val0 - val) > SMALL) {
  459.                 if (i_f < n_tickmarks)
  460.                     val = start_log_f * tvals_f[i_f];
  461.                 else {
  462.                     val = (start_log_f *= 10.0);
  463.                     i_f = 0;
  464.                 }
  465.                 i_f++;
  466.             }
  467.             i_f--;
  468.             start_val_f = val;
  469.             end_val = val1;
  470.             if (anottype == 2) {    /* 10 ^ pow annotations */
  471.                 sprintf (text_l, "10%d\0", (int)floor (v0));
  472.                 sprintf (text_u, "10%d\0", (int)ceil (v1));
  473.             }
  474.             else {
  475.                 sprintf (text_l, "%d\0", (int)floor ((anottype == 0) ? val0 : v0));
  476.                 sprintf (text_u, "%d\0", (int)ceil ((anottype == 0) ? val1 : v1));
  477.             }
  478.             break;
  479.         case LINEAR:
  480.             v0 = val0;
  481.             v1 = val1;
  482.             val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  483.             if (fabs (val - val0) > SMALL) val += anotation_int;
  484.             start_val_a = val;    end_val = val1;
  485.             val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  486.             if (fabs (val - val0) > SMALL) val += tickmark_int;
  487.             start_val_f = val;
  488.             sprintf (text_l, "%d\0", (int)floor (fabs (val0)));
  489.             sprintf (text_u, "%d\0", (int)ceil (fabs (val1)));
  490.             break;
  491.     }
  492.     
  493.     /* Find offset based on no of digits before and after a period, if any */
  494.     
  495.     off = ((MAX (strlen (text_l), strlen (text_u)) + ndec) * 0.55 + ((ndec > 0) ? 0.3 : 0.0) + ((val0 < 0.0) ? 0.3 : 0.0))
  496.         * gmtdefs.anot_font_size / gmt_ppu[gmtdefs.measure_unit];
  497.  
  498.     axis_scale = length / (v1 - v0);
  499.     
  500.     if (gmtdefs.y_axis_type == 0) {    /* Horizontal anotations */
  501.         if (left_side) {
  502.             anot_off = len + gmtdefs.anot_offset;
  503.             label_off = anot_off + 1.5 * gmtdefs.anot_offset + off;
  504.             anot_justify = 7;
  505.             label_justify = 2;
  506.         }
  507.         else {
  508.             anot_off = -(len + gmtdefs.anot_offset + off);
  509.             label_off = anot_off - 1.5 * gmtdefs.anot_offset;
  510.             anot_justify = 7;
  511.             label_justify = 10;
  512.         }
  513.         angle = -90.0;
  514.     }
  515.     else {
  516.         anot_off = sign * (len + gmtdefs.anot_offset);
  517.         label_off = sign * (len + 2.5 * gmtdefs.anot_offset + (gmtdefs.anot_font_size / gmt_ppu[gmtdefs.measure_unit]) * font_height[gmtdefs.anot_font]);
  518.         anot_justify = label_justify = (left_side) ? 2 : 10;
  519.         angle = 0.0;
  520.     }
  521.     
  522.     /* Do anotations */
  523.     
  524.     val = (anotation_int == 0.0) ? end_val + 1.0 : start_val_a;
  525.     
  526.     while (do_anot && val <= end_val) {
  527.     
  528.         i_a++;
  529.         
  530.         switch (anottype) {
  531.             case 0:
  532.                 sprintf (annotation, format, val);
  533.                 (*y_forward) (val, &tmp);
  534.                 xx = (tmp - v0) * axis_scale;
  535.                 break;
  536.             case 1:
  537.                 sprintf (annotation, "%d\0", (int) d_log10 (val));
  538.                 xx = (d_log10 (val) - v0) * axis_scale;
  539.                 break;
  540.             case 2:
  541.                 if (axistype == POW) {
  542.                     xx = (val - v0) * axis_scale;
  543.                     (*y_inverse) (&tmp, val);
  544.                     sprintf (annotation, format, tmp);
  545.                 }
  546.                 else {
  547.                     xx = (d_log10 (val) - v0) * axis_scale;
  548.                     sprintf (annotation, "10@+%d@+\0", (int) d_log10 (val));
  549.                 }
  550.                 break;
  551.         }
  552.                 
  553.         if (left) xx = length - xx;
  554.         
  555.         ps_plot (xx, 0.0, 3);
  556.         ps_plot (xx, dy, -2);
  557.         if (anotate) ps_text (xx, anot_off, gmtdefs.anot_font_size, annotation, angle, anot_justify, 0);
  558.         
  559.         if (axistype == LOG10) {
  560.             if (i_a < n_anotations)
  561.                 val = start_log_a * tvals_a[i_a];
  562.             else {
  563.                 val = (start_log_a *= 10.0);
  564.                 i_a = 0;
  565.             }
  566.         }
  567.         else
  568.             val = start_val_a + i_a * anotation_int;
  569.             
  570.     }
  571.  
  572.  
  573.     /* Now do frame tickmarks */
  574.     
  575.     dy *= 0.5;
  576.     
  577.     val = (tickmark_int == 0.0) ? end_val + 1.0 : start_val_f;
  578.  
  579.     while (do_tick && val <= end_val) {
  580.     
  581.         i_f++;
  582.         
  583.         switch (anottype) {
  584.             case 0:
  585.                 (*y_forward) (val, &tmp);
  586.                 xx = (tmp - v0) * axis_scale;
  587.                 break;
  588.             case 1:
  589.                 xx = (d_log10 (val) - v0) * axis_scale;
  590.                 break;
  591.             case 2:
  592.                 xx =  ((axistype == POW) ? (val - v0) : (d_log10 (val) - v0)) * axis_scale;
  593.                 break;
  594.         }
  595.                 
  596.         if (left) xx = length - xx;
  597.         
  598.         ps_plot (xx, 0.0, 3);
  599.         ps_plot (xx, dy, -2);
  600.         
  601.         if (axistype == LOG10) {
  602.             if (i_f < n_tickmarks)
  603.                 val = start_log_f * tvals_f[i_f];
  604.             else {
  605.                 val = start_log_f *= 10.0;
  606.                 i_f = 0;
  607.             }
  608.         }
  609.         else
  610.             val = start_val_f + i_f * tickmark_int;
  611.     }
  612.  
  613.     /* Finally do label */
  614.     
  615.     ps_setfont (gmtdefs.label_font);
  616.     if (label && anotate) ps_text (0.5 * length, label_off, gmtdefs.label_font_size, label, 0.0, label_justify, 0);
  617.     ps_rotatetrans  (-x0, -y0, -90.0);
  618.     ps_comment ("End of y-axis");
  619. }
  620.  
  621. int logx_grid (w, e, s, n, dval)
  622. double w, e, s, n, dval; {
  623.     int i, nticks, test;
  624.     double val, v0, end_val, start_log, tvals[9];
  625.     
  626.     test = rint (fabs (dval)) - 1;
  627.     if (test < 0 || test > 2) test = 0;
  628.     if (test == 0) {
  629.         tvals[0] = 1.0;
  630.         nticks = 1;
  631.     }
  632.     if (test == 1) {
  633.         tvals[0] = 1.0;
  634.         tvals[1] = 2.0;
  635.         tvals[2] = 5.0;
  636.         nticks = 3;
  637.     }
  638.     else if (test == 2) {
  639.         nticks = 9;
  640.         for (i = 0; i < nticks; i++) tvals[i] = i + 1;
  641.     }
  642.     
  643.     v0 = d_log10 (w);
  644.     start_log = val = pow (10.0, floor (v0));
  645.     i = 0;
  646.     while ((w - val) > SMALL) {
  647.         if (i < nticks)
  648.             val = start_log * tvals[i];
  649.         else {
  650.             val = (start_log *= 10.0);
  651.             i = 0;
  652.         }
  653.         i++;
  654.     }
  655.     i--;
  656.     end_val = e;
  657.     
  658.     while (val <= end_val) {
  659.         i++;
  660.         geoplot (val, s, 3);
  661.         geoplot (val, n, 2);
  662.         if (i < nticks) 
  663.             val = start_log * tvals[i];
  664.         else {
  665.             start_log *= 10;
  666.             val = start_log;
  667.             i = 0;
  668.         }
  669.     }
  670. }
  671.  
  672. int logy_grid (w, e, s, n, dval)
  673. double w, e, s, n, dval; {
  674.     int i, nticks, test;
  675.     double val, v0, end_val, start_log, tvals[9];
  676.     
  677.     test = rint (fabs (dval)) - 1;
  678.     if (test < 0 || test > 2) test = 0;
  679.     if (test == 0) {
  680.         tvals[0] = 1.0;
  681.         nticks = 1;
  682.     }
  683.     if (test == 1) {
  684.         tvals[0] = 1.0;
  685.         tvals[1] = 2.0;
  686.         tvals[2] = 5.0;
  687.         nticks = 3;
  688.     }
  689.     else if (test == 2) {
  690.         nticks = 9;
  691.         for (i = 0; i < nticks; i++) tvals[i] = i + 1;
  692.     }
  693.     v0 = d_log10 (s);
  694.     start_log = val = pow (10.0, floor (v0));
  695.     i = 0;
  696.     while ((s - val) > SMALL) {
  697.         if (i < nticks)
  698.             val = start_log * tvals[i];
  699.         else {
  700.             val = (start_log *= 10.0);
  701.             i = 0;
  702.         }
  703.         i++;
  704.     }
  705.     i--;
  706.     end_val = n;
  707.  
  708.     while (val <= end_val) {
  709.         i++;
  710.         geoplot (w, val, 3);
  711.         geoplot (e, val, 2);
  712.         if (i < nticks)
  713.             val = start_log * tvals[i];
  714.         else {
  715.             start_log *= 10;
  716.             val = start_log;
  717.             i = 0;
  718.         }
  719.     }
  720. }
  721.  
  722. /*    MERCATOR PROJECTION MAP FUNCTIONS    */
  723.  
  724. int merc_map_boundary (w, e, s, n)
  725. double w, e, s, n; {
  726.     double x1, x2, x3, y1, y2, y3, s1, w1, val, v2;
  727.     int shade, i, nx, ny, fat_pen, thin_pen;
  728.     
  729.     if (gmtdefs.basemap_type == 1) {    /* Draw plain boundary and return */
  730.         wesn_map_boundary (w, e, s, n);
  731.         return;
  732.     }
  733.     
  734.     fat_pen = rint (gmtdefs.frame_width * gmtdefs.dpi);
  735.     thin_pen = rint (0.1 * gmtdefs.frame_width * gmtdefs.dpi);
  736.     
  737.     ps_setline (thin_pen);
  738.     if (frame_info.side[3]) {    /* Draw western boundary */
  739.         geo_to_xy (w, s, &x1, &y1);
  740.         y1 -= gmtdefs.frame_width;
  741.         geo_to_xy (w, n, &x2, &y2);
  742.         y2 += gmtdefs.frame_width;
  743.         ps_plot (x1, y1, 3);
  744.         ps_plot (x1, y2, -2);
  745.         x1 -= gmtdefs.frame_width;
  746.         ps_plot (x1, y1, 3);
  747.         ps_plot (x1, y2, -2);
  748.     }
  749.     if (frame_info.side[1]) {    /* Draw eastern boundary */
  750.         geo_to_xy (e, s, &x2, &y1);
  751.         y1 -= gmtdefs.frame_width;
  752.         geo_to_xy (e, n, &x1, &y2);
  753.         y2 += gmtdefs.frame_width;
  754.         ps_plot (x2, y1, 3);
  755.         ps_plot (x2, y2, -2);
  756.         x2 += gmtdefs.frame_width;
  757.         ps_plot (x2, y1, 3);
  758.         ps_plot (x2, y2, -2);
  759.     }
  760.     if (frame_info.side[0]) {    /* Draw southern boundary */
  761.         geo_to_xy (w, s, &x1, &y1);
  762.         x1 -= gmtdefs.frame_width;
  763.         geo_to_xy (e, s, &x2, &y2);
  764.         x2 += gmtdefs.frame_width;
  765.         ps_plot (x1, y1, 3);
  766.         ps_plot (x2, y1, -2);
  767.         y1 -= gmtdefs.frame_width;
  768.         ps_plot (x1, y1, 3);
  769.         ps_plot (x2, y1, -2);
  770.     }
  771.     if (frame_info.side[2]) {    /* Draw northern boundary */
  772.         geo_to_xy (w, n, &x1, &y1);
  773.         x1 -= gmtdefs.frame_width;
  774.         geo_to_xy (e, n, &x2, &y2);
  775.         x2 += gmtdefs.frame_width;
  776.         ps_plot (x1, y2, 3);
  777.         ps_plot (x2, y2, -2);
  778.         y2 += gmtdefs.frame_width;
  779.         ps_plot (x1, y2, 3);
  780.         ps_plot (x2, y2, -2);
  781.     }
  782.     
  783.     /* Draw frame grid for W/E boundaries */
  784.     
  785.     if (frame_info.frame_int[1] != 0.0) {
  786.         shade = TRUE;
  787.         ps_setline(fat_pen);
  788.         s1 = floor (s / frame_info.frame_int[1]) * frame_info.frame_int[1];
  789.         if (fabs (s1 - s) > SMALL) s1 += frame_info.frame_int[1];
  790.         ny = (s1 > n) ? -1 : (n - s1) / frame_info.frame_int[1];
  791.         for (i = 0; i <= ny; i++) {
  792.             val = s1 + i * frame_info.frame_int[1];
  793.             geo_to_xy (w, val, &x1, &y1);
  794.             geo_to_xy (e, val, &x2, &y2);
  795.             if (shade) {
  796.                 v2 = val + frame_info.frame_int[1];
  797.                 if (v2 > n) v2 = n;
  798.                 if (frame_info.side[3]) {
  799.                     geo_to_xy (w, v2, &x3, &y3);
  800.                     ps_plot (x1-0.5*gmtdefs.frame_width, y1, 3);
  801.                     ps_plot (x3-0.5*gmtdefs.frame_width, y3, -2);
  802.                 }
  803.                 if (frame_info.side[1]) {
  804.                     geo_to_xy (e, v2, &x3, &y3);
  805.                     ps_plot (x2+0.5*gmtdefs.frame_width, y2, 3);
  806.                     ps_plot (x3+0.5*gmtdefs.frame_width, y3, -2);
  807.                 }
  808.                 shade = FALSE;
  809.             }
  810.             else
  811.                 shade = TRUE;
  812.         }
  813.     }
  814.     
  815.     /* Draw Frame grid for N and S boundaries */
  816.     
  817.     if (frame_info.frame_int[0] != 0.0) {
  818.         shade = TRUE;
  819.         w1 = floor (w / frame_info.frame_int[0]) * frame_info.frame_int[0];
  820.         if (fabs (w1 - w) > SMALL) w1 += frame_info.frame_int[0];
  821.         nx = (w1 > e) ? -1 : (e - w1) / frame_info.frame_int[0];
  822.         for (i = 0; i <= nx; i++) {
  823.             val = w1 + i * frame_info.frame_int[0];
  824.             geo_to_xy (val, s, &x1, &y1);
  825.             geo_to_xy (val, n, &x2, &y2);
  826.             if (shade) {
  827.                 v2 = val + frame_info.frame_int[0];
  828.                 if (v2 > e) v2 = e;
  829.                 if (frame_info.side[0]) {
  830.                     geo_to_xy (v2, s, &x3, &y3);
  831.                     ps_plot (x1, y1-0.5*gmtdefs.frame_width, 3);
  832.                     ps_plot (x3, y3-0.5*gmtdefs.frame_width, -2);
  833.                 }
  834.                 if (frame_info.side[2]) {
  835.                     geo_to_xy (v2, n, &x3, &y3);
  836.                     ps_plot (x2, y2+0.5*gmtdefs.frame_width, 3);
  837.                     ps_plot (x3, y3+0.5*gmtdefs.frame_width, -2);
  838.                 }
  839.                 shade = FALSE;
  840.             }
  841.             else
  842.                 shade = TRUE;
  843.         }
  844.     }
  845.     ps_setline (thin_pen);
  846. }
  847.  
  848. /*    POLAR STEREOGRAPHIC PROJECTION MAP FUNCTIONS    */
  849.  
  850. int polar_map_boundary (w, e, s, n)
  851. double w, e, s, n; {
  852.     int i, nx, ny, shade, thin_pen, fat_pen;
  853.     double anglew, dx2w, dy2w, anglee, dx2e, dy2e;
  854.     double y0, x0, radiuss, radiusn, da, da0, az1, az2, psize;
  855.     double x1, x2, x3, y1, y2, y3, s1, w1, val, v2, dummy, r2, dr;
  856.     
  857.     if (!project_info.region) { /* Draw rectangular boundary and return */
  858.         rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  859.         return;
  860.     }
  861.     
  862.     if (!project_info.north_pole && s == -90.0) /* Cannot have southern boundary */
  863.         frame_info.side[0] = FALSE;
  864.     if (project_info.north_pole && n == 90.0) /* Cannot have northern boundary */
  865.         frame_info.side[2] = FALSE;
  866.     if ((fabs(w-e) == 360.0 || w == e)) {
  867.         frame_info.side[1] = FALSE;
  868.         frame_info.side[3] = FALSE;
  869.     }
  870.     
  871.     if (gmtdefs.basemap_type == 1) { /* Draw plain boundary and return */
  872.         wesn_map_boundary (w, e, s, n);
  873.         return;
  874.     }
  875.     
  876.     /* Here draw fancy map boundary */
  877.     
  878.     fat_pen = rint (gmtdefs.frame_width * gmtdefs.dpi);
  879.     thin_pen = rint (0.1 * gmtdefs.frame_width * gmtdefs.dpi);
  880.     ps_setline (thin_pen);
  881.     
  882.     psize = gmtdefs.frame_width;
  883.     
  884.     /* Angle of western boundary:  */
  885.     
  886.     geo_to_xy (w, n, &x1, &y1);
  887.     geo_to_xy (w, s, &x2, &y2);
  888.     anglew = d_atan2 (y1 - y2, x1 - x2);
  889.     dx2w = -psize * sin (anglew);
  890.     dy2w = psize * cos (anglew);
  891.     
  892.     /* Angle of eastern boundary:  */
  893.     
  894.     geo_to_xy (e, n, &x1, &y1);
  895.     geo_to_xy (e, s, &x2, &y2);
  896.     anglee = d_atan2 (y1 - y2, x1 - x2);
  897.     dx2e = -psize * cos (anglee);
  898.     dy2e = psize * sin (anglee);
  899.     
  900.     geo_to_xy (project_info.central_meridian, project_info.pole, &x0, &y0);
  901.     geo_to_xy (project_info.central_meridian, s, &dummy, &y1);
  902.     geo_to_xy (project_info.central_meridian, n, &dummy, &y2);
  903.     radiuss = fabs(y1 - y0);
  904.     radiusn = fabs(y2 - y0);
  905.     dr = 0.5 * psize;
  906.     
  907.     if (frame_info.side[3]) {    /* Draw western boundary */
  908.         geo_to_xy (w, n, &x1, &y1);
  909.         geo_to_xy (w, s, &x2, &y2);
  910.         ps_plot (x1+dy2w, y1-dx2w, 3);
  911.         ps_plot (x2-dy2w, y2+dx2w, -2);
  912.         x1 += dx2w;
  913.         y1 += dy2w;
  914.         x2 += dx2w;
  915.         y2 += dy2w;
  916.         ps_plot (x1+dy2w, y1-dx2w, 3);
  917.         ps_plot (x2-dy2w, y2+dx2w, -2);
  918.     }
  919.     if (frame_info.side[1]) {    /* Draw eastern boundary */
  920.         geo_to_xy (e, n, &x1, &y1);
  921.         geo_to_xy (e, s, &x2, &y2);
  922.         ps_plot (x1-dx2e, y1+dy2e, 3);
  923.         ps_plot (x2+dx2e, y2-dy2e, -2);
  924.         x1 += dy2e;
  925.         y1 += dx2e;
  926.         x2 += dy2e;
  927.         y2 += dx2e;
  928.         ps_plot (x1-dx2e, y1+dy2e, 3);
  929.         ps_plot (x2+dx2e, y2-dy2e, -2);
  930.     }
  931.     if (frame_info.side[0]) {    /* Draw southern boundary */
  932.         da0 = R2D * psize /radiuss;
  933.         geo_to_xy (e, s, &x1, &y1);
  934.         geo_to_xy (w, s, &x2, &y2);
  935.         az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  936.         az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  937.         if (project_info.north_pole) {
  938.             r2 = radiuss + psize;
  939.             da = R2D * psize / r2;
  940.             if (az1 <= az2) az1 += 360.0;
  941.             ps_arc (x0, y0, radiuss, az2-da0, az1+da0, 3);
  942.             ps_arc (x0, y0, r2, az2-da, az1+da, 3);
  943.         }
  944.         else {
  945.             r2 = radiuss - psize;
  946.             da = R2D * psize / r2;
  947.             if (az2 <= az1) az2 += 360.0;
  948.             ps_arc (x0, y0, radiuss, az1-da0, az2+da0, 3);
  949.             ps_arc (x0, y0, r2, az1-da, az2+da, 3);
  950.         }
  951.     }
  952.     if (frame_info.side[2]) {    /* Draw northern boundary */
  953.         da0 = R2D * psize / radiusn;
  954.         geo_to_xy (e, n, &x1, &y1);
  955.         geo_to_xy (w, n, &x2, &y2);
  956.         az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  957.         az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  958.         if (project_info.north_pole) {
  959.             r2 = radiusn - psize;
  960.             da = R2D * psize / r2;
  961.             if (az1 <= az2) az1 += 360.0;
  962.             ps_arc (x0, y0, radiusn, az2-da0, az1+da0, 3);
  963.             ps_arc (x0, y0, r2, az2-da, az1+da, 3);
  964.         }
  965.         else {
  966.             r2 = radiusn + psize;
  967.             da = R2D * psize / r2;
  968.             if (az2 <= az1) az2 += 360.0;
  969.             ps_arc (x0, y0, radiusn, az1-da0, az2+da0, 3);
  970.             ps_arc (x0, y0, r2, az1-da, az2+da, 3);
  971.         }
  972.     }
  973.     
  974.     /* Anotate S-N axes */
  975.     
  976.     ps_setline (fat_pen);
  977.     if (frame_info.frame_int[1] != 0.0) {
  978.         shade = TRUE;
  979.         s1 = floor(s/frame_info.frame_int[1])*frame_info.frame_int[1];
  980.         if (fabs(s1-s) > SMALL) s1 += frame_info.frame_int[1];
  981.         ny = (s1 > n) ? -1 : (n-s1) / frame_info.frame_int[1];
  982.         for (i = 0; i <= ny; i++) {
  983.             val = s1 + i*frame_info.frame_int[1];
  984.             geo_to_xy (w, val, &x1, &y1);
  985.             geo_to_xy (e, val, &x2, &y2);
  986.             if (shade) {
  987.                 v2 = val + frame_info.frame_int[1];
  988.                 if (v2 > n) v2 = n;
  989.                 if (frame_info.side[3]) {
  990.                     geo_to_xy (w, v2, &x3, &y3);
  991.                     ps_plot (x1+0.5*dx2w, y1+0.5*dy2w, 3);
  992.                     ps_plot (x3+0.5*dx2w, y3+0.5*dy2w, -2);
  993.                 }
  994.                 if (frame_info.side[1]) {
  995.                     geo_to_xy (e, v2, &x3, &y3);
  996.                     ps_plot (x2+0.5*dy2e, y2+0.5*dx2e, 3);
  997.                     ps_plot (x3+0.5*dy2e, y3+0.5*dx2e, -2);
  998.                 }
  999.                 shade = FALSE;
  1000.             }
  1001.             else
  1002.                 shade = TRUE;
  1003.         }
  1004.     }
  1005.  
  1006.     /* Anotate W-E axes */
  1007.     
  1008.     if (frame_info.frame_int[0] != 0.0) {
  1009.         shade = TRUE;
  1010.         w1 = floor(w/frame_info.frame_int[0])*frame_info.frame_int[0];
  1011.         if (fabs(w1-w) > SMALL) w1 += frame_info.frame_int[0];
  1012.         nx = (w1 > e) ? -1 : (e-w1) / frame_info.frame_int[0];
  1013.         for (i = 0; i <= nx; i++) {
  1014.             val = w1 + i*frame_info.frame_int[0];
  1015.             if (shade) {
  1016.                 v2 = val + frame_info.frame_int[0];
  1017.                 if (v2 > e) v2 = e;
  1018.                 if (frame_info.side[0]) {
  1019.                     geo_to_xy (v2, s, &x1, &y1);
  1020.                     geo_to_xy (val, s, &x2, &y2);
  1021.                     az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1022.                     az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1023.                     if (project_info.north_pole) {
  1024.                         if (az1 < az2) az1 += 360.0;
  1025.                         ps_arc (x0, y0, radiuss+dr, az2, az1, 3);
  1026.                     }
  1027.                     else {
  1028.                         if (az2 < az1) az2 += 360.0;
  1029.                         ps_arc (x0, y0, radiuss-dr, az1, az2, 3);
  1030.                     }
  1031.                 }
  1032.                 if (frame_info.side[2]) {
  1033.                     geo_to_xy (v2, n, &x1, &y1);
  1034.                     geo_to_xy (val, n, &x2, &y2);
  1035.                     az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1036.                     az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1037.                     if (project_info.north_pole) {
  1038.                         if (az1 < az2) az1 += 360.0;
  1039.                         ps_arc (x0, y0, radiusn-dr, az2, az1, 3);
  1040.                     }
  1041.                     else {
  1042.                         if (az2 < az1) az2 += 360.0;
  1043.                         ps_arc (x0, y0, radiusn+dr, az1, az2, 3);
  1044.                     }
  1045.                 }
  1046.                 shade = FALSE;
  1047.             }
  1048.             else
  1049.                 shade = TRUE;
  1050.         }
  1051.     }
  1052.     ps_setline (thin_pen);
  1053. }
  1054.  
  1055. /*    LAMBERT CONFORMAL CONIC PROJECTION MAP FUNCTIONS    */
  1056.  
  1057. int lambert_map_boundary (w, e, s, n)
  1058. double w, e, s, n; {
  1059.     int i, nx, ny, shade, fat_pen, thin_pen;
  1060.     double dx, dy, angle, dx2, dy2, y0, x0, radiuss, radiusn, da, da0, az1, az2, psize;
  1061.     double x1, x2, x3, y1, y2, y3, s1, w1, val, v2, rsize;
  1062.     
  1063.     if (!project_info.region) { /* Draw rectangular boundary and return */
  1064.         rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  1065.         return;
  1066.     }
  1067.     if (gmtdefs.basemap_type == 1) { /* Draw plain boundary and return */
  1068.         wesn_map_boundary (w, e, s, n);
  1069.         return;
  1070.     }
  1071.     
  1072.     fat_pen = rint (gmtdefs.frame_width * gmtdefs.dpi);
  1073.     thin_pen = rint (0.1 * gmtdefs.frame_width * gmtdefs.dpi);
  1074.     ps_setline (thin_pen);
  1075.     
  1076.     psize = (project_info.north_pole) ? gmtdefs.frame_width : -gmtdefs.frame_width;
  1077.     rsize = fabs (psize);
  1078.     geo_to_xy (w, n, &x1, &y1);
  1079.     geo_to_xy (w, s, &x2, &y2);
  1080.     dx = x1 - x2;
  1081.     dy = y1 - y2;
  1082.     angle = R2D*d_atan2 (dy, dx) - 90.0;
  1083.     if (fabs(angle-180.0) < SMALL) angle = 0.0;
  1084.     dx2 = rsize * cos (angle*D2R);
  1085.     dy2 = rsize * sin (angle*D2R);
  1086.     geo_to_xy (project_info.central_meridian, project_info.pole, &x0, &y0);
  1087.     geo_to_xy (e, s, &x1, &y1);
  1088.     geo_to_xy (w, n, &x2, &y2);
  1089.     radiuss = hypot (x1 - x0, y1 - y0);
  1090.     radiusn = hypot (x2 - x0, y2 - y0);
  1091.     
  1092.     if (frame_info.side[3]) {    /* Draw western boundary */
  1093.         geo_to_xy (w, s, &x1, &y1);
  1094.         geo_to_xy (w, n, &x2, &y2);
  1095.         ps_plot (x1+dy2, y1-dx2, 3);
  1096.         ps_plot (x2-dy2, y2+dx2, -2);
  1097.         x1 -= dx2;
  1098.         y1 -= dy2;
  1099.         x2 -= dx2;
  1100.         y2 -= dy2;
  1101.         ps_plot (x1+dy2, y1-dx2, 3);
  1102.         ps_plot (x2-dy2, y2+dx2, -2);
  1103.     }
  1104.     if (frame_info.side[1]) {    /* Draw eastern boundary */
  1105.         geo_to_xy (e, s, &x1, &y1);
  1106.         geo_to_xy (e, n, &x2, &y2);
  1107.         ps_plot (x1-dy2, y1-dx2, 3);
  1108.         ps_plot (x2+dy2, y2+dx2, -2);
  1109.         x1 += dx2;
  1110.         y1 -= dy2;
  1111.         x2 += dx2;
  1112.         y2 -= dy2;
  1113.         ps_plot (x1-dy2, y1-dx2, 3);
  1114.         ps_plot (x2+dy2, y2+dx2, -2);
  1115.     }
  1116.     if (frame_info.side[0]) {    /* Draw southern boundary */
  1117.         da0 = R2D*gmtdefs.frame_width/radiuss;
  1118.         da = R2D*gmtdefs.frame_width/(radiuss+psize);
  1119.         geo_to_xy (e, s, &x1, &y1);
  1120.         geo_to_xy (w, s, &x2, &y2);
  1121.         az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1122.         az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1123.         if (project_info.north_pole) {
  1124.             if (az1 <= az2) az1 += 360.0;
  1125.             ps_arc (x0, y0, radiuss, az2-da0, az1+da0, 3);
  1126.             ps_arc (x0, y0, radiuss + psize, az2-da, az1+da, 3);
  1127.         }
  1128.         else {
  1129.             if (az2 <= az1) az2 += 360.0;
  1130.             ps_arc (x0, y0, radiuss, az1-da0, az2+da0, 3);
  1131.             ps_arc (x0, y0, radiuss + psize, az1-da, az2+da, 3);
  1132.         }
  1133.     }
  1134.     if (frame_info.side[2]) {    /* Draw northern boundary */
  1135.         da0 = R2D*gmtdefs.frame_width/radiusn;
  1136.         da = R2D*gmtdefs.frame_width/(radiusn-psize);
  1137.         geo_to_xy (e, s, &x1, &y1);
  1138.         geo_to_xy (w, s, &x2, &y2);
  1139.         az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1140.         az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1141.         if (project_info.north_pole) {
  1142.             if (az1 <= az2) az1 += 360.0;
  1143.             ps_arc (x0, y0, radiusn, az2-da0, az1+da0, 3);
  1144.             ps_arc (x0, y0, radiusn - psize, az2-da, az1+da, 3);
  1145.         }
  1146.         else {
  1147.             if (az2 <= az1) az2 += 360.0;
  1148.             ps_arc (x0, y0, radiusn, az1-da0, az2+da0, 3);
  1149.             ps_arc (x0, y0, radiusn - psize, az1-da, az2+da, 3);
  1150.         }
  1151.     }
  1152.     
  1153.     /* Anotate S-N axes */
  1154.     
  1155.     ps_setline (fat_pen);
  1156.     if (frame_info.frame_int[1] != 0.0) {
  1157.         shade = TRUE;
  1158.         s1 = floor(s/frame_info.frame_int[1])*frame_info.frame_int[1];
  1159.         if (fabs(s1-s) > SMALL) s1 += frame_info.frame_int[1];
  1160.         ny = (s1 > n) ? -1 : (n-s1) / frame_info.frame_int[1];
  1161.         for (i = 0; i <= ny; i++) {
  1162.             val = s1 + i*frame_info.frame_int[1];
  1163.             geo_to_xy (w, val, &x1, &y1);
  1164.             geo_to_xy (e, val, &x2, &y2);
  1165.             if (shade) {
  1166.                 v2 = val + frame_info.frame_int[1];
  1167.                 if (v2 > n) v2 = n;
  1168.                 if (frame_info.side[3]) {
  1169.                     geo_to_xy (w, v2, &x3, &y3);
  1170.                     ps_plot (x1-0.5*dx2, y1-0.5*dy2, 3);
  1171.                     ps_plot (x3-0.5*dx2, y3-0.5*dy2, -2);
  1172.                 }
  1173.                 if (frame_info.side[1]) {
  1174.                     geo_to_xy (e, v2, &x3, &y3);
  1175.                     ps_plot (x2+0.5*dx2, y2-0.5*dy2, 3);
  1176.                     ps_plot (x3+0.5*dx2, y3-0.5*dy2, -2);
  1177.                 }
  1178.                 shade = FALSE;
  1179.             }
  1180.             else
  1181.                 shade = TRUE;
  1182.         }
  1183.     }
  1184.  
  1185.     /* Anotate W-E axes */
  1186.     
  1187.     if (frame_info.frame_int[0] != 0.0) {
  1188.         shade = TRUE;
  1189.         w1 = floor(w/frame_info.frame_int[0])*frame_info.frame_int[0];
  1190.         if (fabs(w1-w) > SMALL) w1 += frame_info.frame_int[0];
  1191.         nx = (w1 > e) ? -1 : (e-w1) / frame_info.frame_int[0];
  1192.         da = dx;
  1193.         dx = dy;
  1194.         dy = da;
  1195.         for (i = 0; i <= nx; i++) {
  1196.             val = w1 + i*frame_info.frame_int[0];
  1197.             if (shade) {
  1198.                 v2 = val + frame_info.frame_int[0];
  1199.                 if (v2 > e) v2 = e;
  1200.                 if (frame_info.side[0]) {
  1201.                     geo_to_xy (v2, s, &x1, &y1);
  1202.                     geo_to_xy (val, s, &x2, &y2);
  1203.                     az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1204.                     az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1205.                     if (project_info.north_pole) {
  1206.                         if (az1 < az2) az1 += 360.0;
  1207.                         ps_arc (x0, y0, radiuss+0.5*psize, az2, az1, 3);
  1208.                     }
  1209.                     else {
  1210.                         if (az2 < az1) az2 += 360.0;
  1211.                         ps_arc (x0, y0, radiuss+0.5*psize, az1, az2, 3);
  1212.                     }
  1213.                 }
  1214.                 if (frame_info.side[2]) {
  1215.                     geo_to_xy (v2, n, &x1, &y1);
  1216.                     geo_to_xy (val, n, &x2, &y2);
  1217.                     az1 = d_atan2 (y1 - y0, x1 - x0) * R2D;
  1218.                     az2 = d_atan2 (y2 - y0, x2 - x0) * R2D;
  1219.                     if (project_info.north_pole) {
  1220.                         if (az1 < az2) az1 += 360.0;
  1221.                         ps_arc (x0, y0, radiusn-0.5*psize, az2, az1, 3);
  1222.                     }
  1223.                     else {
  1224.                         if (az2 < az1) az2 += 360.0;
  1225.                         ps_arc (x0, y0, radiusn-0.5*psize, az1, az2, 3);
  1226.                     }
  1227.                 }
  1228.                 shade = FALSE;
  1229.             }
  1230.             else
  1231.                 shade = TRUE;
  1232.         }
  1233.     }
  1234.     ps_setline (thin_pen);
  1235. }
  1236.  
  1237. /*    OBLIQUE MERCATOR PROJECTION MAP FUNCTIONS    */
  1238.  
  1239. int oblmrc_map_boundary (w, e, s, n)
  1240. double w, e, s, n; {
  1241.     
  1242.     rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  1243. }
  1244.     
  1245. /*    MOLLWEIDE and HAMMER-AITOFF EQUAL AREA PROJECTION MAP FUNCTIONS    */
  1246.  
  1247. int ellipse_map_boundary (w, e, s, n)
  1248. double w, e, s, n; {
  1249.     
  1250.     if (!project_info.region) { /* Draw rectangular boundary and return */
  1251.         rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  1252.         return;
  1253.     }
  1254.     if (project_info.s == -90.0) /* Cannot have southern boundary */
  1255.         frame_info.side[0] = FALSE;
  1256.     if (project_info.n == 90.0) /* Cannot have northern boundary */
  1257.         frame_info.side[2] = FALSE;
  1258.  
  1259.     wesn_map_boundary (w, e, s, n);
  1260.     
  1261. }
  1262.     
  1263. int basic_map_boundary (w, e, s, n)
  1264. double w, e, s, n; {
  1265.     
  1266.     if (!project_info.region) { /* Draw rectangular boundary and return */
  1267.         rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  1268.         return;
  1269.     }
  1270.     wesn_map_boundary (w, e, s, n);
  1271. }
  1272.     
  1273. /*
  1274.  *    GENERIC MAP PLOTTING FUNCTIONS
  1275.  */
  1276.  
  1277. int wesn_map_boundary (w, e, s, n)
  1278. double w, e, s, n; {
  1279.     int i, np = 0;
  1280.     double *xx, *yy;
  1281.     
  1282.     ps_setline (gmtdefs.frame_pen);
  1283.     
  1284.     if (frame_info.side[3]) {    /* West */
  1285.         np = map_path (w, s, w, n, &xx, &yy);
  1286.         for (i = 0; i < np; i++) geoz_to_xy (xx[i], yy[i], project_info.z_level, &xx[i], &yy[i]);
  1287.         ps_line (xx, yy, np, 3, FALSE, TRUE);
  1288.         free ((char *)xx);    free ((char *)yy);
  1289.     }
  1290.     if (frame_info.side[1]) {    /* East */
  1291.         np = map_path (e, s, e, n, &xx, &yy);
  1292.         for (i = 0; i < np; i++) geoz_to_xy (xx[i], yy[i], project_info.z_level, &xx[i], &yy[i]);
  1293.         ps_line (xx, yy, np, 3, FALSE, TRUE);
  1294.         free ((char *)xx);    free ((char *)yy);
  1295.     }
  1296.     if (frame_info.side[0]) {    /* South */
  1297.         np = map_path (w, s, e, s, &xx, &yy);
  1298.         for (i = 0; i < np; i++) geoz_to_xy (xx[i], yy[i], project_info.z_level, &xx[i], &yy[i]);
  1299.         ps_line (xx, yy, np, 3, FALSE, TRUE);
  1300.         free ((char *)xx);    free ((char *)yy);
  1301.     }
  1302.     if (frame_info.side[2]) {    /* North */
  1303.         np = map_path (w, n, e, n, &xx, &yy);
  1304.         for (i = 0; i < np; i++) geoz_to_xy (xx[i], yy[i], project_info.z_level, &xx[i], &yy[i]);
  1305.         ps_line (xx, yy, np, 3, FALSE, TRUE);
  1306.         free ((char *)xx);    free ((char *)yy);
  1307.     }
  1308. }
  1309.     
  1310. int rect_map_boundary (x0, y0, x1, y1)
  1311. double x0, y0, x1, y1; {
  1312.     double xt[4], yt[4];
  1313.     
  1314.     xy_do_z_to_xy (x0, y0, project_info.z_level, &xt[0], &yt[0]);
  1315.     xy_do_z_to_xy (x1, y0, project_info.z_level, &xt[1], &yt[1]);
  1316.     xy_do_z_to_xy (x1, y1, project_info.z_level, &xt[2], &yt[2]);
  1317.     xy_do_z_to_xy (x0, y1, project_info.z_level, &xt[3], &yt[3]);
  1318.  
  1319.     ps_setline (gmtdefs.frame_pen);
  1320.     
  1321.     if (frame_info.side[3]) {    /* West */
  1322.         ps_plot (xt[0], yt[0], 3);
  1323.         ps_plot (xt[3], yt[3], -2);
  1324.     }
  1325.     if (frame_info.side[1]) {    /* East */
  1326.         ps_plot (xt[1], yt[1], 3);
  1327.         ps_plot (xt[2], yt[2], -2);
  1328.     }
  1329.     if (frame_info.side[0]) {    /* South */
  1330.         ps_plot (xt[0], yt[0], 3);
  1331.         ps_plot (xt[1], yt[1], -2);
  1332.     }
  1333.     if (frame_info.side[2]) {    /* North */
  1334.         ps_plot (xt[3], yt[3], 3);
  1335.         ps_plot (xt[2], yt[2], -2);
  1336.     }
  1337. }
  1338.     
  1339. int circle_map_boundary (w, e, s, n)
  1340. double w, e, s, n; {
  1341.     int i, nr;
  1342.     double x0, y0, a, da;
  1343.     
  1344.     if (!project_info.region) { /* Draw rectangular boundary and return */
  1345.         rect_map_boundary (0.0, 0.0, project_info.xmax, project_info.ymax);
  1346.         return;
  1347.     }
  1348.     
  1349.     ps_setline (gmtdefs.frame_pen);
  1350.     
  1351.     nr = gmtdefs.n_lon_nodes + gmtdefs.n_lat_nodes;
  1352.     if (nr >= gmt_n_alloc) get_plot_array ();
  1353.     da = 2.0 * M_PI / (nr - 1);
  1354.     for (i = 0; i < nr; i++) {
  1355.         a = i * da;
  1356.         x0 = project_info.r * cos (a);
  1357.         y0 = project_info.r * sin (a);
  1358.         xy_do_z_to_xy (x0, y0, project_info.z_level, &gmt_x_plot[i], &gmt_y_plot[i]);
  1359.     }
  1360.     geoz_to_xy (project_info.central_meridian, project_info.pole, project_info.z_level, &x0, &y0);
  1361.     ps_transrotate (x0, y0, 0.0);
  1362.     ps_line (gmt_x_plot, gmt_y_plot, nr, 3, FALSE, TRUE);
  1363.     ps_rotatetrans (-x0, -y0, 0.0);
  1364. }
  1365.     
  1366. int theta_r_map_boundary (w, e, s, n)
  1367. double w, e, s, n; {
  1368.     int i, nr;
  1369.     double x0, y0, x, y, a, da;
  1370.     
  1371.     ps_setline (gmtdefs.frame_pen);
  1372.     
  1373.     if (s == 0.0) frame_info.side[0] = 0;    /* No donuts, please */
  1374.     nr = gmtdefs.n_lon_nodes;
  1375.     if (nr >= gmt_n_alloc) get_plot_array ();
  1376.     da = fabs (project_info.e - project_info.w) / (nr - 1);
  1377.     for (i = 0; i < nr; i++) {
  1378.         a = project_info.w + i * da;
  1379.         x = project_info.r * cosd (a);
  1380.         y = project_info.r * sind (a);
  1381.         xy_do_z_to_xy (x, y, project_info.z_level, &gmt_x_plot[i], &gmt_y_plot[i]);
  1382.     }
  1383.     geoz_to_xy (project_info.central_meridian, project_info.pole, project_info.z_level, &x0, &y0);
  1384.     ps_transrotate (x0, y0, 0.0);
  1385.     ps_line (gmt_x_plot, gmt_y_plot, nr, 3, FALSE, TRUE);
  1386.     if (frame_info.side[0]) {
  1387.         double r;
  1388.         geo_to_xy (0.0, project_info.s, &r, &y);
  1389.         r -= x0;
  1390.         for (i = 0; i < nr; i++) {
  1391.             a = project_info.w + i * da;
  1392.             x = r * cosd (a);
  1393.             y = r * sind (a);
  1394.             xy_do_z_to_xy (x, y, project_info.z_level, &gmt_x_plot[i], &gmt_y_plot[i]);
  1395.         }
  1396.         ps_line (gmt_x_plot, gmt_y_plot, nr, 3, FALSE, TRUE);
  1397.     }
  1398.     ps_rotatetrans (-x0, -y0, 0.0);
  1399.     
  1400.     if (project_info.w != project_info.e) {    /* Must draw radial borders */
  1401.         double xx[2], yy[2];
  1402.         if (frame_info.side[1]) {
  1403.             geoz_to_xy (project_info.e, project_info.s, project_info.z_level, &xx[0], &yy[0]);
  1404.             geoz_to_xy (project_info.e, project_info.n, project_info.z_level, &xx[1], &yy[1]);
  1405.             ps_line (xx, yy, 2, 3, FALSE, TRUE);
  1406.         }
  1407.         if (frame_info.side[3]) {
  1408.             geoz_to_xy (project_info.w, project_info.s, project_info.z_level, &xx[0], &yy[0]);
  1409.             geoz_to_xy (project_info.w, project_info.n, project_info.z_level, &xx[1], &yy[1]);
  1410.             ps_line (xx, yy, 2, 3, FALSE, TRUE);
  1411.         }
  1412.     }
  1413. }
  1414.     
  1415. int map_latline (lat, west, east)        /* Draws a line of constant latitude */
  1416. double lat, west, east; {
  1417.     int nn;
  1418.     double *llon, *llat;
  1419.     char text[20];
  1420.     
  1421.     nn = latpath (lat, west, east, &llon, &llat);
  1422.     
  1423.     gmt_n_plot = geo_to_xy_line (llon, llat, nn);
  1424.     sprintf (text, "Lat = %lg\0", lat);
  1425.     ps_comment (text);
  1426.     gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  1427.     
  1428.     free ((char *)llon);
  1429.     free ((char *)llat);
  1430. }
  1431.     
  1432. int map_lonline (lon, south, north)    /* Draws a line of constant longitude */
  1433. double lon, south, north; {
  1434.     int nn;
  1435.     double *llon, *llat;
  1436.     char text[20];
  1437.     
  1438.     nn = lonpath (lon, south, north, &llon, &llat);
  1439.  
  1440.     gmt_n_plot = geo_to_xy_line (llon, llat, nn);
  1441.     sprintf (text, "Lon = %lg\0", lon);
  1442.     ps_comment (text);
  1443.     gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  1444.     
  1445.     free ((char *)llon);
  1446.     free ((char *)llat);
  1447. }
  1448.  
  1449. int map_lontick (lon, south, north)
  1450. double lon, south, north; {
  1451.     int i, nc;
  1452.     struct XINGS *xings;
  1453.     
  1454.     for (i = 0; i < (nc = map_loncross (lon, south, north, &xings)); i++) {
  1455.         map_tick (xings[i].xx, xings[i].yy, xings[i].sides, xings[i].angle, xings[i].nx);
  1456.         free ((char *)xings);
  1457.     }
  1458. }
  1459.  
  1460. int map_lattick (lat, west, east)
  1461. double lat, west, east; {
  1462.     int i, nc;
  1463.     
  1464.     struct XINGS *xings;
  1465.     
  1466.     for (i = 0; i < (nc = map_latcross (lat, west, east, &xings)); i++) {
  1467.         map_tick (xings[i].xx, xings[i].yy, xings[i].sides, xings[i].angle, xings[i].nx);
  1468.         free ((char *)xings);
  1469.     }
  1470. }
  1471.  
  1472. int map_tick (xx, yy, sides, angles, nx)
  1473. double *xx, *yy, *angles;
  1474. int *sides, nx; {
  1475.     double angle, xl, yl, xt, yt;
  1476.     int i;
  1477.     
  1478.     for (i = 0; i < nx; i++) {
  1479.         if (!project_info.edge[sides[i]]) continue;
  1480.         if (!frame_info.side[sides[i]]) continue;
  1481.         angle = angles[i];
  1482.         xl = 0.5 * gmtdefs.tick_length * cos (angle * D2R);
  1483.         yl = 0.5 * gmtdefs.tick_length * sin (angle * D2R);
  1484.         xy_do_z_to_xy (xx[i], yy[i], project_info.z_level, &xt, &yt);
  1485.         ps_plot (xt, yt, 3);
  1486.         xy_do_z_to_xy (xx[i]+xl, yy[i]+yl, project_info.z_level, &xt, &yt);
  1487.         ps_plot (xt, yt, -2);
  1488.     }
  1489. }
  1490.  
  1491. int map_symbol_ew (lat, label, west, east, anot)
  1492. double lat, west, east;
  1493. char *label;
  1494. BOOLEAN anot; {
  1495.     int i, nc;
  1496.     struct XINGS *xings;
  1497.     
  1498.     for (i = 0; i < (nc = map_latcross (lat, west, east, &xings)); i++) {
  1499.         map_symbol (xings[i].xx, xings[i].yy, xings[i].sides, xings[i].angle, label, xings[i].nx, 1, anot);
  1500.         free ((char *)xings);
  1501.     }
  1502. }
  1503.  
  1504. int map_symbol_ns (lon, label, south, north, anot)
  1505. double lon, south, north;
  1506. char *label;
  1507. BOOLEAN anot; {
  1508.     int i, nc;
  1509.     struct XINGS *xings;
  1510.     
  1511.     for (i = 0; i < (nc = map_loncross (lon, south, north, &xings)); i++) {
  1512.         map_symbol (xings[i].xx, xings[i].yy, xings[i].sides, xings[i].angle, label, xings[i].nx, 0, anot);
  1513.         free ((char *)xings);
  1514.     }
  1515. }
  1516.  
  1517. int map_symbol (xx, yy, sides, line_angles, label, nx, type, anot)
  1518. double *xx, *yy, *line_angles;
  1519. int *sides, nx, type;    /* type = 0 for lon and 1 for lat */
  1520. char *label;
  1521. BOOLEAN anot; {
  1522.     double line_angle, text_angle, len, dx, dy, angle, ca, sa, xt1, yt1, zz;
  1523.     int i, justify;
  1524.     
  1525.     len = ((gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0) + gmtdefs.anot_offset;
  1526.     
  1527.     for (i = 0; i < nx; i++) {
  1528.     
  1529.         if (prepare_label (line_angles[i], sides[i], xx[i], yy[i], type, &line_angle, &text_angle, &justify)) continue;
  1530.         
  1531.         angle = line_angle * D2R;
  1532.         ca = cos (angle);    sa = sin (angle);
  1533.         dx = gmtdefs.tick_length * ca;
  1534.         dy = gmtdefs.tick_length * sa;
  1535.         z_to_zz (project_info.z_level, &zz);
  1536.         xyz_to_xy (xx[i], yy[i], zz, &xt1, &yt1);
  1537.         ps_plot (xt1, yt1, 3);
  1538.         xyz_to_xy (xx[i]+dx, yy[i]+dy, zz, &xt1, &yt1);
  1539.         ps_plot (xt1, yt1, -2);
  1540.         xx[i] += len * ca;
  1541.         yy[i] += len * sa;
  1542.         xyz_to_xy (xx[i], yy[i], zz, &xt1, &yt1);
  1543.             
  1544.         if (project_info.three_D) {
  1545.             int upside = FALSE, k;
  1546.             double xp[2], yp[2], xt2, xt3, yt2, yt3, del_y;
  1547.             double size, xsize, ysize, xshrink, yshrink, tilt, baseline_shift, cb, sb, a;
  1548.                 
  1549.             upside = (z_project.quadrant == 1 || z_project.quadrant == 4);
  1550.             cb = cosd (text_angle);    sb = sind (text_angle);
  1551.             if (sides[i]%2 == 0 && (justify%2 == 0)) {
  1552.                 if (upside) {
  1553.                     k = (sides[i] == 0) ? 2 : 10;
  1554.                     text_angle += 180.0;
  1555.                 }
  1556.                 else
  1557.                     k = justify;
  1558.                 del_y = 0.5 * gmtdefs.anot_font_size * 0.732 * (k/4) / gmt_ppu[gmtdefs.measure_unit];
  1559.                 justify = 2;
  1560.                 xx[i] += del_y * ca;    yy[i] += del_y * sa;
  1561.                 xyz_to_xy (xx[i], yy[i], zz, &xt1, &yt1);
  1562.             }
  1563.             else {
  1564.                 del_y = -0.5 * gmtdefs.anot_font_size * 0.732 * (justify/4) / gmt_ppu[gmtdefs.measure_unit];
  1565.                 if (upside) {
  1566.                     if (sides[i]%2) del_y = -del_y;
  1567.                     text_angle += 180.0;
  1568.                     justify = (justify == 5) ? 7 : 5;
  1569.                 }
  1570.                 justify -= 4;
  1571.                 switch (sides[i]) {
  1572.                     case 0:
  1573.                         a = (justify == 1) ? line_angle + 90.0 : line_angle - 90.0;
  1574.                         break;
  1575.                     case 1:
  1576.                         a = line_angle + 90.0;
  1577.                         break;
  1578.                     case 2:
  1579.                         a = (justify == 1) ? line_angle + 90.0 : line_angle - 90.0;
  1580.                         break;
  1581.                     case 3:
  1582.                         a = line_angle - 90.0;
  1583.                         break;
  1584.                 }
  1585.                 ca = cosd (a);    sa = sind (a);
  1586.                 xx[i] += del_y * ca;    yy[i] += del_y * sa;
  1587.                 xyz_to_xy (xx[i], yy[i], zz, &xt1, &yt1);
  1588.             }
  1589.             xp[0] = xx[i] + len * cb;    yp[0] = yy[i] + len * sb;
  1590.             xp[1] = xx[i] - len * sb;    yp[1] = yy[i] + len * cb;
  1591.             xyz_to_xy (xp[0], yp[0], zz, &xt2, &yt2);
  1592.             xyz_to_xy (xp[1], yp[1], zz, &xt3, &yt3);
  1593.             xshrink = hypot (xt2-xt1, yt2-yt1) / hypot (xp[0]-xx[i], yp[0]-yy[i]);
  1594.             yshrink = hypot (xt3-xt1, yt3-yt1) / hypot (xp[1]-xx[i], yp[1]-yy[i]);
  1595.             baseline_shift = d_atan2 (yt2 - yt1, xt2 - xt1) - d_atan2 (yp[0] - yy[i], xp[0] - xx[i]);
  1596.             tilt = 90.0 - R2D * (d_atan2 (yt3 - yt1, xt3 - xt1) - d_atan2 (yt2 - yt1, xt2 - xt1));
  1597.             tilt = tand (tilt);
  1598.             size = gmtdefs.anot_font_size * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit];
  1599.             xsize = size * xshrink;
  1600.             ysize = size * yshrink;
  1601.             /* Temporarely modify meaning of F0 */
  1602.             printf ("/F0 {pop /%s findfont [%lg 0 %lg %lg 0 0] makefont setfont} bind def\n",
  1603.             font_name[gmtdefs.anot_font], xsize, ysize * tilt, ysize);
  1604.             ps_setfont (0);
  1605.             text_angle += (R2D * baseline_shift);
  1606.         }
  1607.         if (anot) ps_text (xt1, yt1, gmtdefs.anot_font_size, label, text_angle, justify, 0);
  1608.     }
  1609. }
  1610.  
  1611. int map_latcross (lat, west, east, xings)
  1612. double lat, east, west;
  1613. struct XINGS *xings[]; {
  1614.     int i, go = FALSE, nx, nc = 0, n_alloc = 50;
  1615.     double lon, lon_old, this_x, this_y, last_x, last_y, xlon[2], xlat[2];
  1616.     double get_angle ();
  1617.     struct XINGS *X;
  1618.     
  1619.     X = (struct XINGS *) memory (CNULL, n_alloc, sizeof (struct XINGS), "map_loncross");
  1620.         
  1621.     lon_old = west - SMALL;
  1622.     map_outside (lon_old, lat);
  1623.     geo_to_xy (lon_old, lat, &last_x, &last_y);
  1624.     for (i = 1; i <= gmtdefs.n_lon_nodes; i++) {
  1625.         lon = (i == gmtdefs.n_lon_nodes) ? east + SMALL : west + i * gmtdefs.dlon;
  1626.         map_outside (lon, lat);
  1627.         geo_to_xy (lon, lat, &this_x, &this_y);
  1628.         nx = 0;
  1629.         if ( break_through (lon_old, lat, lon, lat) ) {    /* Crossed map boundary */
  1630.             nx = map_crossing (lon_old, lat, lon, lat, xlon, xlat, X[nc].xx, X[nc].yy, X[nc].sides);
  1631.             if (nx == 1) X[nc].angle[0] = get_angle (lon_old, lat, lon, lat);
  1632.             if (nx == 2) X[nc].angle[1] = X[nc].angle[0] + 180.0;
  1633.             if (gmt_corner > 0) {
  1634.                 X[nc].sides[0] = (gmt_corner%4 > 1) ? 1 : 3;
  1635.                 gmt_corner = 0;
  1636.             }
  1637.         }
  1638.         if (gmt_world_map) wrap_around_check (X[nc].angle, last_x, last_y, this_x, this_y, X[nc].xx, X[nc].yy, X[nc].sides, &nx);
  1639.         if (nx == 2 && (fabs (X[nc].xx[1] - X[nc].xx[0]) - gmt_map_width) < SMALL && !gmt_world_map)
  1640.             go = FALSE;
  1641.         else if (nx > 0)
  1642.             go = TRUE;
  1643.         if (go) {
  1644.             X[nc].nx = nx;
  1645.             nc++;
  1646.             if (nc == n_alloc) {
  1647.                 n_alloc += 50;
  1648.                 X = (struct XINGS *) memory ((char *)X, n_alloc, sizeof (struct XINGS), "map_loncross");
  1649.             }
  1650.             go = FALSE;
  1651.         }
  1652.         lon_old = lon;
  1653.         last_x = this_x;    last_y = this_y;
  1654.     }
  1655.     
  1656.     if (nc > 0) {
  1657.         X = (struct XINGS *) memory ((char *)X, nc, sizeof (struct XINGS), "map_loncross");
  1658.         *xings = X;
  1659.     }
  1660.     else
  1661.         free ((char *)X);
  1662.     
  1663.     return (nc);
  1664. }
  1665.  
  1666. int map_loncross (lon, south, north, xings)
  1667. double lon, south, north;
  1668. struct XINGS *xings[]; {
  1669.     int go = FALSE, j, nx, nc = 0, n_alloc = 50;
  1670.     double lat, lat_old, this_x, this_y, last_x, last_y, xlon[2], xlat[2];
  1671.     double get_angle ();
  1672.     struct XINGS *X;
  1673.     
  1674.     X = (struct XINGS *) memory (CNULL, n_alloc, sizeof (struct XINGS), "map_loncross");
  1675.     
  1676.     lat_old = ((south - SMALL) >= -90.0) ? south - SMALL : south;    /* Outside */
  1677.     if ((north + SMALL) <= 90.0) north += SMALL;
  1678.     map_outside (lon, lat_old);
  1679.     geo_to_xy (lon, lat_old, &last_x, &last_y);
  1680.     for (j = 1; j <= gmtdefs.n_lat_nodes; j++) {
  1681.         lat = (j == gmtdefs.n_lat_nodes) ? north: south + j * gmtdefs.dlat;
  1682.         map_outside (lon, lat);
  1683.         geo_to_xy (lon, lat, &this_x, &this_y);
  1684.         nx = 0;
  1685.         if ( break_through (lon, lat_old, lon, lat) ) {    /* Crossed map boundary */
  1686.             nx = map_crossing (lon, lat_old, lon, lat, xlon, xlat, X[nc].xx, X[nc].yy, X[nc].sides);
  1687.             if (nx == 1) X[nc].angle[0] = get_angle (lon, lat_old, lon, lat);
  1688.             if (nx == 2) X[nc].angle[1] = X[nc].angle[0] + 180.0;
  1689.             if (gmt_corner > 0) {
  1690.                 X[nc].sides[0] = (gmt_corner < 3) ? 0 : 2;
  1691.                 gmt_corner = 0;
  1692.             }
  1693.         }
  1694.         if (gmt_world_map) wrap_around_check (X[nc].angle, last_x, last_y, this_x, this_y, X[nc].xx, X[nc].yy, X[nc].sides, &nx);
  1695.         if (nx == 2 && (fabs (X[nc].xx[1] - X[nc].xx[0]) - gmt_map_width) < SMALL && !gmt_world_map)
  1696.             go = FALSE;
  1697.         else if (nx > 0)
  1698.             go = TRUE;
  1699.         if (go) {
  1700.             X[nc].nx = nx;
  1701.             nc++;
  1702.             if (nc == n_alloc) {
  1703.                 n_alloc += 50;
  1704.                 X = (struct XINGS *) memory ((char *)X, n_alloc, sizeof (struct XINGS), "map_loncross");
  1705.             }
  1706.             go = FALSE;
  1707.         }
  1708.         lat_old = lat;
  1709.         last_x = this_x;    last_y = this_y;
  1710.     }
  1711.     
  1712.     if (nc > 0) {
  1713.         X = (struct XINGS *) memory ((char *)X, nc, sizeof (struct XINGS), "map_loncross");
  1714.         *xings = X;
  1715.     }
  1716.     else
  1717.         free ((char *)X);
  1718.     
  1719.     return (nc);
  1720. }
  1721.  
  1722. int map_gridlines (w, e, s, n)
  1723. double w, e, s, n; {
  1724.     int i, nx, ny;
  1725.     BOOLEAN dash;
  1726.     double dx, dy, w1, s1;
  1727.     
  1728.     if (gmtdefs.grid_cross_size > 0.0) return;
  1729.     
  1730.     dx = fabs (frame_info.grid_int[0]);
  1731.     dy = fabs (frame_info.grid_int[1]);
  1732.     
  1733.     dash = (!MAPPING && (dx > 0.0 || dy > 0.0));
  1734.     
  1735.     ps_setline (gmtdefs.grid_pen);
  1736.     if (dash) ps_setdash ("1 4", 0);
  1737.     if (dx > 0.0 && project_info.xyz_projection[0] == LOG10)
  1738.         logx_grid (w, e, s, n, dx);
  1739.     else if (dx > 0.0) {    /* Draw grid lines that go E to W */
  1740.         w1 = floor (w / dx) * dx;
  1741.         if ((w - w1) > SMALL) w1 += dx;
  1742.         nx = (w1 > e) ? -1 : (e - w1) / dx;
  1743.         for (i = 0; i <= nx; i++) map_lonline (w1 + i * dx, s, n);
  1744.     }
  1745.     
  1746.     if (dy > 0.0 && project_info.xyz_projection[1] == LOG10)
  1747.         logy_grid (w, e, s, n, dy);
  1748.     else if (dy > 0.0) {    /* Draw grid lines that go S to N */
  1749.         s1 = floor (s / dy) * dy;
  1750.         if ((s - s1) > SMALL) s1 += dy;
  1751.         ny = (s1 > n) ? -1 : (n - s1) / dy;
  1752.         for (i = 0; i <= ny; i++) map_latline (s1 + i * dy, w, e);
  1753.     }
  1754.     if (dash) ps_setdash ((char *)NULL, 0);
  1755. }
  1756.  
  1757. int map_gridcross (w, e, s, n)
  1758. double w, e, s, n; {
  1759.     int i, j, nx, ny;
  1760.     double dx, dy, w1, s1, lon, lat, x0, y0, x1, y1, xa, xb, ya, yb;
  1761.     double x_angle, y_angle, e1, n1, xt1, xt2, yt1, yt2, C, S, L;
  1762.     
  1763.     if (!MAPPING) return;
  1764.     if (gmtdefs.grid_cross_size <= 0.0) return;
  1765.     
  1766.     dx = fabs (frame_info.grid_int[0]);
  1767.     dy = fabs (frame_info.grid_int[1]);
  1768.     
  1769.     if (dx <= 0.0 || dy <= 0.0) return;
  1770.     
  1771.     ps_setline (gmtdefs.grid_pen);
  1772.     
  1773.     w1 = floor (w / dx) * dx;
  1774.     e1 = ceil (e / dx) * dx;
  1775.     nx = (e1 - w1) / dx;
  1776.     s1 = floor (s / dy) * dy;
  1777.     n1 = ceil (n / dy) * dy;
  1778.     ny = (n1 - s1) / dy;
  1779.     
  1780.     L = 0.5 * gmtdefs.grid_cross_size;
  1781.     
  1782.     for (i = 0; i <= nx; i++) {
  1783.         lon = w1 + i * dx;
  1784.         for (j = 0; j <= ny; j++) {
  1785.             lat = s1 + j * dy;
  1786.             
  1787.             if (!map_outside (lon, lat)) {    /* Inside map */
  1788.             
  1789.                 geo_to_xy (lon, lat, &x0, &y0);
  1790.                 geo_to_xy (lon + gmtdefs.dlon, lat, &x1, &y1);
  1791.                 x_angle = d_atan2 (y1-y0, x1-x0);
  1792.                 y_angle = x_angle + M_PI_2;
  1793.                 C = cos (x_angle);
  1794.                 S = sin (x_angle);
  1795.                 xa = x0 - L * C;
  1796.                 xb = x0 + L * C;
  1797.                 ya = y0 - L * S;
  1798.                 yb = y0 + L * S;
  1799.                 
  1800.                 /* Clip to map */
  1801.                 
  1802.                 if (xa < 0.0) xa = 0.0;
  1803.                 if (xb < 0.0) xb = 0.0;
  1804.                 if (ya < 0.0) ya = 0.0;
  1805.                 if (yb < 0.0) yb = 0.0;
  1806.                 if (xa > gmt_map_width) xa = gmt_map_width;
  1807.                 if (xb > gmt_map_width) xb = gmt_map_width;
  1808.                 if (ya > gmt_map_height) ya = gmt_map_height;
  1809.                 if (yb > gmt_map_height) yb = gmt_map_height;
  1810.                 
  1811.                 /* 3-D projection */
  1812.                 
  1813.                 xy_do_z_to_xy (xa, ya, project_info.z_level, &xt1, &yt1);
  1814.                 xy_do_z_to_xy (xb, yb, project_info.z_level, &xt2, &yt2);
  1815.                 ps_plot (xt1, yt1, 3);
  1816.                 ps_plot (xt2, yt2, -2);
  1817.                 
  1818.                 C = cos (y_angle);
  1819.                 S = sin (y_angle);
  1820.                 xa = x0 - L * C;
  1821.                 xb = x0 + L * C;
  1822.                 ya = y0 - L * S;
  1823.                 yb = y0 + L * S;
  1824.                 
  1825.                 /* Clip to map */
  1826.                 
  1827.                 if (xa < 0.0) xa = 0.0;
  1828.                 if (xb < 0.0) xb = 0.0;
  1829.                 if (ya < 0.0) ya = 0.0;
  1830.                 if (yb < 0.0) yb = 0.0;
  1831.                 if (xa > gmt_map_width) xa = gmt_map_width;
  1832.                 if (xb > gmt_map_width) xb = gmt_map_width;
  1833.                 if (ya > gmt_map_height) ya = gmt_map_height;
  1834.                 if (yb > gmt_map_height) yb = gmt_map_height;
  1835.                 
  1836.                 /* 3-D projection */
  1837.                 
  1838.                 xy_do_z_to_xy (xa, ya, project_info.z_level, &xt1, &yt1);
  1839.                 xy_do_z_to_xy (xb, yb, project_info.z_level, &xt2, &yt2);
  1840.                 ps_plot (xt1, yt1, 3);
  1841.                 ps_plot (xt2, yt2, -2);
  1842.             }
  1843.         }
  1844.     }
  1845. }
  1846.  
  1847. int map_tickmarks (w, e, s, n)
  1848. double w, e, s, n; {
  1849.     int i, nx, ny;
  1850.     double dx, dy, w1, s1;
  1851.     
  1852.     if (!MAPPING || gmtdefs.basemap_type == 0) return;
  1853.     
  1854.     ps_setline (gmtdefs.tick_pen);
  1855.     dx = fabs (frame_info.frame_int[0]);
  1856.     dy = fabs (frame_info.frame_int[1]);
  1857.     
  1858.     on_border_is_outside = TRUE;    /* Temporarily, points on the border are outside */
  1859.     
  1860.     if (dx > 0.0 && dx != fabs (frame_info.anot_int[0])) {    /* Draw grid lines that go E to W */
  1861.         w1 = floor (w / dx) * dx;
  1862.         if (fabs (w1 - w) > SMALL) w1 += dx;
  1863.         nx = (w1 > e) ? -1 : (e - w1) / dx;
  1864.         for (i = 0; i <= nx; i++) map_lontick (w1 + i * dx, s, n);
  1865.     }
  1866.     
  1867.     if (dy > 0.0 && dy != fabs (frame_info.anot_int[1])) {    /* Draw grid lines that go S to N */
  1868.         s1 = floor (s / dy) * dy;
  1869.         if (fabs (s1 - s) > SMALL) s1 += dy;
  1870.         ny = (s1 > n) ? -1 : (n - s1) / dy;
  1871.         for (i = 0; i <= ny; i++) map_lattick (s1 + i * dy, w, e);
  1872.     }
  1873.     
  1874.     on_border_is_outside = FALSE;    /* Reset back to default */
  1875. }
  1876.  
  1877. int map_anotate (w, e, s, n)
  1878. double w, e, s, n; {
  1879.     double s1, w1, val, dx, dy, x, y;
  1880.     int do_minutes, do_seconds, move_up, i, nx, ny, done_zero = FALSE, anot, gmt_world_map_save;
  1881.     char label[80];
  1882.     
  1883.     if (frame_info.header[0]) {    /* Make plot header */
  1884.         move_up = (MAPPING || frame_info.side[2] == 2);
  1885.         ps_setfont (gmtdefs.header_font);
  1886.         x = project_info.xmax * 0.5;
  1887.         y = project_info.ymax + ((gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0) + 2.5 * gmtdefs.anot_offset;
  1888.         y += ((move_up) ? (gmtdefs.anot_font_size + gmtdefs.label_font_size) / gmt_ppu[gmtdefs.measure_unit] : 0.0) + 2.5 * gmtdefs.anot_offset;
  1889.         if (project_info.three_D && project_info.z_scale == 0.0) {    /* Only do this if flat 2-D plot */
  1890.             double size, xsize, ysize;
  1891.             
  1892.             ps_setfont (0);
  1893.             xy_do_z_to_xy (x, y, project_info.z_level, &x, &y);
  1894.             size = gmtdefs.header_font_size * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit];
  1895.             xsize = size * z_project.xshrink[0];
  1896.             ysize = size * z_project.yshrink[0];
  1897.             printf ("/F0 {pop /%s findfont [%lg 0 %lg %lg 0 0] makefont setfont} bind def\n",
  1898.                 font_name[gmtdefs.header_font], xsize, ysize * z_project.tilt[0], ysize);
  1899.             
  1900.             ps_text (x, y, gmtdefs.header_font_size, frame_info.header, z_project.phi[0], -2, 0);
  1901.             printf ("/F0 {/Helvetica Y} bind def\n");    /* Reset definition of F0 */
  1902.             ps_setfont (gmtdefs.header_font);
  1903.         }
  1904.         else if (!project_info.three_D)
  1905.             ps_text (x, y, gmtdefs.header_font_size, frame_info.header, 0.0, -2, 0);
  1906.     }
  1907.     
  1908.     if (!MAPPING) return;    /* Annotation already done by linear_axis */
  1909.     
  1910.     ps_setfont (gmtdefs.anot_font);
  1911.     ps_setline (gmtdefs.tick_pen);
  1912.     dx = (project_info.edge[0] || project_info.edge[2]) ? fabs (frame_info.anot_int[0]) : 0.0;
  1913.     dy = (project_info.edge[1] || project_info.edge[3]) ? fabs (frame_info.anot_int[1]) : 0.0;
  1914.     
  1915.     on_border_is_outside = TRUE;    /* Temporarily, points on the border are outside */
  1916.     gmt_world_map_save = gmt_world_map;
  1917.     if (project_info.region) gmt_world_map = FALSE;
  1918.     
  1919.     if (dx > 0.0) {    /* Anotate the S and N boundaries */
  1920.         BOOLEAN full_lat_range, proj_A, proj_B, anot_0_and_360;
  1921.         
  1922.         /* Determine if we should annotate both 0 and 360 degrees */
  1923.         
  1924.         full_lat_range = (fabs (180.0 - fabs (project_info.n - project_info.s)) < SMALL);
  1925.         proj_A = (project_info.projection == MERCATOR || project_info.projection == OBLIQUE_MERC ||
  1926.             project_info.projection == WINKEL || project_info.projection == ECKERT ||
  1927.             project_info.projection == ROBINSON || project_info.projection == CYL_EQ ||
  1928.             project_info.projection == CYL_EQDIST || project_info.projection == LINEAR);
  1929.         proj_B = (project_info.projection == HAMMER || project_info.projection == MOLLWEIDE ||
  1930.             project_info.projection == SINUSOIDAL);
  1931. /*        anot_0_and_360 = (gmt_world_map_save && ((full_lat_range && proj_A) || (!full_lat_range && proj_B))); */
  1932.         anot_0_and_360 = (gmt_world_map_save && (proj_A || (!full_lat_range && proj_B)));
  1933.         
  1934.         do_minutes = (fabs (fmod (dx, 1.0)) > SMALL);
  1935.         do_seconds = (fabs (60.0 * fmod (fmod (dx, 1.0) * 60.0, 1.0)) >= 1.0);
  1936.         w1 = floor (w / dx) * dx;
  1937.         if (fabs (w1 - w) > SMALL) w1 += dx;
  1938.         nx = (w1 > e) ? -1 : (e - w1) / dx;
  1939.         for (i = 0; i <= nx; i++) {
  1940.             val = w1 + i * dx;
  1941.             if (val == 0.0) done_zero = TRUE;
  1942.             get_anot_label (val, label, do_minutes, do_seconds, 0);
  1943.             anot = anot_0_and_360 || !(done_zero && val == 360.0);
  1944.             map_symbol_ns (val, label, s, n, anot);
  1945.         }
  1946.     }
  1947.     
  1948.     if (dy > 0.0) {    /* Anotate W and E boundaries */
  1949.         do_minutes = (fabs (fmod (dy, 1.0)) > SMALL);
  1950.         do_seconds = (fabs (60.0 * fmod (fmod (dy, 1.0) * 60.0, 1.0)) >= 1.0);
  1951.         s1 = floor (s / dy) * dy;
  1952.         if (fabs (s1 - s) > SMALL) s1 += frame_info.anot_int[1];
  1953.         ny = (s1 > n) ? -1: (n - s1) / dy;
  1954.         for (i = 0; i <= ny; i++) {
  1955.             val = s1 + i * frame_info.anot_int[1];
  1956.             if (project_info.polar && fabs (val) == 90.0) continue;
  1957.             get_anot_label (val, label, do_minutes, do_seconds, 1);
  1958.             map_symbol_ew (val, label, w, e, TRUE);
  1959.         }
  1960.     }
  1961.     
  1962.     on_border_is_outside = FALSE;    /* Reset back to default */
  1963.     if (project_info.region) gmt_world_map = gmt_world_map_save;
  1964. }
  1965.  
  1966. int map_boundary (w, e, s, n)
  1967. double w, e, s, n; {
  1968.     switch (project_info.projection) {
  1969.         case LINEAR:
  1970.             if (MAPPING)    /* xy is lonlat */
  1971.                 merc_map_boundary (w, e, s, n);
  1972.             else if (project_info.three_D)
  1973.                 basemap_3D (3);
  1974.             else
  1975.                 linear_map_boundary (w, e, s, n);
  1976.             break;
  1977.         case POLAR:
  1978.             theta_r_map_boundary (w, e, s, n);
  1979.             break;
  1980.         case MERCATOR:
  1981.         case CYL_EQ:
  1982.         case CYL_EQDIST:
  1983.             merc_map_boundary (w, e, s, n);
  1984.             break;
  1985.         case LAMBERT:
  1986.             lambert_map_boundary (w, e, s, n);
  1987.             break;
  1988.         case OBLIQUE_MERC:
  1989.             oblmrc_map_boundary (w, e, s, n);
  1990.             break;
  1991.         case STEREO:
  1992.         case ORTHO:
  1993.         case LAMB_AZ_EQ:
  1994.         case AZ_EQDIST:
  1995.             if (project_info.polar)
  1996.                 polar_map_boundary (w, e, s, n);
  1997.             else
  1998.                 circle_map_boundary (w, e, s, n);
  1999.             break;
  2000.         case HAMMER:
  2001.         case MOLLWEIDE:
  2002.         case SINUSOIDAL:
  2003.             ellipse_map_boundary (w, e, s, n);
  2004.             break;
  2005.         case TM:
  2006.         case UTM:
  2007.         case ALBERS:
  2008.         case CASSINI:
  2009.         case WINKEL:
  2010.         case ECKERT:
  2011.         case ROBINSON:
  2012.             basic_map_boundary (w, e, s, n);
  2013.             break;
  2014.     }
  2015.     
  2016.     if (project_info.three_D) vertical_axis (3);
  2017. }
  2018.  
  2019. /* map_basemap will create a basemap for the given area. 
  2020.  * Scaling and wesn are assumed to be passed throught the project_info-structure (see gmt_project.h)
  2021.  * Tickmark info are passed through the frame_info-structure
  2022.  *
  2023.  */
  2024.  
  2025. int map_basemap () {
  2026.     double w, e, s, n;
  2027.  
  2028.     w = project_info.w;    e = project_info.e;    s = project_info.s;    n = project_info.n;
  2029.     
  2030.     ps_comment ("Start of basemap");
  2031.     
  2032.     ps_comment ("Map gridlines");
  2033.     map_gridlines (w, e, s, n);
  2034.     map_gridcross (w, e, s, n);
  2035.     
  2036.     ps_comment ("Map tickmarks");
  2037.     map_tickmarks (w, e, s, n);
  2038.     
  2039.     ps_comment ("Map anotations");
  2040.     map_anotate (w, e, s, n);
  2041.     
  2042.     ps_comment ("Map boundaries");
  2043.     map_boundary (w, e, s, n);
  2044.     
  2045.     ps_comment ("End of basemap");
  2046.     
  2047.     return (0);
  2048. }
  2049.  
  2050. int basemap_3D (mode)
  2051. int mode; {
  2052.     /* Mode means: 1 = background axis, 2 = foreground axis, 3 = all */
  2053.     BOOLEAN go[4], back;
  2054.     int i;
  2055.     
  2056.     back = (mode % 2);
  2057.     for (i = 0; i < 4; i++) go[i] = (mode == 3) ? TRUE : ((back) ? z_project.draw[i] : !z_project.draw[i]);
  2058.     
  2059.     if (go[0] && frame_info.side[0])    /* South or lower x-axis */
  2060.         xyz_axis3D (0, 'x', frame_info.anot_int[0], frame_info.frame_int[0],
  2061.         frame_info.label[0], project_info.xyz_projection[0], frame_info.anot_type[0], frame_info.side[0]-1);
  2062.     
  2063.     if (go[2] && frame_info.side[2])    /* North or upper x-axis */
  2064.         xyz_axis3D (2, 'x', frame_info.anot_int[0], frame_info.frame_int[0],
  2065.         frame_info.label[0], project_info.xyz_projection[0], frame_info.anot_type[0], frame_info.side[2]-1);
  2066.     
  2067.     if (go[3] && frame_info.side[3])    /* West or left y-axis */
  2068.         xyz_axis3D (3, 'y', frame_info.anot_int[1], frame_info.frame_int[1],
  2069.         frame_info.label[1], project_info.xyz_projection[1], frame_info.anot_type[1], frame_info.side[3]-1);
  2070.             
  2071.     if (go[1] && frame_info.side[1])    /* East or right y-axis */
  2072.         xyz_axis3D (1, 'y', frame_info.anot_int[1], frame_info.frame_int[1],
  2073.         frame_info.label[1], project_info.xyz_projection[1], frame_info.anot_type[1], frame_info.side[1]-1);
  2074.         
  2075. }
  2076.  
  2077. int vertical_axis (mode)
  2078. int mode; {
  2079.     /* Mode means: 1 = background axis, 2 = foreground axis, 3 = all */
  2080.     BOOLEAN go[4], fore, back;
  2081.     int i, j;
  2082.     double xp[2], yp[2];
  2083.     
  2084.     if (frame_info.anot_int[2] == 0.0) return;
  2085.  
  2086.     fore = (mode > 1);    back = (mode % 2);
  2087.     for (i = 0; i < 4; i++) go[i] = (mode == 3) ? TRUE : ((back) ? z_project.draw[i] : !z_project.draw[i]);
  2088.     
  2089.     /* Vertical */
  2090.         
  2091.     if (fore && frame_info.side[4]) xyz_axis3D (z_project.z_axis, 'z', frame_info.anot_int[2], frame_info.frame_int[2], frame_info.label[2],
  2092.         project_info.xyz_projection[2], frame_info.anot_type[2], frame_info.side[4]-1);
  2093.             
  2094.     if (frame_info.draw_box) {
  2095.         go[0] = ( (back && z_project.quadrant == 1) || (fore && z_project.quadrant != 1) );
  2096.         go[1] = ( (back && z_project.quadrant == 4) || (fore && z_project.quadrant != 4) );
  2097.         go[2] = ( (back && z_project.quadrant == 3) || (fore && z_project.quadrant != 3) );
  2098.         go[3] = ( (back && z_project.quadrant == 2) || (fore && z_project.quadrant != 2) );
  2099.         for (i = 0; i < 4; i++) {
  2100.             if (!go[i]) continue;
  2101.             geoz_to_xy (z_project.corner_x[i], z_project.corner_y[i], project_info.z_bottom, &xp[0], &yp[0]);
  2102.             geoz_to_xy (z_project.corner_x[i], z_project.corner_y[i], project_info.z_top, &xp[1], &yp[1]);
  2103.             ps_line (xp, yp, 2, 3, FALSE, TRUE);
  2104.         }
  2105.         go[0] = ( (back && (z_project.quadrant == 1 || z_project.quadrant == 4)) || (fore && (z_project.quadrant == 2 || z_project.quadrant == 3)) );
  2106.         go[1] = ( (back && (z_project.quadrant == 3 || z_project.quadrant == 4)) || (fore && (z_project.quadrant == 1 || z_project.quadrant == 2)) );
  2107.         go[2] = ( (back && (z_project.quadrant == 2 || z_project.quadrant == 3)) || (fore && (z_project.quadrant == 1 || z_project.quadrant == 4)) );
  2108.         go[3] = ( (back && (z_project.quadrant == 1 || z_project.quadrant == 2)) || (fore && (z_project.quadrant == 3 || z_project.quadrant == 4)) );
  2109.         for (i = 0; i < 4; i++) {
  2110.             if (!go[i]) continue;
  2111.             j = (i + 1) % 4;
  2112.             geoz_to_xy (z_project.corner_x[i], z_project.corner_y[i], project_info.z_top, &xp[0], &yp[0]);
  2113.             geoz_to_xy (z_project.corner_x[j], z_project.corner_y[j], project_info.z_top, &xp[1], &yp[1]);
  2114.             ps_line (xp, yp, 2, 3, FALSE, TRUE);
  2115.         }
  2116.     }
  2117.     if (back && frame_info.header[0]) {
  2118.         ps_setfont (gmtdefs.header_font);
  2119.         xp[0] = 0.5 * (z_project.xmin + z_project.xmax);
  2120.         yp[0] = z_project.ymax + 0.5;
  2121.         ps_text (xp[0], yp[0], gmtdefs.header_font_size, frame_info.header, 0.0, -2, 0);
  2122.     }
  2123. }
  2124.  
  2125. int xyz_axis3D (axis_no, axis, anotation_int, tickmark_int, label, axistype, anottype, anotate)
  2126. double anotation_int, tickmark_int;
  2127. char *label, axis;
  2128. int axis_no, axistype, anottype, anotate; {
  2129.     int i, j, i_a, i_f, k, id, justify, n_anotations = 0, n_tickmarks = 0, test;
  2130.     
  2131.     BOOLEAN do_anot, do_tick;
  2132.     
  2133.     double val, v0, v1, anot_off, label_off, start_val_a, start_val_f, end_val;
  2134.     double tvals_a[9], tvals_f[9], sign, dy, tmp, xyz[3][2], len, x0, x1, y0, y1;
  2135.     double pp[3], w[3], xp, yp, del_y, val_xyz[3], phi, size, xsize, ysize;
  2136.     double start_log_a, start_log_f, val0, val1;
  2137.     
  2138.     PFI xyz_forward, xyz_inverse;
  2139.     
  2140.     char annotation[80], format[20];
  2141.     
  2142.     id = (axis == 'x') ? 0 : ((axis == 'y') ? 1 : 2);
  2143.     j = (id == 0) ? 1 : ((id == 1) ? 0 : z_project.k);
  2144.     xyz_forward = (PFI) ((id == 0) ? x_to_xx : ((id == 1) ? y_to_yy : z_to_zz));
  2145.     xyz_inverse = (PFI) ((id == 0) ? xx_to_x : ((id == 1) ? yy_to_y : zz_to_z));
  2146.     phi = (id < 2 && axis_no > 1) ? z_project.phi[id] + 180.0 : z_project.phi[id];
  2147.     
  2148.     /* Get projected anchor point */
  2149.     
  2150.     if (id == 2) {
  2151.         geoz_to_xy (z_project.corner_x[axis_no], z_project.corner_y[axis_no], project_info.z_bottom, &x0, &y0);
  2152.         k = axis_no;
  2153.         geoz_to_xy (z_project.corner_x[axis_no], z_project.corner_y[axis_no], project_info.z_top, &x1, &y1);
  2154.         if (j == 0)
  2155.             sign = z_project.sign[z_project.z_axis];
  2156.         else
  2157.             sign = (z_project.z_axis%2) ? -z_project.sign[z_project.z_axis] : z_project.sign[z_project.z_axis];
  2158.     }
  2159.     else {
  2160.         geoz_to_xy (z_project.corner_x[axis_no], z_project.corner_y[axis_no], project_info.z_level, &x0, &y0);
  2161.         k = (axis_no + 1) % 4;
  2162.         geoz_to_xy (z_project.corner_x[k], z_project.corner_y[k], project_info.z_level, &x1, &y1);
  2163.         sign = z_project.sign[axis_no];
  2164.     }
  2165.     xyz[0][0] = project_info.w;        xyz[0][1] = project_info.e;
  2166.     xyz[1][0] = project_info.s;        xyz[1][1] = project_info.n;
  2167.     xyz[2][0] = project_info.z_bottom;    xyz[2][1] = project_info.z_top;
  2168.     
  2169.     size = gmtdefs.anot_font_size * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit];
  2170.     xsize = size * z_project.xshrink[id];
  2171.     ysize = size * z_project.yshrink[id];
  2172.     ps_comment ("Start of xyz-axis3D");
  2173.     /* Temporarely modify meaning of F0 */
  2174.     printf ("/F0 {pop /%s findfont [%lg 0 %lg %lg 0 0] makefont setfont} bind def\n",
  2175.         font_name[gmtdefs.anot_font], xsize, ysize * z_project.tilt[id], ysize);
  2176.     ps_setfont (0);
  2177.     justify = (id == 2) ? 2 : 10;
  2178.     dy = sign * gmtdefs.tick_length;
  2179.     len = (gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0;
  2180.     anotation_int = fabs (anotation_int);
  2181.     tickmark_int = fabs (tickmark_int);
  2182.     
  2183.     do_anot = (anotation_int > 0.0);
  2184.     do_tick = (tickmark_int > 0.0);
  2185.     val0 = xyz[id][0];
  2186.     val1 = xyz[id][1];
  2187.     
  2188.     
  2189.     /* Find number of decimals needed, if any */
  2190.     
  2191.     get_format (anotation_int, format);
  2192.  
  2193.     anot_off = sign * (len + gmtdefs.anot_offset);
  2194.     label_off = sign * (len + 2.5 * gmtdefs.anot_offset + (gmtdefs.anot_font_size / gmt_ppu[gmtdefs.measure_unit]) * font_height[gmtdefs.anot_font]);
  2195.     
  2196.     /* Ready to draw axis */
  2197.     
  2198.     ps_setline (gmtdefs.frame_pen);
  2199.     ps_plot (x0, y0, 3);
  2200.     ps_plot (x1, y1, -2);
  2201.     ps_setline (gmtdefs.tick_pen);
  2202.     
  2203.     i_a = i_f = 0;
  2204.     
  2205.     switch (axistype) {
  2206.         case POW:    /* Anotate in pow(x) */
  2207.             (*xyz_forward) (val0, &v0);
  2208.             (*xyz_forward) (val1, &v1);
  2209.             if (anottype == 2) {
  2210.                 val = (anotation_int == 0.0) ? 0.0 : floor (v0 / anotation_int) * anotation_int;
  2211.                 if (fabs (val - v0) > SMALL) val += anotation_int;
  2212.                 start_val_a = val;    end_val = v1;
  2213.                 val = (tickmark_int == 0.0) ? 0.0 : floor (v0 / tickmark_int) * tickmark_int;
  2214.                 if (fabs (val - v0) > SMALL) val += tickmark_int;
  2215.                 start_val_f = val;
  2216.             }
  2217.             else {
  2218.                 val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  2219.                 if (fabs (val - val0) > SMALL) val += anotation_int;
  2220.                 start_val_a = val;    end_val = val1;
  2221.                 val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  2222.                 if (fabs (val - val0) > SMALL) val += tickmark_int;
  2223.                 start_val_f = val;
  2224.             }
  2225.             break;
  2226.         case LOG10:    /* Anotate in d_log10 (x) */
  2227.             v0 = d_log10 (val0);
  2228.             v1 = d_log10 (val1);
  2229.             val = pow (10.0, floor (v0));
  2230.             test = rint (anotation_int) - 1;
  2231.             if (test < 0 || test > 2) test = 0;
  2232.             if (test == 0) {
  2233.                 n_anotations = 1;
  2234.                 tvals_a[0] = 10.0;
  2235.                 if (fabs (val - val0) > SMALL) val *= 10.0;
  2236.             }
  2237.             else if (test == 1) {
  2238.                 tvals_a[0] = 1.0;
  2239.                 tvals_a[1] = 2.0;
  2240.                 tvals_a[2] = 5.0;
  2241.                 n_anotations = 3;
  2242.             }
  2243.             else if (test == 2) {
  2244.                 n_anotations = 9;
  2245.                 for (i = 0; i < n_anotations; i++) tvals_a[i] = i + 1;
  2246.             }
  2247.             test = rint (tickmark_int) - 1;
  2248.             if (test < 0 || test > 2) test = 0;
  2249.             if (test == 0) {
  2250.                 n_tickmarks = 1;
  2251.                 tvals_f[0] = 10.0;
  2252.             }
  2253.             else if (test == 1) {
  2254.                 tvals_f[0] = 1.0;
  2255.                 tvals_f[1] = 2.0;
  2256.                 tvals_f[2] = 5.0;
  2257.                 n_tickmarks = 3;
  2258.             }
  2259.             else if (test == 2) {
  2260.                 n_tickmarks = 9;
  2261.                 for (i = 0; i < n_tickmarks; i++) tvals_f[i] = i + 1;
  2262.             }
  2263.             i_a = 0;
  2264.             start_log_a = val = pow (10.0, floor (v0));
  2265.             while ((val0 - val) > SMALL) {
  2266.                 if (i_a < n_anotations)
  2267.                     val = start_log_a * tvals_a[i_a];
  2268.                 else {
  2269.                     val = (start_log_a *= 10.0);
  2270.                     i_a = 0;
  2271.                 }
  2272.                 i_a++;
  2273.             }
  2274.             i_a--;
  2275.             start_val_a = val;
  2276.             i_f = 0;
  2277.             start_log_f = val = pow (10.0, floor (v0));
  2278.             while ((val0 - val) > SMALL) {
  2279.                 if (i_f < n_tickmarks)
  2280.                     val = start_log_f * tvals_f[i_f];
  2281.                 else {
  2282.                     val = (start_log_f *= 10.0);
  2283.                     i_f = 0;
  2284.                 }
  2285.                 i_f++;
  2286.             }
  2287.             i_f--;
  2288.             start_val_f = val;
  2289.             end_val = val1;
  2290.             break;
  2291.         case LINEAR:
  2292.             v0 = val0;
  2293.             v1 = val1;
  2294.             val = (anotation_int == 0.0) ? 0.0 : floor (val0 / anotation_int) * anotation_int;
  2295.             if (fabs (val - val0) > SMALL) val += anotation_int;
  2296.             start_val_a = val;    end_val = val1;
  2297.             val = (tickmark_int == 0.0) ? 0.0 : floor (val0 / tickmark_int) * tickmark_int;
  2298.             if (fabs (val - val0) > SMALL) val += tickmark_int;
  2299.             start_val_f = val;
  2300.             break;
  2301.     }
  2302.     
  2303.     del_y = 0.5 * sign * gmtdefs.anot_font_size * 0.732 * (justify/4) / gmt_ppu[gmtdefs.measure_unit];
  2304.     
  2305.     /* Do anotations with tickmarks */
  2306.     
  2307.     val_xyz[0] = z_project.corner_x[axis_no];
  2308.     val_xyz[1] = z_project.corner_y[axis_no];
  2309.     val_xyz[2] = project_info.z_level;
  2310.     val = (anotation_int == 0.0) ? end_val + 1.0 : start_val_a;
  2311.     while (do_anot && val <= end_val) {
  2312.     
  2313.         i_a++;
  2314.         
  2315.         val_xyz[id] = val;
  2316.         
  2317.         switch (anottype) {
  2318.             case 0:
  2319.                 sprintf (annotation, format, val_xyz[id]);
  2320.                 break;
  2321.             case 1:
  2322.                 sprintf (annotation, "%d\0", (int) d_log10 (val_xyz[id]));
  2323.                 break;
  2324.             case 2:
  2325.                 if (axistype == POW) {
  2326.                     (*xyz_inverse) (&tmp, val_xyz[id]);
  2327.                     val_xyz[id] = tmp;
  2328.                     sprintf (annotation, format, val_xyz[id]);
  2329.                 }
  2330.                 else
  2331.                     sprintf (annotation, "10@+%d@+\0", (int) d_log10 (val_xyz[id]));
  2332.                 break;
  2333.         }
  2334.         
  2335.         project3D (val_xyz[0], val_xyz[1], val_xyz[2], &w[0], &w[1], &w[2]);
  2336.         pp[0] = w[0];
  2337.         pp[1] = w[1];
  2338.         pp[2] = w[2];
  2339.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2340.         ps_plot (xp, yp, 3);
  2341.         pp[j] += dy;
  2342.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2343.         ps_plot (xp, yp, -2);
  2344.         pp[j] += anot_off -dy + del_y;
  2345.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2346.         if (anotate) {
  2347.             if (id < 2)
  2348.                 ps_text (xp, yp, gmtdefs.anot_font_size, annotation, phi, 2, 0);
  2349.             else if (val != project_info.z_level)
  2350.                 ps_text (xp, yp, gmtdefs.anot_font_size, annotation, phi, 2, 0);
  2351.         }
  2352.         
  2353.         if (axistype == LOG10) {
  2354.             if (i_a < n_anotations)
  2355.                 val = start_log_a * tvals_a[i_a];
  2356.             else {
  2357.                 val = (start_log_a *= 10.0);
  2358.                 i_a = 0;
  2359.             }
  2360.         }
  2361.         else
  2362.             val = start_val_a + i_a * anotation_int;
  2363.             
  2364.     }
  2365.  
  2366.     /* Now do frame tickmarks */
  2367.     
  2368.     dy *= 0.5;
  2369.     
  2370.     val_xyz[0] = z_project.corner_x[axis_no];
  2371.     val_xyz[1] = z_project.corner_y[axis_no];
  2372.     val_xyz[2] = project_info.z_level;
  2373.     val = (tickmark_int == 0.0) ? end_val + 1.0 : start_val_f;
  2374.     while (do_tick && val <= end_val) {
  2375.     
  2376.         i_f++;
  2377.         
  2378.         val_xyz[id] = val;
  2379.         if (anottype == 2 && axistype == POW) {
  2380.             (*xyz_inverse) (&tmp, val_xyz[id]);
  2381.             val_xyz[id] = tmp;
  2382.         }
  2383.         project3D (val_xyz[0], val_xyz[1], val_xyz[2], &w[0], &w[1], &w[2]);
  2384.                 
  2385.         pp[0] = w[0];
  2386.         pp[1] = w[1];
  2387.         pp[2] = w[2];
  2388.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2389.         ps_plot (xp, yp, 3);
  2390.         pp[j] += dy;
  2391.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2392.         ps_plot (xp, yp, -2);
  2393.         
  2394.         if (axistype == LOG10) {
  2395.             if (i_f < n_tickmarks)
  2396.                 val = start_log_f * tvals_f[i_f];
  2397.             else {
  2398.                 val = start_log_f *= 10.0;
  2399.                 i_f = 0;
  2400.             }
  2401.         }
  2402.         else
  2403.             val = start_val_f + i_f * tickmark_int;
  2404.     }
  2405.  
  2406.     /* Finally do label */
  2407.     
  2408.     if (label && anotate) {
  2409.         val_xyz[0] = z_project.corner_x[axis_no];
  2410.         val_xyz[1] = z_project.corner_y[axis_no];
  2411.         val_xyz[2] = project_info.z_level;
  2412.         size = gmtdefs.label_font_size * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit];
  2413.         xsize = size * z_project.xshrink[id];
  2414.         ysize = size * z_project.yshrink[id];
  2415.         printf ("/F0 {pop /%s findfont [%lg 0 %lg %lg 0 0] makefont setfont} bind def\n",
  2416.             font_name[gmtdefs.label_font], xsize, ysize * z_project.tilt[id], ysize);
  2417.         project3D (val_xyz[0], val_xyz[1], val_xyz[2], &w[0], &w[1], &w[2]);
  2418.         x0 = w[id];
  2419.         val_xyz[id] = (val_xyz[id] == xyz[id][0]) ? xyz[id][1] : xyz[id][0];
  2420.         project3D (val_xyz[0], val_xyz[1], val_xyz[2], &w[0], &w[1], &w[2]);
  2421.         x1 = w[id];
  2422.         pp[0] = w[0];
  2423.         pp[1] = w[1];
  2424.         pp[2] = w[2];
  2425.         pp[id] = 0.5 * (x1 + x0);
  2426.         pp[j] += label_off + del_y;
  2427.         xyz_to_xy (pp[0], pp[1], pp[2], &xp, &yp);
  2428.     
  2429.         ps_text (xp, yp, gmtdefs.label_font_size, label, phi, 2, 0);
  2430.         printf ("/F0 {/Helvetica Y} bind def\n");    /* Reset definition of F0 */
  2431.     }
  2432.     ps_comment ("End of xyz-axis3D");
  2433. }
  2434.  
  2435. int map_clip_on (r, g, b, flag)
  2436. int r, g, b, flag; {
  2437.     /* This function sets up a clip path so that only plotting
  2438.      * inside the map area will be drawn on paper. map_setup
  2439.      * must have been called first.  If r >= 0, the map area will
  2440.      * first be painted in the r,g,b colors specified.  flag can
  2441.      * be 0-3, as described in ps_clipon().
  2442.      */
  2443.      
  2444.     double *work_x, *work_y;
  2445.     int np;
  2446.     BOOLEAN donut;
  2447.     
  2448.     np = map_clip_path (&work_x, &work_y, &donut);
  2449.         
  2450.     if (donut) {
  2451.         ps_clipon (work_x, work_y, np, r, g, b, 1);
  2452.         ps_clipon (&work_x[np], &work_y[np], np, r, g, b, 2);
  2453.     }
  2454.     else
  2455.         ps_clipon (work_x, work_y, np, r, g, b, flag);
  2456.     
  2457.     free ((char *)work_x);
  2458.     free ((char *)work_y);
  2459. }
  2460.  
  2461. int map_clip_off () {
  2462.     /* Restores the original clipping path (or pours white paint) */
  2463.     
  2464.     ps_clipoff ();
  2465. }
  2466.  
  2467. int geoplot (lon, lat, pen)
  2468. double lon, lat;
  2469. int pen; {
  2470.     /* Computes x/y from lon/lat, then calls plot */
  2471.     double x, y;
  2472.     
  2473.     geo_to_xy (lon, lat, &x, &y);
  2474.     ps_plot (x, y, pen);
  2475. }
  2476.  
  2477. int timestamp (argc, argv)
  2478. int argc;
  2479. char **argv; {
  2480.     time_t right_now, time();
  2481.     int i, echo_command = FALSE;
  2482.     char *ctime(), label[512];
  2483.     double x, y, dim[5];
  2484.     
  2485.     dim[0] = 0.365; dim[1] = 0.7215; dim[2] = 0.15; dim[3] = 0.075; dim[4] = 0.1;
  2486.     if (gmtdefs.measure_unit == 0) for (i = 0; i < 5; i++) dim[i] *= 2.54;
  2487.     x = gmtdefs.unix_time_pos[0];
  2488.     y = gmtdefs.unix_time_pos[1];
  2489.     right_now = time ((time_t *)0);
  2490.     strcpy (label, ctime (&right_now));
  2491.     label[16] = 0;
  2492.     for (i = 1; i < argc && argv[i][1] != 'J'; i++);
  2493.     ps_comment (" ");
  2494.     ps_comment ("Begin time-stamp");
  2495.     ps_transrotate (x, y, 0.0);
  2496.     ps_setline (1);
  2497.     ps_rect (0.0, 0.0, dim[0]+dim[1], dim[2], gmtdefs.foreground_rgb[0], gmtdefs.foreground_rgb[1], gmtdefs.foreground_rgb[2], TRUE);
  2498.     ps_rect (0.0, 0.0, dim[0], dim[2], gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2], TRUE);
  2499.     ps_setfont (1);
  2500.     ps_setpaint (gmtdefs.foreground_rgb[0], gmtdefs.foreground_rgb[1], gmtdefs.foreground_rgb[2]);
  2501.     ps_text (0.5*dim[0], dim[3], 10, "GMT", 0.0, 6, 0);
  2502.     ps_setfont (0);
  2503.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  2504.     ps_text (dim[0]+0.5*dim[1], dim[3], 8, &label[4], 0.0, 6, 0);
  2505.     ps_setfont (1);
  2506.     label[0] = 0;
  2507.     if (gmtdefs.unix_time_label[0] == 'c' && gmtdefs.unix_time_label[1] == 0) {
  2508.         echo_command = TRUE;
  2509.         gmtdefs.unix_time_label[0] = 0;
  2510.     }
  2511.     if (echo_command) {
  2512.         strcpy (label, argv[0]);
  2513.         for (i = 1; i < argc; i++) {
  2514.             if (argv[i][0] != '-') continue;
  2515.             strcat (label, " ");
  2516.             strcat (label, argv[i]);
  2517.         }
  2518.     }
  2519.     else if (gmtdefs.unix_time_label[0])
  2520.         strcpy (label, gmtdefs.unix_time_label);
  2521.         
  2522.     if (label[0]) ps_text (dim[0]+dim[1]+dim[4], dim[3], 7, label, 0.0, 5, 0);
  2523.     ps_rotatetrans  (-x, -y, 0.0);
  2524.     ps_comment ("End time-stamp\n");
  2525. }
  2526.  
  2527. int echo_command (argc, argv)
  2528. int argc;
  2529. char **argv; {
  2530.     /* This routine will echo the command and its arguments to the
  2531.      * PostScript output file so that the user can see what scales
  2532.      * etc was used to produce this plot
  2533.      */
  2534.     int i, length = 0;
  2535.     char outstring[512];
  2536.     
  2537.     ps_comment ("PostScript produced by:\n");
  2538.     outstring[0] = 0;
  2539.     for (i = 0; i < argc; i++) {
  2540.         strcat (outstring, argv[i]);
  2541.         strcat (outstring, " ");
  2542.         length += (strlen (argv[i]) + 1);
  2543.         if (length >= 80) {
  2544.             ps_comment (outstring);
  2545.             length = 0;
  2546.             outstring[0] = 0;
  2547.         }
  2548.     }
  2549.     if (length > 0) ps_comment (outstring);
  2550.     ps_comment ("\n");
  2551. }
  2552.  
  2553. int gmt_plot_line (x, y, pen, n)
  2554. double x[], y[];
  2555. int pen[], n; {
  2556.     int i, j, i1, way, stop, close;
  2557.     double x_cross[2], y_cross[2], *xx, *yy, xt1, yt1, xt2, yt2;
  2558.     
  2559.     if (n < 2) return;
  2560.     
  2561.     i = 0;
  2562.     while (i < (n-1) && pen[i+1] == 3) i++;    /* Skip repeating pen == 3 in beginning */
  2563.     if ((n-i) < 2) return;
  2564.     while (n > 1 && pen[n-1] == 3) n--;    /* Cut off repeating pen == 3 at end */
  2565.     if ((n-i) < 2) return;
  2566.     
  2567.     for (j = i + 1; j < n && pen[j] == 2; j++);    /* j == n means no moveto's present */
  2568.     close = (j == n) ? (hypot (x[n-1] - x[i], y[n-1] - y[i]) < SMALL) : FALSE;
  2569.     
  2570.     /* First see if we can use the ps_line call directly to save points */
  2571.     
  2572.     for (j = i + 1, stop = FALSE; !stop && j < n; j++) stop = (pen[j] == 3 || map_jump (x[j-1], y[j-1], x[j], y[j]));
  2573.     if (!stop) {
  2574.         if (project_info.three_D) {    /* Must project first */
  2575.             xx = (double *) memory (CNULL, n-i, sizeof (double), "gmt_plot_line");
  2576.             yy = (double *) memory (CNULL, n-i, sizeof (double), "gmt_plot_line");
  2577.             for (j = i; j < n; j++) xy_do_z_to_xy (x[j], y[j], project_info.z_level, &xx[j], &yy[j]);
  2578.             ps_line (&xx[i], &yy[i], n - i, 3, close, TRUE);
  2579.             free ((char *)xx);
  2580.             free ((char *)yy);
  2581.         }
  2582.         else
  2583.             ps_line (&x[i], &y[i], n - i, 3, close, TRUE);
  2584.         return;
  2585.     }
  2586.     
  2587.     /* Here we must check for jumps, pen changes etc */
  2588.     
  2589.     if (project_info.three_D) {
  2590.         xy_do_z_to_xy (x[i], y[i], project_info.z_level, &xt1, &yt1);
  2591.         ps_plot (xt1, yt1, pen[i]);
  2592.     }
  2593.     else
  2594.         ps_plot (x[i], y[i], pen[i]);
  2595.  
  2596.     i++;
  2597.     while (i < n) {
  2598.         i1 = i - 1;
  2599.         if (pen[i] == pen[i1] && (way = map_jump (x[i1], y[i1], x[i], y[i]))) {    /* Jumped across the map */
  2600.             get_crossings (x_cross, y_cross, x[i1], y[i1], x[i], y[i]);
  2601.             xy_do_z_to_xy (x_cross[0], y_cross[0], project_info.z_level, &xt1, &yt1);
  2602.             xy_do_z_to_xy (x_cross[1], y_cross[1], project_info.z_level, &xt2, &yt2);
  2603.             if (project_info.three_D) {
  2604.                 xy_do_z_to_xy (xt1, yt1, project_info.z_level, &xt1, &yt1);
  2605.                 xy_do_z_to_xy (xt2, yt2, project_info.z_level, &xt2, &yt2);
  2606.             }
  2607.             if (way == -1) {    /* Add left border point */
  2608.                 ps_plot (xt1, yt1, 2);
  2609.                 ps_plot (xt2, yt2, 3);
  2610.             }
  2611.             else {
  2612.                 ps_plot (xt2, yt2, 2);
  2613.                 ps_plot (xt1, yt1, 3);
  2614.             }
  2615.             close = FALSE;
  2616.         }
  2617.         if (project_info.three_D) {
  2618.             xy_do_z_to_xy (x[i], y[i], project_info.z_level, &xt1, &yt1);
  2619.             ps_plot (xt1, yt1, pen[i]);
  2620.         }
  2621.         else
  2622.             ps_plot (x[i], y[i], pen[i]);
  2623.         i++;
  2624.     }
  2625.     if (close) ps_command ("P S") ; else ps_command ("S");
  2626. }
  2627.  
  2628. int color_image (x0, y0, x_side, y_side, image, nx, ny)
  2629. double x0, y0;        /* Lower left corner in inches */
  2630. double x_side, y_side;    /* Size of cell in inches */
  2631. unsigned char *image;    /* color image  */
  2632. int nx, ny; {        /* image size */
  2633.     
  2634.     /* Call the appropriate image filler (see pslib) */
  2635.     
  2636.     switch (gmtdefs.color_image) {
  2637.         case 0:
  2638.             ps_colorimage (x0, y0, x_side, y_side, image, nx, ny);
  2639.             break;
  2640.         case 1:
  2641.             ps_colortiles (x0, y0, x_side, y_side, image, nx, ny);
  2642.             break;
  2643.         case 2:
  2644.             ps_imagergb (x0, y0, x_side, y_side, image, nx, ny);
  2645.             break;
  2646.     }
  2647. }
  2648.  
  2649. int gmt_text3d (x, y, z, fsize, fontno, text, angle, justify, form)
  2650. double x, y, z, angle;
  2651. int fsize, fontno, justify, form;
  2652. char *text; {
  2653.     double xb, yb, xt, yt, xt1, xt2, xt3, yt1, yt2, yt3, del_y;
  2654.     double ca, sa, xshrink, yshrink, tilt, baseline_shift, size, xsize, ysize;
  2655.         if (project_info.three_D) {
  2656.                 ps_setfont (0);
  2657.                 justify = abs (justify);
  2658.                 del_y = 0.5 * fsize * 0.732 * (justify / 4) / gmt_ppu[gmtdefs.measure_unit];
  2659.                 justify %= 4;
  2660.                 ca = cosd (angle);
  2661.                 sa = sind (angle);
  2662.                 x += del_y * sa;    /* Move anchor point down on baseline */
  2663.                 y -= del_y * ca;
  2664.                 xb = x + ca;        /* Point a distance of 1.0 along baseline */
  2665.                 yb = y + sa;
  2666.                 xt = x - sa;        /* Point a distance of 1.0 normal to baseline */
  2667.                 yt = y + ca;
  2668.                 xyz_to_xy (x, y, z, &xt1, &yt1);
  2669.                 xyz_to_xy (xb, yb, z, &xt2, &yt2);
  2670.                 xyz_to_xy (xt, yt, z, &xt3, &yt3);
  2671.         xshrink = hypot (xt2-xt1, yt2-yt1) / hypot (xb-x, yb-y);    /* How lines in baseline-direction shrink */
  2672.         yshrink = hypot (xt3-xt1, yt3-yt1) / hypot (xt-x, yt-y);    /* How lines _|_ to baseline-direction shrink */
  2673.         baseline_shift = R2D * (d_atan2 (yt2 - yt1, xt2 - xt1) - d_atan2 (yb - y, xb - x));    /* Rotation of baseline */
  2674.         tilt = 90.0 - R2D * (d_atan2 (yt3 - yt1, xt3 - xt1) - d_atan2 (yt2 - yt1, xt2 - xt1));
  2675.         tilt = tand (tilt);
  2676.         size = fsize * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit];
  2677.         xsize = size * xshrink;
  2678.         ysize = size * yshrink;
  2679.         /* Temporarely modify meaning of F0 */
  2680.         printf ("/F0 {pop /%s findfont [%lg 0 %lg %lg 0 0] makefont setfont} bind def\n",
  2681.             font_name[fontno], xsize, ysize * tilt, ysize);
  2682.  
  2683.                 ps_text (xt1, yt1, fsize, text, angle + baseline_shift, justify, form);
  2684.                 printf ("/F0 {/Helvetica Y} bind def\n");       /* Reset definition of F0 */
  2685.                 ps_setfont (fontno);
  2686.         }
  2687.         else {
  2688.         ps_setfont (fontno);
  2689.         ps_text (x, y, fsize, text, angle, justify, form);
  2690.     }
  2691. }
  2692.  
  2693. int gmt_textbox3d (x, y, z, size, font, label, angle, just, outline, dx, dy, r, g, b)
  2694. double x, y, z, angle, dx, dy;
  2695. int size, font, just, r, g, b;
  2696. char *label;
  2697. BOOLEAN outline; {
  2698.         if (project_info.three_D) {
  2699.             int i, len, ndig = 0, ndash = 0, nperiod = 0;
  2700.             double xx[4], yy[4], h, w, xa, ya, cosa, sina;
  2701.         len = strlen (label);
  2702.         for (i = 0; label[i]; i++) {
  2703.             if (isdigit ((int)label[i])) ndig++;
  2704.             if (strchr (label, '.')) nperiod++;
  2705.             if (strchr (label, '-')) ndash++;
  2706.         }
  2707.         len -= (ndig + nperiod + ndash);
  2708.         w = ndig * 0.78 + nperiod * 0.38 + ndash * 0.52 + len;
  2709.         
  2710.         h = 0.58 * font_height[font] * size / gmt_ppu[gmtdefs.measure_unit];
  2711.         w *= (0.81 * h);
  2712.         just = abs (just);
  2713.         y -= (((just/4) - 1) * h);
  2714.         x -= (((just-1)%4 - 1) * w);
  2715.         xx[0] = xx[3] = -w - dx;
  2716.         xx[1] = xx[2] = w + dx;
  2717.         yy[0] = yy[1] = -h - dy;
  2718.         yy[2] = yy[3] = h + dy;
  2719.         cosa = cosd (angle);    sina = sind (angle);
  2720.         for (i = 0; i < 4; i++) {
  2721.             xa = xx[i] * cosa - yy[i] * sina;
  2722.             ya = xx[i] * sina + yy[i] * cosa;
  2723.             xx[i] = x + xa;    yy[i] = y + ya;
  2724.         }
  2725.         d_swap (z, project_info.z_level); 
  2726.         two_d_to_three_d (xx, yy, 4);
  2727.         d_swap (z, project_info.z_level);
  2728.         if (r < 0)
  2729.             ps_clipon (xx, yy, 4, r, g, b, 0);
  2730.         else
  2731.             ps_patch (xx, yy, 4, r, g, b, outline);
  2732.     }
  2733.     else
  2734.         ps_textbox (x, y, size, label, angle, just, outline, dx, dy, r, g, b);
  2735. }
  2736.  
  2737. int gmt_vector3d (x0, y0, x1, y1, z0, tailwidth, headlength, headwidth, shape, r, g, b, outline)
  2738. double x0, y0, x1, y1, z0, tailwidth, headlength, headwidth, shape;
  2739. int r, g, b;
  2740. BOOLEAN outline; {
  2741.     int i;
  2742.     double xx[7], yy[7], dx, dy, angle, length;
  2743.     
  2744.     if (project_info.three_D) {
  2745.         angle = atan2 (y1 - y0, x1 - x0);
  2746.         length = hypot (y1 - y0, x1 - x0);
  2747.         xx[3] = x0 + length * cos (angle);
  2748.         yy[3] = y0 + length * sin (angle);
  2749.         dx = 0.5 * tailwidth * sin (angle);
  2750.         dy = 0.5 * tailwidth * cos (angle);
  2751.         xx[0] = x0 + dx;    xx[6] = x0 - dx;
  2752.         yy[0] = y0 - dy;    yy[6] = y0 + dy;
  2753.         dx = (length - (1.0 - 0.5 * shape) * headlength) * cos (angle);
  2754.         dy = (length - (1.0 - 0.5 * shape) * headlength) * sin (angle);
  2755.         xx[1] = xx[0] + dx;    xx[5] = xx[6] + dx;
  2756.         yy[1] = yy[0] + dy;    yy[5] = yy[6] + dy;
  2757.         x0 += (length - headlength) * cos (angle);
  2758.         y0 += (length - headlength) * sin (angle);
  2759.         dx = headwidth * sin (angle);
  2760.         dy = headwidth * cos (angle);
  2761.         xx[2] = x0 + dx;    xx[4] = x0 - dx;
  2762.         yy[2] = y0 - dy;    yy[4] = y0 + dy;
  2763.         for (i = 0; i < 7; i++) {
  2764.             xyz_to_xy (xx[i], yy[i], z0, &x0, &y0);
  2765.             xx[i] = x0;
  2766.             yy[i] = y0;
  2767.         }
  2768.         ps_polygon (xx, yy, 7, r, g, b, outline);
  2769.     }
  2770.     else
  2771.         ps_vector (x0, y0, x1, y1, tailwidth, headlength, headwidth, gmtdefs.vector_shape, r, g, b, outline);
  2772. }
  2773.  
  2774. int prepare_label (angle, side, x, y, type, line_angle, text_angle, justify)
  2775. double angle, x, y, *line_angle, *text_angle;
  2776. int side, type, *justify; {
  2777.     BOOLEAN set_angle;
  2778.         
  2779.     if (!project_info.edge[side]) return -1;        /* Side doesnt exist */
  2780.     if (frame_info.side[side] < 2) return -1;    /* Dont want labels here */
  2781.         
  2782.     if (frame_info.check_side == TRUE) {
  2783.         if (type == 0 && side%2) return -1;
  2784.         if (type == 1 && !(side%2)) return -1;
  2785.     }
  2786.     
  2787.     if (frame_info.horizontal == 2 && !(side%2)) angle = -90.0;    /* get_label_parameters will make this 0 */
  2788.     
  2789.     if (angle < 0.0) angle += 360.0;
  2790.     
  2791.     set_angle = ((project_info.region && !(AZIMUTHAL || CONICAL)) || !project_info.region);
  2792.     if (set_angle) {
  2793.         if (side == 0 && angle < 180.0) angle -= 180.0;
  2794.         if (side == 1 && (angle > 90.0 && angle < 270.0)) angle -= 180.0;
  2795.         if (side == 2 && angle > 180.0) angle -= 180.0;
  2796.         if (side == 3 && (angle < 90.0 || angle > 270.0)) angle -= 180.0;
  2797.     }
  2798.     
  2799.     if (!get_label_parameters (side, angle, text_angle, justify)) return -1;
  2800.     *line_angle = angle;
  2801.     
  2802.     if (!set_angle) *justify = polar_adjust (side, angle, x, y);
  2803.     
  2804.     return 0;
  2805. }
  2806.  
  2807. int get_anot_label (val, label, do_minutes, do_seconds, lonlat)
  2808. double val;    /* Degree value of anotation */
  2809. char *label;    /* String to hold the final anotation */
  2810. int do_minutes;    /* TRUE if degree and minutes are desired, FALSE for just integer degrees */
  2811. int do_seconds;    /* TRUE if degree, minutes, and seconds are desired */
  2812. int lonlat; {    /* 0 = longitudes, 1 = latitudes */
  2813.     int ival, minutes, seconds;
  2814.     BOOLEAN zero_fix = FALSE;
  2815.     char letter = 0;
  2816.     
  2817.     if (lonlat == 0) {    /* Fix longitudes to 0.0 <= lon <= 360 first */
  2818.         while (val > 360.0) val -= 360.0;
  2819.         while (val < 0.0) val += 360.0;
  2820.     }
  2821.  
  2822.     if (val == 360.0 && !gmt_world_map) val = 0.0;
  2823.     if (val == 360.0 && gmt_world_map && project_info.projection == OBLIQUE_MERC) val = 0.0;
  2824.  
  2825.     switch (gmtdefs.degree_format) {
  2826.         case 0:    /* Use 0 to 360 for longitudes and -90 to +90 for latitudes */
  2827.             break;
  2828.         case 1:    /* Use -180 to +180 for longitudes and -90 to +90 for latitudes */
  2829.             if (lonlat == 0 && val > 180.0) val -= 360.0;
  2830.             break;
  2831.         case 2:    /* Use unsigned 0 to 180 for longitudes and 0 to 90 for latitudes */
  2832.             if (lonlat == 0 && val > 180.0) val -= 360.0;
  2833.             val = fabs (val);
  2834.             break;
  2835.         case 3:    /* Use 0 to 180E and 0 to 180W for longitudes and 90S to 90N for latitudes */
  2836.             if (lonlat == 0) {
  2837.                 if (val > 180.0) val -= 360.0;
  2838.                 letter = (val == 0.0 || val == 180.0) ? 0 : ((val < 0.0) ? 'W' : 'E');
  2839.             }
  2840.             else
  2841.                 letter = (val == 0.0) ? 0 : ((val < 0.0) ? 'S' : 'N');
  2842.             val = fabs (val);
  2843.             break;
  2844.         case 4:    /* Use decimal 0 to 360 for longitudes and -90 to +90 for latitudes */
  2845.             break;
  2846.         case 5:    /* Use decimal -180 to +180 for longitudes and -90 to +90 for latitudes */
  2847.             if (lonlat == 0 && val > 180.0) val -= 360.0;
  2848.             break;
  2849.     }
  2850.     
  2851.     if (gmtdefs.degree_format == 6) {    /* theta-r */
  2852.         (lonlat) ? sprintf (label, "%lg\0", val) : sprintf (label, "%lg\\312\0", val);
  2853.         return;
  2854.     }
  2855.     
  2856.     if (gmtdefs.degree_format > 3) {
  2857.         sprintf (label, "%lg\\312\0", val);
  2858.         return;
  2859.     }
  2860.         
  2861.     ival = val;    /* Truncate to integer in the direction toward 0 */
  2862.     minutes = seconds = 0;
  2863.     if (fabs (val - (double) ival) > SMALL) {
  2864.         minutes = floor (fabs ((val - ival) * 60.0) + SMALL);
  2865.         if (minutes == 60) {
  2866.             minutes = 0;
  2867.             ival = rint (val);
  2868.         }
  2869.         seconds = rint (fabs ((val - ival - minutes / 60.0) * 3600.0));
  2870.         if (seconds == 60) {
  2871.             seconds = 0;
  2872.             minutes++;
  2873.             if (minutes == 60) {
  2874.                 minutes = 0;
  2875.                 ival = rint (val);
  2876.             }
  2877.         }
  2878.     }
  2879.     
  2880.     if (do_minutes) {
  2881.         if (ival == 0 && val < 0.0) {    /* Must write out -0 degrees, do so by writing -1 and change 1 to 0 */
  2882.             ival = -1;
  2883.             zero_fix = TRUE;
  2884.         }
  2885.         if (do_seconds) {
  2886.             if (letter)
  2887.                 sprintf (label, "%d\\312 %.2d\\251 %.2d\\042%c\0", ival, minutes, seconds, letter);
  2888.             else
  2889.                 sprintf (label, "%d\\312 %.2d\\251 %.2d\\042\0", ival, minutes, seconds);
  2890.         }
  2891.         else {
  2892.             if (letter)
  2893.                 sprintf (label, "%d\\312 %.2d\\251%c\0", ival, minutes, letter);
  2894.             else
  2895.                 sprintf (label, "%d\\312 %.2d\\251\0", ival, minutes);
  2896.         }
  2897.         if (zero_fix) label[1] = '0';    /* Undo the fix above */
  2898.     }
  2899.     else {
  2900.         if (letter)
  2901.             sprintf (label, "%d\\312%c\0", ival, letter);
  2902.         else
  2903.             sprintf (label, "%d\\312\0", ival);
  2904.     }
  2905.     return;
  2906. }
  2907.  
  2908. int polar_adjust (side, angle, x, y)
  2909. int side;
  2910. double angle, x, y; {
  2911.     int justify, left, right, top, bottom, low;
  2912.     double x0, y0;
  2913.     
  2914.     geo_to_xy (project_info.central_meridian, project_info.pole, &x0, &y0);
  2915.     if (project_info.north_pole) {
  2916.         low = 0;
  2917.         left = 7;
  2918.         right = 5;
  2919.     }
  2920.     else {
  2921.         low = 2;
  2922.         left = 5;
  2923.         right = 7;
  2924.     }
  2925.     if (y >= y0) {
  2926.         top = 2;
  2927.         bottom = 10;
  2928.     }
  2929.     else {
  2930.         top = 10;
  2931.         bottom = 2;
  2932.     }
  2933.     
  2934.     if (side%2) {    /* W and E border */
  2935.         if (y >= y0)
  2936.             justify = (side == 1) ? left : right;
  2937.         else
  2938.             justify = (side == 1) ? right : left;
  2939.     }
  2940.     else {
  2941.         if (frame_info.horizontal) {
  2942.             if (side == low)
  2943.                 justify = (angle == 180.0) ? bottom : top;
  2944.             else
  2945.                 justify = (angle == 0.0) ? top : bottom;
  2946.         }    
  2947.         else {
  2948.             if (x >= x0)
  2949.                 justify = (side == 2) ? left : right;
  2950.             else
  2951.                 justify = (side == 2) ? right : left;
  2952.         }
  2953.     }
  2954.     return (justify);
  2955. }
  2956.  
  2957. double get_angle (lon1, lat1, lon2, lat2)
  2958. double lon1, lat1, lon2, lat2; {
  2959.     double x1, y1, x2, y2, angle, direction;
  2960.     
  2961.     geo_to_xy (lon1, lat1, &x1, &y1);
  2962.     geo_to_xy (lon2, lat2, &x2, &y2);
  2963.     angle = d_atan2 (y2-y1, x2-x1) * R2D;
  2964.     
  2965.     if (abs (gmt_x_status_old) == 2 && abs (gmt_y_status_old) == 2)    /* Last point outside */
  2966.         direction = angle + 180.0;
  2967.     else if (gmt_x_status_old == 0 && gmt_y_status_old == 0)        /* Last point inside */
  2968.         direction = angle;
  2969.     else {
  2970.         if (abs (gmt_x_status_new) == 2 && abs (gmt_y_status_new) == 2)    /* This point outside */
  2971.             direction = angle;
  2972.         else if (gmt_x_status_new == 0 && gmt_y_status_new == 0)        /* This point inside */
  2973.             direction = angle + 180.0;
  2974.         else {    /* Special case of corners and sides only */
  2975.             if (gmt_x_status_old == gmt_x_status_new)
  2976.                 direction = (gmt_y_status_old == 0) ? angle : angle + 180.0;
  2977.             else if (gmt_y_status_old == gmt_y_status_new)
  2978.                 direction = (gmt_x_status_old == 0) ? angle : angle + 180.0;
  2979.             else
  2980.                 direction = angle;
  2981.             
  2982.         }
  2983.     }
  2984.     
  2985.     if (direction < 0.0) direction += 360.0;
  2986.     if (direction >= 360.0) direction -= 360.0;
  2987.     return (direction);
  2988. }
  2989.  
  2990.  
  2991. int draw_map_scale (x0, y0, lat, length, miles, gave_xy, fancy)
  2992. double x0, y0, lat, length;
  2993. BOOLEAN miles, gave_xy, fancy; {
  2994.     int i, j, k, *rgb, n_a_ticks[9], n_f_ticks[9];
  2995.     double dlon, x1, x2, dummy, a, b, tx, ty, off, f_len, a_len, x_left, bar_length;
  2996.     double xx[4], yy[4], bx[4], by[4], base, d_base, width, half, bar_width, dx_f, dx_a;
  2997.     char txt[80];
  2998.     
  2999.     if (!MAPPING) return;    /* Only for geographic projections */
  3000.     
  3001.     if (!gave_xy) {
  3002.         geo_to_xy (x0, y0, &a, &b);
  3003.         x0 = a;
  3004.         y0 = b;
  3005.     }
  3006.     
  3007.     bar_length = (miles) ? 1.609344 * length : length;
  3008.     dlon = 0.5 * bar_length * 1000.0 / (M_PR_DEG * cosd (lat));
  3009.     
  3010.     geoz_to_xy (project_info.central_meridian - dlon, lat, project_info.z_level, &x1, &dummy);
  3011.     geoz_to_xy (project_info.central_meridian + dlon, lat, project_info.z_level, &x2, &dummy);
  3012.     width = x2 - x1;
  3013.     half = 0.5 * width;
  3014.     a_len = fabs (gmtdefs.tick_length);
  3015.     off = a_len + 0.75 * gmtdefs.anot_offset;
  3016.     
  3017.         ps_setline (gmtdefs.tick_pen);
  3018.     if (fancy) {    /* Fancy scale */
  3019.         n_f_ticks[8] = 3;
  3020.         n_f_ticks[1] = n_f_ticks[3] = n_f_ticks[7] = 4;
  3021.         n_f_ticks[0] = n_f_ticks[4] = 5;
  3022.         n_f_ticks[2] = n_f_ticks[5] = 6;
  3023.         n_f_ticks[6] = 7;
  3024.         n_a_ticks[4] = n_a_ticks[6] = n_a_ticks[8] = 1;
  3025.         n_a_ticks[0] = n_a_ticks[1] = n_a_ticks[3] = n_a_ticks[7] = 2;
  3026.         n_a_ticks[2] = n_a_ticks[5] = 3;
  3027.         base = pow (10.0, floor (d_log10 (length)));
  3028.         i = rint (length / base) - 1;
  3029.         d_base = length / n_a_ticks[i];
  3030.         dx_f = width / n_f_ticks[i];
  3031.         dx_a = width / n_a_ticks[i];
  3032.         bar_width = 0.5 * fabs (gmtdefs.tick_length);
  3033.         f_len = 0.75 * fabs (gmtdefs.tick_length);
  3034.         yy[2] = yy[3] = y0;
  3035.         yy[0] = yy[1] = y0 - bar_width;
  3036.         x_left = x0 - half;
  3037.         xyz_to_xy (x_left, y0 - f_len, project_info.z_level, &a, &b);
  3038.         ps_plot (a, b, 3);
  3039.         xyz_to_xy (x_left, y0, project_info.z_level, &a, &b);
  3040.         ps_plot (a, b, 2);
  3041.         for (j = 0; j < n_f_ticks[i]; j++) {
  3042.             xx[0] = xx[3] = x_left + j * dx_f;
  3043.             xx[1] = xx[2] = xx[0] + dx_f;
  3044.             for (k = 0; k < 4; k++) xyz_to_xy (xx[k], yy[k], project_info.z_level, &bx[k], &by[k]);
  3045.             rgb = (j%2) ? gmtdefs.foreground_rgb : gmtdefs.background_rgb;
  3046.             ps_polygon (bx, by, 4, rgb[0], rgb[1], rgb[2], TRUE);
  3047.             xyz_to_xy (xx[1], y0 - f_len, project_info.z_level, &a, &b);
  3048.             ps_plot (a, b, 3);
  3049.             xyz_to_xy (xx[1], y0, project_info.z_level, &a, &b);
  3050.             ps_plot (a, b, 2);
  3051.         }
  3052.         ty = y0 - off;
  3053.         for (j = 0; j <= n_a_ticks[i]; j++) {
  3054.             tx = x_left + j * dx_a;
  3055.             xyz_to_xy (tx, y0 - a_len, project_info.z_level, &a, &b);
  3056.             ps_plot (a, b, 3);
  3057.             xyz_to_xy (tx, y0, project_info.z_level, &a, &b);
  3058.             ps_plot (a, b, 2);
  3059.             sprintf (txt, "%lg\0", j * d_base);
  3060.             gmt_text3d (tx, ty, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, txt, 0.0, 10, 0);
  3061.         }
  3062.         if (miles)
  3063.             strcpy (txt, "miles");
  3064.         else
  3065.             strcpy (txt, "km");
  3066.         xyz_to_xy (x0, y0 + f_len, project_info.z_level, &tx, &ty);
  3067.         gmt_text3d (tx, ty, project_info.z_level, gmtdefs.label_font_size, gmtdefs.label_font, txt, 0.0, 2, 0);
  3068.     }
  3069.     else {    /* Simple scale */
  3070.     
  3071.         if (miles)
  3072.             sprintf (txt, "%lg miles\0", length);
  3073.         else
  3074.             sprintf (txt, "%lg km\0", length);
  3075.         xyz_to_xy (x0 - half, y0 - gmtdefs.tick_length, project_info.z_level, &a, &b);
  3076.         ps_plot (a, b, 3);
  3077.         xyz_to_xy (x0 - half, y0, project_info.z_level, &a, &b);
  3078.         ps_plot (a, b, 2);
  3079.         xyz_to_xy (x0 + half, y0, project_info.z_level, &a, &b);
  3080.         ps_plot (a, b, 2);
  3081.         xyz_to_xy (x0 + half, y0 - gmtdefs.tick_length, project_info.z_level, &a, &b);
  3082.         ps_plot (a, b, 2);
  3083.         gmt_text3d (x0, y0 - off, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, txt, 0.0, 10, 0);
  3084.     }
  3085. }
  3086.