home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_support.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-28  |  72.5 KB  |  2,571 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)gmt_support.c    2.49  28 Jul 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  *
  9.  *            G M T _ S U P P O R T . C
  10.  *
  11.  *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  12.  * gmt_support.c contains code used by most GMT programs
  13.  *
  14.  * Author:    Paul Wessel
  15.  * Date:    27-JUL-1992
  16.  * Version:    2.1
  17.  *
  18.  * Modules in this file:
  19.  *
  20.  *    akima            Akima's 1-D spline
  21.  *    check_rgb        Check rgb for valid range
  22.  *    check_region        Check region in -R for valid range
  23.  *    comp_double_asc        Used when sorting doubles into ascending order [checks for NaN]
  24.  *    comp_float_asc        Used when sorting floats into ascending order [checks for NaN]
  25.  *    contours        Subroutine for contouring
  26.  *    cspline/csplint        Natural cubic 1-D spline
  27.  *    bcr_init        Initialize structure for bicubic interpolation
  28.  *    delaunay        Performs a Delaunay triangulation
  29.  *    get_bcr_z        Get bicubic interpolated value
  30.  *      get_bcr_nodal_values    Supports -"-
  31.  *      get_bcr_cardinals       "
  32.  *      get_bcr_ij           "
  33.  *      get_bcr_xy           "
  34.  *    get_index        Return color table entry for given z
  35.  *    get_format :        Find # of decimals and create format string
  36.  *    get_rgb24        Return rgb for given z
  37.  *    get_plot_array        Allocate memory for plotting arrays
  38.  *    free_plot_array        Free such memory
  39.  *    gmt_getfill        Decipher and check fill argument
  40.  *    gmt_getinc        Decipher and check increment argument
  41.  *    gmt_getpen        Decipher and check pen argument
  42.  *    gmt_getrgb        Dechipher and check color argument
  43.  *    gmt_init_fill        Initialize fill attributes
  44.  *    gmt_init_pen        Initialize pen attributes
  45.  *    grd_init :        Initialize grd header structure
  46.  *    grd_shift :        Rotates grdfiles in x-direction
  47.  *    grd_setregion :        Determines subset coordinates for grdfiles
  48.  *    gmt_hsv_to_rgb        Convert HSV to RGB
  49.  *    illuminate        Add illumination effects to rgb
  50.  *    intpol            1-D interpolation
  51.  *    memory            Memory allocation/reallocation
  52.  *    median            Finds the median without sorting
  53.  *    non_zero_winding    Finds if a point is inside/outside a polygon
  54.  *    plane            Find best-fitting plane to grd matrix
  55.  *    read_cpt        Read color palette file
  56.  *    gmt_rgb_to_hsv        Convert RGB to HSV
  57.  *    set_bcr_boundaries    Sets natural boundary conditions at grid edges
  58.  *    smooth_contour        Use Akima's spline to smooth contour
  59.  *    start_trace        Subfunction used by trace_contour
  60.  *    trace_contour        Function that trace the contours in contours
  61.  */
  62.  
  63. #include "gmt.h"
  64. #include "gmt_bcr.h"
  65.  
  66. #define I_255    (1.0 / 255.0)
  67.  
  68. int check_rgb (r, g, b)
  69. int r, g, b; {
  70.     return (( (r < 0 || r > 255) || (g < 0 || g > 255) || (b < 0 || b > 255) ));
  71. }
  72.  
  73. int check_region (w, e, s, n)
  74. double w, e, s, n; {    /* If region is given then we must have w < e and s < n */
  75.     return ((w >= e || s >= n) && project_info.region);
  76. }
  77.  
  78. int gmt_init_fill (fill, r, g, b)
  79. struct FILL *fill;
  80. int r, g, b; {    /* Initialize FILL structure */
  81.     fill->use_pattern = fill->inverse = FALSE;
  82.     fill->pattern[0] = 0;
  83.     fill->pattern_no = 0;
  84.     fill->icon_size = 0.0;
  85.     fill->r = r;    fill->g = g; fill->b = b;
  86. }
  87.  
  88. int gmt_getfill (line, fill)
  89. char *line;
  90. struct FILL *fill; {
  91.     int n, count, error = FALSE, slash_count();
  92.     
  93.     if (line[0] == 'p' || line[0] == 'P') {    /* Pattern specified */
  94.         n = sscanf (&line[1], "%lf/%s", &fill->icon_size, fill->pattern);
  95.         if (n != 2) error = TRUE;
  96.         fill->pattern_no = atoi (fill->pattern);
  97.         if (fill->pattern_no == 0) fill->pattern_no = -1;
  98.         fill->inverse = !(line[0] == 'P');
  99.         fill->use_pattern = TRUE;
  100.     }
  101.     else {    /* Plain color or shade */
  102.         if ((count = slash_count (line)) == 2)
  103.             n = sscanf (line, "%d/%d/%d", &fill->r, &fill->g, &fill->b);
  104.         else if (count == 0) {
  105.             n = sscanf (line, "%d", &fill->r);
  106.             fill->g = fill->b = fill->r;
  107.         }
  108.         else
  109.             n = 0;
  110.         fill->use_pattern = FALSE;
  111.         error = (n < 1 || n > 3) ? TRUE : check_rgb (fill->g, fill->b, fill->r);
  112.     }
  113.     return (error);
  114. }
  115.  
  116. int slash_count (txt)
  117. char *txt; {
  118.     int i = 0, n = 0;
  119.     while (txt[i]) if (txt[i++] == '/') n++;
  120.     return (n);
  121. }
  122.  
  123. int gmt_getrgb (line, r, g, b)
  124. char *line;
  125. int *r, *g, *b; {
  126.     int n, count, slash_count();
  127.     
  128.     if ((count = slash_count (line)) == 2)
  129.         n = sscanf (line, "%d/%d/%d", r, g, b);
  130.     else if (count == 0) {
  131.         n = sscanf (line, "%d", r);
  132.         *g = *b = *r;
  133.     }
  134.     else
  135.         n = 0;
  136.     return (n < 1 || n > 3 || check_rgb (*r, *g, *b));
  137. }
  138.  
  139. int gmt_init_pen (pen, width)
  140. struct PEN *pen;
  141. int width; {
  142.     pen->width = width;
  143.     pen->r = gmtdefs.basemap_frame_rgb[0];
  144.     pen->g = gmtdefs.basemap_frame_rgb[1];
  145.     pen->b = gmtdefs.basemap_frame_rgb[2];
  146.     pen->texture[0] = 0;
  147.     pen->offset = 0;
  148. }
  149.  
  150. int gmt_getpen (line, pen)
  151. char *line;
  152. struct PEN *pen; {
  153.     int i, n_slash, t_pos, c_pos, s_pos, width;
  154.     BOOLEAN get_pen, points = FALSE;
  155.     double val = 0.0;
  156.     
  157.     /* Check for p which means pen thickness is given in points */
  158.     
  159.     points = (BOOLEAN) strchr (line, 'p');
  160.     
  161.     /* First check if a pen thickness has been entered */
  162.     
  163.     get_pen = ((line[0] == '-' && (line[1] >= '0' && line[1] <= '9')) || (line[0] >= '0' && line[0] <= '9'));
  164.     
  165.     /* Then look for slash which indicate start of color information */
  166.     
  167.     for (i = n_slash = 0, s_pos = -1; line[i]; i++) if (line[i] == '/') {
  168.         n_slash++;
  169.         if (s_pos < 0) s_pos = i;    /* First slash position */
  170.     }
  171.     
  172.     /* Finally see if a texture is given */
  173.     
  174.     for (i = 0, t_pos = -1; line[i] && t_pos == -1; i++) if (line[i] == 't') t_pos = i;
  175.     
  176.     /* Now decode values */
  177.     
  178.     if (get_pen) {    /* Decode pen width */
  179.         val = atof (line);
  180.         pen->width = (points) ? rint (val * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit]) : rint (val);
  181.     }
  182.     if (s_pos >= 0) {    /* Get color of pen */
  183.         s_pos++;
  184.         if (n_slash == 1) {    /* Gray shade is given */
  185.             sscanf (&line[s_pos], "%d", &pen->r);
  186.             pen->g = pen->b = pen->r;
  187.         }
  188.         else if (n_slash == 3)    /* Color is given */
  189.             sscanf (&line[s_pos], "%d/%d/%d", &pen->r, &pen->g, &pen->b);
  190.     }
  191.     
  192.     if (t_pos >= 0) {    /* Get texture */
  193.         t_pos++;
  194.         if (line[t_pos] == 'o') {    /* Dotted */
  195.             width = (pen->width == 0) ? 1 : pen->width;
  196.             sprintf (pen->texture, "%d %d\0", width, 4 * width);
  197.             pen->offset = 0;
  198.         }
  199.         else if (line[t_pos] == 'a') {    /* Dashed */
  200.             width = (pen->width == 0) ? 1 : pen->width;
  201.             sprintf (pen->texture, "%d %d\0", 8 * width, 8 * width);
  202.             pen->offset = 4 * width;
  203.         }
  204.         else {    /* Specified pattern */
  205.             for (i = t_pos+1, c_pos = 0; line[i] && c_pos == 0; i++) if (line[i] == ':') c_pos = i;
  206.             if (c_pos == 0) return (pen->width < 0 || check_rgb (pen->r, pen->g, pen->b));
  207.             line[c_pos] = ' ';
  208.             sscanf (&line[t_pos], "%s %d", pen->texture, &pen->offset);
  209.             line[c_pos] = ':';
  210.             for (i = 0; pen->texture[i]; i++) if (pen->texture[i] == '_') pen->texture[i] = ' ';
  211.             if (points) {    /* Must convert points to current dpi */
  212.                 char tmp[10], string[BUFSIZ], *ptr;
  213.                 
  214.                 ptr = strtok (pen->texture, " ");
  215.                 memset (string, 0, BUFSIZ);
  216.                 while (ptr) {
  217.                     sprintf (tmp, "%d \0", (int) rint (atof (ptr) * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit]));
  218.                     strcat (string, tmp);
  219.                     ptr = strtok (CNULL, " ");
  220.                 }
  221.                 string[strlen (string) - 1] = 0;
  222.                 strcpy (pen->texture, string);
  223.                 pen->offset = rint (pen->offset * gmtdefs.dpi / gmt_ppu[gmtdefs.measure_unit]);
  224.             }
  225.                     
  226.         }
  227.     }
  228.     return (pen->width < 0 || check_rgb (pen->r, pen->g, pen->b));
  229. }
  230.  
  231. int gmt_getinc (line, dx, dy)
  232. char *line;
  233. double *dx, *dy; {
  234.     int t_pos, i;
  235.     char xstring[80], ystring[80];
  236.     
  237.     for (t_pos = -1, i = 0; line[i] && t_pos < 0; i++) if (line[i] == '/') t_pos = i;
  238.     
  239.     if (t_pos != -1) {
  240.         strcpy (xstring, line);
  241.         xstring[t_pos] = 0;
  242.         if (t_pos > 0 && (xstring[t_pos-1] == 'm' || xstring[t_pos-1] == 'M') ) {
  243.             xstring[t_pos-1] = 0;
  244.             if ( (sscanf(xstring, "%lf", dx)) != 1) return(1);
  245.             (*dx) /= 60.0;
  246.         }
  247.         else if (t_pos > 0 && (xstring[t_pos-1] == 'c' || xstring[t_pos-1] == 'C') ) {
  248.             xstring[t_pos-1] = 0;
  249.             if ( (sscanf(xstring, "%lf", dx)) != 1) return(1);
  250.             (*dx) /= 3600.0;
  251.         }
  252.         else {
  253.             if ( (sscanf(xstring, "%lf", dx)) != 1) return(1);
  254.         }
  255.  
  256.         strcpy (ystring, &line[t_pos+1]);
  257.         t_pos = strlen(ystring);
  258.         if (t_pos > 0 && (ystring[t_pos-1] == 'm' || ystring[t_pos-1] == 'M') ) {
  259.             ystring[t_pos-1] = 0;
  260.             if ( (sscanf(ystring, "%lf", dy)) != 1) return(1);
  261.             (*dy) /= 60.0;
  262.         }
  263.         else if (t_pos > 0 && (ystring[t_pos-1] == 'c' || ystring[t_pos-1] == 'C') ) {
  264.             ystring[t_pos-1] = 0;
  265.             if ( (sscanf(ystring, "%lf", dy)) != 1) return(1);
  266.             (*dy) /= 3600.0;
  267.         }
  268.         else {
  269.             if ( (sscanf(ystring, "%lf", dy)) != 1) return(1);
  270.         }
  271.     }
  272.     else {
  273.         t_pos = strlen(line);
  274.         if (t_pos > 0 && (line[t_pos-1] == 'm' || line[t_pos-1] == 'M') ) {
  275.             line[t_pos-1] = 0;
  276.             if ( (sscanf(line, "%lf", dx)) != 1) return(1);
  277.             (*dx) /= 60.0;
  278.         }
  279.         else if (t_pos > 0 && (line[t_pos-1] == 'c' || line[t_pos-1] == 'C') ) {
  280.             line[t_pos-1] = 0;
  281.             if ( (sscanf(line, "%lf", dx)) != 1) return(1);
  282.             (*dx) /= 3600.0;
  283.         }
  284.         else {
  285.             if ( (sscanf(line, "%lf", dx)) != 1) return(1);
  286.         }
  287.         *dy = (*dx);
  288.     }
  289.     return(0);
  290. }
  291.  
  292. int read_cpt (cpt_file)
  293. char *cpt_file; {
  294.     /* Opens and reads a color palette file in RGB or HSV of arbitrary length */
  295.     
  296.     int n = 0, i, nread, n_alloc = GMT_SMALL_CHUNK;
  297.     double r0, r1, g0, g1, b0, b1;
  298.     BOOLEAN gap, error = FALSE;
  299.     char line[512], c;
  300.     FILE *fp = NULL;
  301.     
  302.     if ((fp = fopen (cpt_file, "r")) == NULL) {
  303.         fprintf (stderr, "GMT Fatal Error: Cannot open color palette table %s\n", cpt_file);
  304.         exit (-1);
  305.     }
  306.     
  307.     gmt_lut = (struct GMT_LUT *) memory (CNULL, n_alloc, sizeof (struct GMT_LUT), "read_cpt");
  308.     
  309.     gmt_b_and_w = gmt_gray = TRUE;
  310.     gmt_continuous = FALSE;
  311.     
  312.     while (!error && fgets (line, 512, fp)) {
  313.     
  314.         if (line[0] == '#') continue;    /* Comment */
  315.         
  316.         if (line[0] == 'F' || line[0] == 'f') {    /* Foreground color */
  317.             nread = sscanf (&line[1], "%lf %lf %lf", &r0, &g0, &b0);
  318.             if (!(nread == 1 || nread == 3)) error = TRUE;
  319.             if (nread == 1 || gmtdefs.color_model == 0) {    /* Read r,g,b or gray */
  320.                 if (nread == 1)    /* Grayshades given */
  321.                     gmtdefs.foreground_rgb[0] = gmtdefs.foreground_rgb[1] = gmtdefs.foreground_rgb[2] = r0;
  322.                 else {
  323.                     gmtdefs.foreground_rgb[0] = r0;
  324.                     gmtdefs.foreground_rgb[1] = g0;
  325.                     gmtdefs.foreground_rgb[2] = b0;
  326.                 }
  327.             }
  328.             else
  329.                 gmt_hsv_to_rgb (&gmtdefs.foreground_rgb[0], &gmtdefs.foreground_rgb[1], &gmtdefs.foreground_rgb[2], r0, g0, b0);
  330.             if (!(gmtdefs.foreground_rgb[0] == gmtdefs.foreground_rgb[1] && gmtdefs.foreground_rgb[0] == gmtdefs.foreground_rgb[2])) gmt_gray = FALSE;
  331.             if (gmt_gray && !(gmtdefs.foreground_rgb[0] == 0 || gmtdefs.foreground_rgb[0] == 255)) gmt_b_and_w = FALSE;
  332.             continue;
  333.         }
  334.         else if (line[0] == 'B' || line[0] == 'b') {    /* Background color */
  335.             nread = sscanf (&line[1], "%lf %lf %lf", &r0, &g0, &b0);
  336.             if (!(nread == 1 || nread == 3)) error = TRUE;
  337.             if (nread == 1 || gmtdefs.color_model == 0) {    /* Read r,g,b  or gray*/
  338.                 if (nread == 1)    /* Grayshades given */
  339.                     gmtdefs.background_rgb[0] = gmtdefs.background_rgb[1] = gmtdefs.background_rgb[2] = r0;
  340.                 else {
  341.                     gmtdefs.background_rgb[0] = r0;
  342.                     gmtdefs.background_rgb[1] = g0;
  343.                     gmtdefs.background_rgb[2] = b0;
  344.                 }
  345.             }
  346.             else
  347.                 gmt_hsv_to_rgb (&gmtdefs.background_rgb[0], &gmtdefs.background_rgb[1], &gmtdefs.background_rgb[2], r0, g0, b0);
  348.             if (!(gmtdefs.background_rgb[0] == gmtdefs.background_rgb[1] && gmtdefs.background_rgb[0] == gmtdefs.background_rgb[2])) gmt_gray = FALSE;
  349.             if (gmt_gray && !(gmtdefs.background_rgb[0] == 0 || gmtdefs.background_rgb[0] == 255)) gmt_b_and_w = FALSE;
  350.             continue;
  351.         }
  352.         else if (line[0] == 'N' || line[0] == 'n') {    /* NaN color */
  353.             nread = sscanf (&line[1], "%lf %lf %lf", &r0, &g0, &b0);
  354.             if (!(nread == 1 || nread == 3)) error = TRUE;
  355.             if (nread == 1 || gmtdefs.color_model == 0) {    /* Read r,g,b */
  356.                 if (nread == 1)    /* Grayshades given */
  357.                     gmtdefs.nan_rgb[0] = gmtdefs.nan_rgb[1] = gmtdefs.nan_rgb[2] = r0;
  358.                 else {
  359.                     gmtdefs.nan_rgb[0] = r0;
  360.                     gmtdefs.nan_rgb[1] = g0;
  361.                     gmtdefs.nan_rgb[2] = b0;
  362.                 }
  363.             }
  364.             else
  365.                 gmt_hsv_to_rgb (&gmtdefs.nan_rgb[0], &gmtdefs.nan_rgb[1], &gmtdefs.nan_rgb[2], r0, g0, b0);
  366.             if (!(gmtdefs.nan_rgb[0] == gmtdefs.nan_rgb[1] && gmtdefs.nan_rgb[0] == gmtdefs.nan_rgb[2])) gmt_gray = FALSE;
  367.             if (gmt_gray && !(gmtdefs.nan_rgb[0] == 0 || gmtdefs.nan_rgb[0] == 255)) gmt_b_and_w = FALSE;
  368.             continue;
  369.         }
  370.  
  371.         nread = sscanf (line, "%lf %lf %lf %lf %lf %lf %lf %lf", 
  372.             &gmt_lut[n].z_low, &r0, &g0, &b0, &gmt_lut[n].z_high, &r1, &g1, &b1);
  373.         if (!(nread == 4 || nread == 8)) error = TRUE;
  374.             
  375.         if (nread == 4 || gmtdefs.color_model == 0) {    /* Read r,b,g or gray */
  376.             if (nread == 4)    {    /* Grayshades given */
  377.                 gmt_lut[n].z_high = g0;
  378.                 r1 = b0;
  379.                 gmt_lut[n].rgb_low[0] = gmt_lut[n].rgb_low[1] = gmt_lut[n].rgb_low[2] = r0;
  380.                 gmt_lut[n].rgb_high[0] = gmt_lut[n].rgb_high[1] = gmt_lut[n].rgb_high[2] = r1;
  381.             }
  382.             else {
  383.                 gmt_lut[n].rgb_low[0] = r0;
  384.                 gmt_lut[n].rgb_low[1] = g0;
  385.                 gmt_lut[n].rgb_low[2] = b0;
  386.                 gmt_lut[n].rgb_high[0] = r1;
  387.                 gmt_lut[n].rgb_high[1] = g1;
  388.                 gmt_lut[n].rgb_high[2] = b1;
  389.             }
  390.         }
  391.         else {    /* Read hue, saturation, value */
  392.             gmt_hsv_to_rgb (&gmt_lut[n].rgb_low[0], &gmt_lut[n].rgb_low[1], &gmt_lut[n].rgb_low[2], r0, g0, b0);
  393.             gmt_hsv_to_rgb (&gmt_lut[n].rgb_high[0], &gmt_lut[n].rgb_high[1], &gmt_lut[n].rgb_high[2], r1, g1, b1);
  394.         }
  395.         
  396.         if (gmt_lut[n].rgb_low[0] != -1) {    /* Ok to check these colors */
  397.             gmt_lut[n].skip = FALSE;
  398.             if (!(gmt_lut[n].rgb_low[0] == gmt_lut[n].rgb_low[1] && 
  399.                 gmt_lut[n].rgb_low[0] == gmt_lut[n].rgb_low[2])) gmt_gray = FALSE;
  400.             if (!(gmt_lut[n].rgb_high[0] == gmt_lut[n].rgb_high[1] && 
  401.                 gmt_lut[n].rgb_high[0] == gmt_lut[n].rgb_high[2])) gmt_gray = FALSE;
  402.             if (gmt_gray && !(gmt_lut[n].rgb_low[0] == 0 || gmt_lut[n].rgb_low[0] == 255)) gmt_b_and_w = FALSE;
  403.             if (gmt_gray && !(gmt_lut[n].rgb_high[0] == 0 || gmt_lut[n].rgb_high[0] == 255)) gmt_b_and_w = FALSE;
  404.             for (i = 0; !gmt_continuous && i < 3; i++)
  405.                 if (gmt_lut[n].rgb_low[i] != gmt_lut[n].rgb_high[i]) gmt_continuous = TRUE;
  406.         }
  407.         else {
  408.             gmt_lut[n].skip = TRUE;    /* Dont paint this slice */
  409.             for (i = 0; i < 3; i++) gmt_lut[n].rgb_low[i] = gmt_lut[n].rgb_high[i] = gmtdefs.page_rgb[i];
  410.         }
  411.  
  412.         /* Determine if psscale need to label these steps */
  413.  
  414.         c = line[strlen(line)-2]; 
  415.         if (c == 'L' || c == 'l')
  416.             gmt_lut[n].anot = 1;
  417.         else if (c == 'U' || c == 'u')
  418.             gmt_lut[n].anot = 2;
  419.         else if (c == 'B' || c == 'b')
  420.             gmt_lut[n].anot = 3;
  421.         n++;
  422.         if (n == n_alloc) {
  423.             n_alloc += GMT_SMALL_CHUNK;
  424.             gmt_lut = (struct GMT_LUT *) memory ((char *)gmt_lut, n_alloc, sizeof (struct GMT_LUT), "read_cpt");
  425.         }
  426.     }
  427.     
  428.     if (error) {
  429.         fprintf (stderr, "GMT Fatal Error: Error when decoding %s - aborts!\n", cpt_file);
  430.         exit (-1);
  431.     }
  432.     gmt_lut = (struct GMT_LUT *) memory ((char *)gmt_lut, n, sizeof (struct GMT_LUT), "read_cpt");
  433.     gmt_n_colors = n;
  434.     for (i = 0, gap = FALSE; !gap && i < gmt_n_colors - 1; i++) if (gmt_lut[i].z_high != gmt_lut[i+1].z_low) gap = TRUE;
  435.     if (gap) {
  436.         fprintf (stderr, "GMT Fatal Error: Color palette table %s has gaps - aborts!\n", cpt_file);
  437.         exit (-1);
  438.     }
  439.     if (!gmt_gray) gmt_b_and_w = FALSE;
  440. }
  441.  
  442.  
  443. int get_index (value)
  444. double value; {
  445.     int index;
  446.     
  447.     if (bad_float (value))    /* Set to NaN color */
  448.         return (-1);
  449.     if (value < gmt_lut[0].z_low)    /* Set to background color */
  450.         return (-2);
  451.     if (value > gmt_lut[gmt_n_colors-1].z_high)    /* Set to foreground color */
  452.         return (-3);
  453.     
  454.     /* Must search for correct index */
  455.     
  456.     index = 0;
  457.     while (index < gmt_n_colors && ! (value >= gmt_lut[index].z_low && value < gmt_lut[index].z_high) ) index++;
  458.     if (index == gmt_n_colors) index--;    /* Because we use <= for last range */
  459.     return (index);
  460. }
  461.  
  462. int get_rgb24 (value, r, g, b)
  463. double value;
  464. int *r, *g, *b; {
  465.     int index;
  466.     double rel;
  467.     
  468.     index = get_index (value);
  469.     
  470.     if (index == -1) {
  471.         *r = gmtdefs.nan_rgb[0];
  472.         *g = gmtdefs.nan_rgb[1];
  473.         *b = gmtdefs.nan_rgb[2];
  474.     }
  475.     else if (index == -2) {
  476.         *r = gmtdefs.background_rgb[0];
  477.         *g = gmtdefs.background_rgb[1];
  478.         *b = gmtdefs.background_rgb[2];
  479.     }
  480.     else if (index == -3) {
  481.         *r = gmtdefs.foreground_rgb[0];
  482.         *g = gmtdefs.foreground_rgb[1];
  483.         *b = gmtdefs.foreground_rgb[2];
  484.     }
  485.     else if (gmt_lut[index].skip) {        /* Set to page color for now */
  486.         *r = gmtdefs.page_rgb[0];
  487.         *g = gmtdefs.page_rgb[1];
  488.         *b = gmtdefs.page_rgb[2];
  489.     }
  490.     else {
  491.         rel = (value - gmt_lut[index].z_low) / (gmt_lut[index].z_high - gmt_lut[index].z_low);
  492.         *r = gmt_lut[index].rgb_low[0] + rint (rel * (gmt_lut[index].rgb_high[0] - gmt_lut[index].rgb_low[0]));
  493.         *g = gmt_lut[index].rgb_low[1] + rint (rel * (gmt_lut[index].rgb_high[1] - gmt_lut[index].rgb_low[1]));
  494.         *b = gmt_lut[index].rgb_low[2] + rint (rel * (gmt_lut[index].rgb_high[2] - gmt_lut[index].rgb_low[2]));
  495.     }
  496.     return (index);
  497. }
  498.  
  499. int gmt_rgb_to_hsv (r, g, b, h, s, v)
  500. int r, g, b;
  501. double *h, *s, *v; {
  502.     double xr, xg, xb, r_dist, g_dist, b_dist, max_v, min_v, diff, idiff;
  503.     
  504.     xr = r * I_255;
  505.     xg = g * I_255;
  506.     xb = b * I_255;
  507.     max_v = MAX (MAX (xr, xg), xb);
  508.     min_v = MIN (MIN (xr, xg), xb);
  509.     diff = max_v - min_v;
  510.     *v = max_v;
  511.     *s = (max_v == 0.0) ? 0.0 : diff / max_v;
  512.     *h = 0.0;
  513.     if ((*s) == 0.0) return;    /* Hue is undefined */
  514.     idiff = 1.0 / diff;
  515.     r_dist = (max_v - xr) * idiff;
  516.     g_dist = (max_v - xg) * idiff;
  517.     b_dist = (max_v - xb) * idiff;
  518.     if (xr == max_v)
  519.         *h = b_dist - g_dist;
  520.     else if (xg == max_v)
  521.         *h = 2.0 + r_dist - b_dist;
  522.     else
  523.         *h = 4.0 + g_dist - r_dist;
  524.     (*h) *= 60.0;
  525.     if ((*h) < 0.0) (*h) += 360.0;
  526. }
  527.     
  528. int gmt_hsv_to_rgb (r, g, b, h, s, v)
  529. int *r, *g, *b;
  530. double h, s, v; {
  531.     int i;
  532.     double f, p, q, t, rr, gg, bb;
  533.     
  534.     if (s == 0.0)
  535.         *r = *g = *b = floor (255.999 * v);
  536.     else {
  537.         while (h >= 360.0) h -= 360.0;
  538.         h /= 60.0;
  539.         i = h;
  540.         f = h - i;
  541.         p = v * (1.0 - s);
  542.         q = v * (1.0 - (s * f));
  543.         t = v * (1.0 - (s * (1.0 - f)));
  544.         switch (i) {
  545.             case 0:
  546.                 rr = v;    gg = t;    bb = p;
  547.                 break;
  548.             case 1:
  549.                 rr = q;    gg = v;    bb = p;
  550.                 break;
  551.             case 2:
  552.                 rr = p;    gg = v;    bb = t;
  553.                 break;
  554.             case 3:
  555.                 rr = p;    gg = q;    bb = v;
  556.                 break;
  557.             case 4:
  558.                 rr = t;    gg = p;    bb = v;
  559.                 break;
  560.             case 5:
  561.                 rr = v;    gg = p;    bb = q;
  562.                 break;
  563.         }
  564.         
  565.         *r = floor (rr * 255.999);
  566.         *g = floor (gg * 255.999);
  567.         *b = floor (bb * 255.999);
  568.     }
  569. }
  570.  
  571. int illuminate (intensity, r, g, b)
  572. double intensity;
  573. int *r, *g, *b; {
  574.     double h, s, v;
  575.     
  576.     if (bad_float (intensity)) return;
  577.     if (intensity == 0.0) return;
  578.     if (fabs (intensity) > 1.0) intensity = copysign (1.0, intensity);
  579.     
  580.     gmt_rgb_to_hsv (*r, *g, *b, &h, &s, &v);
  581.     if (intensity > 0) {
  582.         if (s != 0.0) s = (1.0 - intensity) * s + intensity * gmtdefs.hsv_max_saturation;
  583.         v = (1.0 - intensity) * v + intensity * gmtdefs.hsv_max_value;
  584.     }
  585.     else {
  586.         if (s != 0.0) s = (1.0 + intensity) * s - intensity * gmtdefs.hsv_min_saturation;
  587.         v = (1.0 + intensity) * v - intensity * gmtdefs.hsv_min_value;
  588.     }
  589.     if (v < 0.0) v = 0.0;
  590.     if (s < 0.0) s = 0.0;
  591.     if (v > 1.0) v = 1.0;
  592.     if (s > 1.0) s = 1.0;
  593.     gmt_hsv_to_rgb (r, g, b, h, s, v);
  594. }
  595.  
  596. /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  597.  * AKIMA computes the coefficients for a quasi-cubic hermite spline.
  598.  * Same algorithm as in the IMSL library.
  599.  * Programmer:    Paul Wessel
  600.  * Date:    16-JAN-1987
  601.  * Ver:        v.1-pc
  602.  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  603.  */
  604.  
  605. int akima (x,y,nx,c)
  606. int nx;
  607. double x[], y[], c[]; {
  608.     int i, no;
  609.     double t1, t2, b, rm1, rm2, rm3, rm4;
  610.  
  611.     if (nx < 4) return (-1);
  612.     for (i = 1; i < nx; i++) if (x[i] <= x[i-1]) return (-2);
  613.   
  614.     rm3 = (y[1] - y[0])/(x[1] - x[0]);
  615.     t1 = rm3 - (y[1] - y[2])/(x[1] - x[2]);
  616.     rm2 = rm3 + t1;
  617.     rm1 = rm2 + t1;
  618.     
  619.     /* get slopes */
  620.     
  621.     no = nx - 2;
  622.     for (i = 0; i < nx; i++) {
  623.         if (i >= no)
  624.             rm4 = rm3 - rm2 + rm3;
  625.         else
  626.             rm4 = (y[i+2] - y[i+1])/(x[i+2] - x[i+1]);
  627.         t1 = fabs(rm4 - rm3);
  628.         t2 = fabs(rm2 - rm1);
  629.         b = t1 + t2;
  630.         c[3*i] = (b != 0.0) ? (t1*rm2 + t2*rm3) / b : 0.5*(rm2 + rm3);
  631.         rm1 = rm2;
  632.         rm2 = rm3;
  633.         rm3 = rm4;
  634.     }
  635.     no = nx - 1;
  636.     
  637.     /* compute the coefficients for the nx-1 intervals */
  638.     
  639.     for (i = 0; i < no; i++) {
  640.         t1 = 1.0 / (x[i+1] - x[i]);
  641.         t2 = (y[i+1] - y[i])*t1;
  642.         b = (c[3*i] + c[3*i+3] - t2 - t2)*t1;
  643.         c[3*i+2] = b*t1;
  644.         c[3*i+1] = -b + (t2 - c[3*i])*t1;
  645.     }
  646.     return (0);
  647. }
  648.  
  649. /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  650.  * CSPLINE computes the coefficients for a natural cubic spline.
  651.  * To evaluate, call csplint
  652.  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  653.  */
  654.  
  655. int cspline (x, y, n, c)
  656. double x[], y[], c[];
  657. int n;
  658. {
  659.     int i, k;
  660.     double p, s, *u;
  661.  
  662.     u = (double *) memory (CNULL, n, sizeof (double), "cspline");
  663.     c[1] = c[n] = u[1] = 0.0;
  664.     for (i = 1; i < n-1; i++) {
  665.         s = (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
  666.         p = s * c[i-1] + 2.0;
  667.         c[i] = (s - 1.0) / p;
  668.         u[i] = (y[i+1] - y[i]) / (x[i+1] - x[i]) - (y[i] - y[i-1]) / (x[i] - x[i-1]);
  669.         u[i] = (6.0 * u[i] / (x[i+1] - x[i-1]) - s * u[i-1]) / p;
  670.     }
  671.     for (k = n-2; k >= 0; k--) c[k] = c[k] * c[k+1] + u[k];
  672.     free ((char *)u);
  673.     
  674.     return (0);
  675. }
  676.  
  677. double csplint (x, y, c, n, xp, klo)
  678. double x[], y[], c[], xp;
  679. int n, klo;
  680. {
  681.     int khi;
  682.     double h, ih, b, a, yp;
  683.  
  684.     khi = klo + 1;
  685.     h = x[khi] - x[klo];
  686.     ih = 1.0 / h;
  687.     a = (x[khi] - xp) * ih;
  688.     b = (xp - x[klo]) * ih;
  689.     yp = a * y[klo] + b * y[khi] + ((a*a*a - a) * c[klo] + (b*b*b - b) * c[khi]) * (h*h) / 6.0;
  690.     
  691.     return (yp);
  692. }
  693.  
  694. /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  695.  * INTPOL will interpolate from the dataset <x,y> onto a new set <u,v>
  696.  * where <x,y> and <u> is supplied by the user. <v> is returned. The 
  697.  * parameter mode governs what interpolation scheme that will be used.
  698.  * If u[i] is outside the range of x, then v[i] will contain NaN.
  699.  *
  700.  * input:  x = x-values of input data
  701.  *       y = y-values "    "     "
  702.  *       n = number of data pairs
  703.  *       m = number of desired interpolated values
  704.  *       u = x-values of these points
  705.  *    mode = type of interpolation
  706.  *         mode = 0 : Linear interpolation
  707.  *         mode = 1 : Quasi-cubic hermite spline (akima)
  708.  *         mode = 2 : Natural cubic spline (cubspl)
  709.  * output: v = y-values at interpolated points
  710.  *
  711.  * Programmer:    Paul Wessel
  712.  * Date:    16-JAN-1987
  713.  * Ver:        v.2
  714.  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  715.  */
  716.  
  717. int intpol (x,y,n,m,u,v,mode)
  718. double x[], y[], u[], v[];
  719. int n, m, mode; {
  720.     int i, j, err_flag = 0;
  721.     double dx, *c, csplint ();
  722.     
  723.     if (n < 4) mode = 0;
  724.     if (mode > 0) c = (double *) memory (CNULL, 3*n, sizeof(double), "intpol");
  725.  
  726.     if (mode == 0) {    /* Linear interpolation */
  727.         for (i = 1; i < n; i++) if (x[i] <= x[i-1]) return (-1);
  728.     }
  729.     else if (mode == 1)     /* Akimas spline */
  730.         err_flag = akima (x,y,n,c);
  731.     else if (mode == 2)    /* Natural cubic spline */
  732.         err_flag = cspline (x,y,n,c);
  733.     else
  734.         err_flag = -4;
  735.     if (err_flag != 0) {
  736.         free ((char *)c);
  737.         return (err_flag);
  738.     }
  739.     
  740.     /* Compute the interpolated values from the coefficients */
  741.     
  742.     j = 0;
  743.     for (i = 0; i < m; i++) {
  744.         if (u[i] < x[0] || u[i] > x[n-1]) {    /* Desired point outside data range */
  745.             v[i] = gmt_NaN;
  746.             continue;
  747.         } 
  748.         while (x[j] > u[i] && j > 0) j--;    /* In case u is not sorted */
  749.         while (j < n && x[j] <= u[i]) j++;
  750.         if (j == n) j--;
  751.         if (j > 0) j--;
  752.         switch (mode) {
  753.             case 0:
  754.                 dx = u[i] - x[j];
  755.                 v[i] = (y[j+1]-y[j])*dx/(x[j+1]-x[j]) + y[j];
  756.                 break;
  757.             case 1:
  758.                 dx = u[i] - x[j];
  759.                 v[i] = ((c[3*j+2]*dx + c[3*j+1])*dx + c[3*j])*dx + y[j];
  760.                 break;
  761.             case 2:
  762.                 v[i] = csplint (x, y, c, n, u[i], j);
  763.                 break;
  764.         }
  765.     }
  766.     if (mode > 0) free ((char *)c);
  767.     return (0);
  768. }
  769.  
  770. char *memory (prev_addr, n, size, progname)
  771. char *prev_addr, *progname;
  772. int n, size; {
  773.     char *tmp;
  774.  
  775.     if (n == 0) return((char *) NULL); /* Take care of n = 0 */
  776.     
  777.     if (prev_addr) {
  778.         if ((tmp = realloc ((char *) prev_addr, (unsigned) (n * size))) == CNULL) {
  779.             fprintf (stderr, "GMT Fatal Error: %s could not reallocate more memory, n = %d\n", progname, n);
  780.             exit (-1);
  781.         }
  782.     }
  783.     else {
  784.         if ((tmp = calloc ((unsigned) n, (unsigned) size)) == CNULL) {
  785.             fprintf (stderr, "GMT Fatal Error: %s could not allocate memory, n = %d\n", progname, n);
  786.             exit (-1);
  787.         }
  788.     }
  789.     return (tmp);
  790. }
  791.  
  792. #define CONTOUR_IJ(i,j) ((i) + (j) * nx)
  793.  
  794. int contours (grd, header, smooth_factor, int_scheme, side, edge, first, x_array, y_array)
  795. float grd[];
  796. struct GRD_HEADER *header;
  797. int smooth_factor, int_scheme, *side, first;
  798. int edge[];
  799. double *x_array[], *y_array[]; {
  800.     /* The routine finds the zero-contour in the grd dataset.  it assumes that
  801.      * no node has a value exactly == 0.0.  If more than max points are found
  802.      * trace_contour will try to allocate more memory in blocks of GMT_CHUNK points
  803.      */
  804.      
  805.     static int i0, j0;
  806.     int i, j, ij, n = 0, n2, n_edges, edge_word, edge_bit, nx, ny, nans = 0;
  807.     BOOLEAN go_on = TRUE;
  808.     double x0, y0, r, west, east, south, north, dx, dy, xinc2, yinc2, *x, *y, *x2, *y2;
  809.     int p[5], i_off[5], j_off[5], k_off[5], offset;
  810.     unsigned int bit[32];
  811.      
  812.     nx = header->nx;    ny = header->ny;
  813.     west = header->x_min;    east = header->x_max;
  814.     south = header->y_min;    north = header->y_max;
  815.     dx = header->x_inc;    dy = header->y_inc;
  816.     xinc2 = (header->node_offset) ? 0.5 * dx : 0.0;
  817.     yinc2 = (header->node_offset) ? 0.5 * dy : 0.0;
  818.     x = *x_array;    y = *y_array;
  819.     
  820.     n_edges = ny * (int) ceil (nx / 16.0);
  821.     offset = n_edges / 2;
  822.      
  823.     /* Reset edge-flags to zero, if necessary */
  824.     if (first) {    /* Set i0,j0 for southern boundary */
  825.         memset ((char *)edge, 0, n_edges * sizeof (int));
  826.         i0 = 0;
  827.         j0 = ny - 1;
  828.     }
  829.     p[0] = p[4] = 0;    p[1] = 1;    p[2] = 1 - nx;    p[3] = -nx;
  830.     i_off[0] = i_off[2] = i_off[3] = i_off[4] = 0;    i_off[1] =  1;
  831.     j_off[0] = j_off[1] = j_off[3] = j_off[4] = 0;    j_off[2] = -1;
  832.     k_off[0] = k_off[2] = k_off[4] = 0;    k_off[1] = k_off[3] = 1;
  833.     for (i = 1, bit[0] = 1; i < 32; i++) bit[i] = bit[i-1] << 1;
  834.  
  835.     switch (*side) {
  836.         case 0:    /* Southern boundary */
  837.             for (i = i0, j = j0, ij = (ny - 1) * nx + i0; go_on && i < nx-1; i++, ij++) {
  838.                 edge_word = ij / 32;
  839.                 edge_bit = ij % 32;
  840.                 if (start_trace (grd[ij+1], grd[ij], edge, edge_word, edge_bit, bit)) { /* Start tracing countour */
  841.                     r = grd[ij+1] - grd[ij];
  842.                     x0 = west + (i - grd[ij]/r)*dx + xinc2;
  843.                     y0 = south + yinc2;
  844.                     edge[edge_word] |= bit[edge_bit];
  845.                     n = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x, &y, i, j, 0, offset, i_off, j_off, k_off, p, bit, &nans);
  846.                     n = smooth_contour (&x, &y, n, smooth_factor, int_scheme);
  847.                     go_on = FALSE;
  848.                     i0 = i + 1;    j0 = j;
  849.                 }
  850.             }
  851.             if (n == 0) {
  852.                 i0 = nx - 2;
  853.                 j0 = ny - 1;
  854.                 (*side)++;
  855.             }
  856.             break;
  857.             
  858.         case 1:        /* Eastern boundary */
  859.         
  860.             for (j = j0, ij = nx * (j0 + 1) - 1, i = i0; go_on && j > 0; j--, ij -= nx) {
  861.                 edge_word = ij / 32 + offset;
  862.                 edge_bit = ij % 32;
  863.                 if (start_trace (grd[ij-nx], grd[ij], edge, edge_word, edge_bit, bit)) { /* Start tracing countour */
  864.                     r = grd[ij] - grd[ij-nx];
  865.                     x0 = east - xinc2;
  866.                     y0 = north - (j - grd[ij]/r) * dy - yinc2;
  867.                     edge[edge_word] |= bit[edge_bit];
  868.                     n = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x, &y, i, j, 1, offset, i_off, j_off, k_off, p, bit, &nans);
  869.                     n = smooth_contour (&x, &y, n, smooth_factor, int_scheme);
  870.                     go_on = FALSE;
  871.                     i0 = i;    j0 = j - 1;
  872.                 }
  873.             }
  874.             if (n == 0) {
  875.                 i0 = nx - 2;
  876.                 j0 = 1;
  877.                 (*side)++;
  878.             }
  879.             break;
  880.  
  881.         case 2:        /* Northern boundary */    
  882.         
  883.             for (i = i0, j = j0, ij = i0; go_on && i >= 0; i--, ij--) {
  884.                 edge_word = ij / 32;
  885.                 edge_bit = ij % 32;
  886.                 if (start_trace (grd[ij], grd[ij+1], edge, edge_word, edge_bit, bit)) { /* Start tracing countour */
  887.                     r = grd[ij+1] - grd[ij];
  888.                     x0 = west + (i - grd[ij]/r)*dx + xinc2;
  889.                     y0 = north - yinc2;
  890.                     edge[edge_word] |= bit[edge_bit];
  891.                     n = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x, &y, i, j, 2, offset, i_off, j_off, k_off, p, bit, &nans);
  892.                     n = smooth_contour (&x, &y, n, smooth_factor, int_scheme);
  893.                     go_on = FALSE;
  894.                     i0 = i - 1;    j0 = j;
  895.                 }
  896.             }
  897.             if (n == 0) {
  898.                 i0 = 0;
  899.                 j0 = 1;
  900.                 (*side)++;
  901.             }
  902.             break;
  903.         
  904.         case 3:        /* Western boundary */
  905.         
  906.             for (j = j0, ij = j * nx, i = i0; go_on && j < ny; j++, ij += nx) {
  907.                 edge_word = ij / 32 + offset;
  908.                 edge_bit = ij % 32;
  909.                 if (start_trace (grd[ij], grd[ij-nx], edge, edge_word, edge_bit, bit)) { /* Start tracing countour */
  910.                     r = grd[ij] - grd[ij-nx];
  911.                     x0 = west + xinc2;
  912.                     y0 = north - (j - grd[ij] / r) * dy - yinc2;
  913.                     edge[edge_word] |= bit[edge_bit];
  914.                     n = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x, &y, i, j, 3, offset, i_off, j_off, k_off, p, bit, &nans);
  915.                     n = smooth_contour (&x, &y, n, smooth_factor, int_scheme);
  916.                     go_on = FALSE;
  917.                     i0 = i;    j0 = j + 1;
  918.                 }
  919.             }
  920.             if (n == 0) {
  921.                 i0 = 1;
  922.                 j0 = 1;
  923.                 (*side)++;
  924.             }
  925.             break;
  926.  
  927.         case 4:    /* Then loop over interior boxes */
  928.         
  929.             for (j = j0; go_on && j < ny; j++) {
  930.                 ij = i0 + j * nx;
  931.                 for (i = i0; go_on && i < nx-1; i++, ij++) {
  932.                     edge_word = ij / 32 + offset;
  933.                     edge_bit = ij % 32;
  934.                     if (start_trace (grd[ij], grd[ij-nx], edge, edge_word, edge_bit, bit)) { /* Start tracing countour */
  935.                         r = grd[ij] - grd[ij-nx];
  936.                         x0 = west + i*dx + xinc2;
  937.                         y0 = north - (j - grd[ij] / r) * dy - yinc2;
  938.                         /* edge[edge_word] |= bit[edge_bit]; */
  939.                         nans = 0;
  940.                         n = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x, &y, i, j, 3, offset, i_off, j_off, k_off, p, bit, &nans);
  941.                         if (nans) {    /* Must trace in other direction, then splice */
  942.                             n2 = trace_contour (grd, header, x0, y0, edge, smooth_factor, &x2, &y2, i-1, j, 1, offset, i_off, j_off, k_off, p, bit, &nans);
  943.                             n = splice_contour (&x, &y, n, &x2, &y2, n2);
  944.                             free ((char *)x2);
  945.                             free ((char *)y2);
  946.                             
  947.                             n = smooth_contour (&x, &y, n, smooth_factor, int_scheme);
  948.                         }
  949.                         i0 = i + 1;
  950.                         go_on = FALSE;
  951.                         j0 = j;
  952.                     }
  953.                 }
  954.                 if (go_on) i0 = 1;
  955.             }
  956.             if (n == 0) (*side)++;
  957.             break;
  958.         default:
  959.             break;
  960.     }
  961.     *x_array = x;    *y_array = y;
  962.     return (n);
  963. }
  964.  
  965. int start_trace (first, second, edge, edge_word, edge_bit, bit)
  966. float first, second;
  967. int edge[], edge_word,edge_bit;
  968. unsigned int bit[]; {
  969.     /* First make sure we havent been here before */
  970.  
  971.     if ( (edge[edge_word] & bit[edge_bit]) ) return (FALSE);
  972.  
  973.     /* Then make sure both points are real numbers */
  974.  
  975.     if (bad_float ((double)first) ) return (FALSE);
  976.     if (bad_float ((double)second) ) return (FALSE);
  977.  
  978.     /* Finally return TRUE if values of opposite sign */
  979.  
  980.     return ( first * second < 0.0);
  981. }
  982.  
  983. int splice_contour (x, y, n, x2, y2, n2)
  984. double *x[], *y[], *x2[], *y2[];
  985. int n, n2; {    /* Basically does a "tail -r" on the array x,y and append arrays x2/y2 */
  986.  
  987.     int i, j, m;
  988.     double *xtmp, *ytmp, *xa, *ya, *xb, *yb;
  989.     
  990.     xa = *x;    ya = *y;
  991.     xb = *x2;    yb = *y2;
  992.     
  993.     m = n + n2 - 1;    /* Total length since one point is shared */
  994.     
  995.     /* First copy and reverse first contour piece */
  996.     
  997.     xtmp = (double *) memory (CNULL, n , sizeof (double), "splice_contour");
  998.     ytmp = (double *) memory (CNULL, n, sizeof (double), "splice_contour");
  999.     
  1000.     memcpy ((char *)xtmp, (char *)xa, n * sizeof (double));
  1001.     memcpy ((char *)ytmp, (char *)ya, n * sizeof (double));
  1002.     
  1003.     xa = (double *) memory ((char *)xa, m, sizeof (double), "splice_contour");
  1004.     ya = (double *) memory ((char *)ya, m, sizeof (double), "splice_contour");
  1005.     
  1006.     for (i = 0; i < n; i++) xa[i] = xtmp[n-i-1];
  1007.     for (i = 0; i < n; i++) ya[i] = ytmp[n-i-1];
  1008.     
  1009.     /* Then append second piece */
  1010.     
  1011.     for (j = 1, i = n; j < n2; j++, i++) xa[i] = xb[j];
  1012.     for (j = 1, i = n; j < n2; j++, i++) ya[i] = yb[j];
  1013.  
  1014.     free ((char *)xtmp);
  1015.     free ((char *)ytmp);
  1016.     
  1017.     *x = xa;
  1018.     *y = ya;
  1019.     
  1020.     return (m);
  1021. }
  1022.     
  1023. int trace_contour (grd, header, x0, y0, edge, smooth_factor, x_array, y_array, i, j, kk, offset, i_off, j_off, k_off, p, bit, nan_flag)
  1024. float grd[];
  1025. struct GRD_HEADER *header;
  1026. int edge[], smooth_factor, i, j, kk, offset, i_off[], j_off[], k_off[], p[];
  1027. unsigned int bit[];
  1028. double x0, y0, *x_array[], *y_array[];
  1029. int *nan_flag; {
  1030.     int side = 0, n = 1, k, k0, ij, ij0, n_cuts, k_index[2], kk_opposite, more;
  1031.     int edge_word, edge_bit, m, n_nan, nx, ny, n_alloc, ij_in;
  1032.     double xk[4], yk[4], r, dr[2];
  1033.     double west, north, dx, dy, xinc2, yinc2, *xx, *yy;
  1034.     
  1035.     west = header->x_min;    north = header->y_max;
  1036.     dx = header->x_inc;    dy = header->y_inc;
  1037.     nx = header->nx;    ny = header->ny;
  1038.     xinc2 = (header->node_offset) ? 0.5 * dx : 0.0;
  1039.     yinc2 = (header->node_offset) ? 0.5 * dy : 0.0;
  1040.     
  1041.     n_alloc = GMT_CHUNK;
  1042.     m = n_alloc - 2;
  1043.  
  1044.     xx = (double *) memory (CNULL, n_alloc, sizeof (double), "trace_contour");
  1045.     yy = (double *) memory (CNULL, n_alloc, sizeof (double), "trace_contour");
  1046.     
  1047.     xx[0] = x0;    yy[0] = y0;
  1048.     ij_in = j * nx + i - 1;
  1049.     
  1050.     more = TRUE;
  1051.     do {
  1052.         ij = i + j * nx;
  1053.         x0 = west + i * dx + xinc2;
  1054.         y0 = north - j * dy - yinc2;
  1055.         n_cuts = 0;
  1056.         k0 = kk;
  1057.  
  1058.         for (k = n_nan = 0; k < 4; k++) {    /* Loop over box sides */
  1059.  
  1060.             /* Skip where we already have a cut (k == k0) */
  1061.             
  1062.             if (k == k0) continue;
  1063.             
  1064.             /* Skip if NAN encountered */
  1065.             
  1066.             if (bad_float ((double)grd[ij+p[k+1]]) || bad_float ( (double)grd[ij+p[k]])) {
  1067.                 n_nan++;
  1068.                 continue;
  1069.             }
  1070.             
  1071.             /* Skip edge already has been used */
  1072.             
  1073.             ij0 = CONTOUR_IJ (i + i_off[k], j + j_off[k]);
  1074.             edge_word = ij0 / 32 + k_off[k] * offset;
  1075.             edge_bit = ij0 % 32;
  1076.             if (edge[edge_word] & bit[edge_bit]) continue;
  1077.             
  1078.             /* Skip if no zero-crossing on this edge */
  1079.             
  1080.             if (grd[ij+p[k+1]] * grd[ij+p[k]] > 0.0) continue;
  1081.             
  1082.             /* Here we have a crossing */
  1083.             
  1084.             r = grd[ij+p[k+1]] - grd[ij+p[k]];
  1085.             
  1086.             if (k%2) {    /* Cutting a S-N line */
  1087.                 if (k == 1) {
  1088.                     xk[1] = x0 + dx;
  1089.                     yk[1] = y0 - grd[ij+p[k]]*dy/r;
  1090.                 }
  1091.                 else {
  1092.                     xk[3] = x0;
  1093.                     yk[3] = y0 + (1.+grd[ij+p[k]]/r)*dy;
  1094.                 }
  1095.             }
  1096.             else {    /* Cutting a E-W line */
  1097.                 if (k == 0) {
  1098.                     xk[0] = x0 - grd[ij+p[k]]*dx/r;
  1099.                     yk[0] = y0;
  1100.                 }
  1101.                 else {
  1102.                     xk[2] = x0 + (1.+grd[ij+p[k]]/r)*dx;
  1103.                     yk[2] = y0 + dy;
  1104.                 }
  1105.             }
  1106.             kk = k;
  1107.             n_cuts++;
  1108.         }
  1109.         
  1110.         if (n > m) {    /* Must try to allocate more memory */
  1111.             n_alloc += GMT_CHUNK;
  1112.             m += GMT_CHUNK;
  1113.             xx = (double *) memory ((char *)xx, n_alloc, sizeof (double), "trace_contour");
  1114.             yy = (double *) memory ((char *)yy, n_alloc, sizeof (double), "trace_contour");
  1115.         }
  1116.         if (n_cuts == 0) {    /* Close interior contour and return */
  1117.         /*    if (n_nan == 0 && (n > 0 && fabs (xx[n-1] - xx[0]) <= dx && fabs (yy[n-1] - yy[0]) <= dy)) { */
  1118.             if (ij == ij_in) {    /* Close iterior contour */
  1119.                 xx[n] = xx[0];
  1120.                 yy[n] = yy[0];
  1121.                 n++;
  1122.             }
  1123.             more = FALSE;
  1124.             *nan_flag = n_nan;
  1125.         }
  1126.         else if (n_cuts == 1) {    /* Draw a line to this point and keep tracing */
  1127.             xx[n] = xk[kk];
  1128.             yy[n] = yk[kk];
  1129.             n++;
  1130.         }
  1131.         else {    /* Saddle point, we decide to connect to the point nearest previous point */
  1132.             kk_opposite = (k0 + 2) % 4;    /* But it cannot be on opposite side */
  1133.             for (k = side = 0; k < 4; k++) {
  1134.                 if (k == k0 || k == kk_opposite) continue;
  1135.                 dr[side] = (xx[n-1]-xk[k])*(xx[n-1]-xk[k]) + (yy[n-1]-yk[k])*(yy[n-1]-yk[k]);
  1136.                 k_index[side++] = k;
  1137.             }
  1138.             kk = (dr[0] < dr[1]) ? k_index[0] : k_index[1];
  1139.             xx[n] = xk[kk];
  1140.             yy[n] = yk[kk];
  1141.             n++;
  1142.         }
  1143.         if (more) {    /* Mark this edge as used */
  1144.             ij0 = CONTOUR_IJ (i + i_off[kk], j + j_off[kk]);
  1145.             edge_word = ij0 / 32 + k_off[kk] * offset;
  1146.             edge_bit = ij0 % 32;
  1147.             edge[edge_word] |= bit[edge_bit];    
  1148.         }
  1149.         
  1150.         if ((kk == 0 && j == ny - 1) || (kk == 1 && i == nx - 2) ||
  1151.             (kk == 2 && j == 1) || (kk == 3 && i == 0))    /* Going out of grid */
  1152.             more = FALSE;
  1153.             
  1154.         /* Get next box (i,j,kk) */
  1155.         
  1156.         i -= (kk-2)%2;
  1157.         j -= (kk-1)%2;
  1158.         kk = (kk+2)%4;
  1159.         
  1160.     } while (more);
  1161.     
  1162.     xx = (double *) memory ((char *)xx, n, sizeof (double), "trace_contour");
  1163.     yy = (double *) memory ((char *)yy, n, sizeof (double), "trace_contour");
  1164.     
  1165.     *x_array = xx;    *y_array = yy;
  1166.     return (n);
  1167. }
  1168.  
  1169. int smooth_contour (x_in, y_in, n, sfactor, stype)
  1170. double *x_in[], *y_in[];    /* Input (x,y) points */
  1171. int n;        /* Number of input points */
  1172. int sfactor;    /* n_out = sfactor * n -1 */
  1173. int stype; {    /* Interpolation scheme used (0 = linear, 1 = Akima, 2 = Cubic spline */
  1174.     int i, j, k, n_out;
  1175.     double ds, t_next, *x, *y;
  1176.     double *t_in, *t_out, *x_tmp, *y_tmp, x0, x1, y0, y1;
  1177.     char *flag;
  1178.     
  1179.     if (sfactor == 0 || n < 4) return (n);    /* Need at least 4 points to call Akima */
  1180.     
  1181.     x = *x_in;    y = *y_in;
  1182.     
  1183.     n_out = sfactor * n - 1;    /* Number of new points */
  1184.     
  1185.     t_in = (double *) memory (CNULL, n, sizeof (double), "smooth_contour");
  1186.     t_out = (double *) memory (CNULL, n_out + n, sizeof (double), "smooth_contour");
  1187.     x_tmp = (double *) memory (CNULL, n_out + n, sizeof (double), "smooth_contour");
  1188.     y_tmp = (double *) memory (CNULL, n_out + n, sizeof (double), "smooth_contour");
  1189.     flag = memory (CNULL, n_out + n, sizeof (char), "smooth_contour");
  1190.     
  1191.     /* Create dummy distance values for input points */
  1192.     
  1193.     t_in[0] = 0.0;
  1194.     for (i = 1; i < n; i++)    t_in[i] = t_in[i-1] + hypot ((x[i]-x[i-1]), (y[i]-y[i-1]));
  1195.     
  1196.     /* Create equidistance output points */
  1197.     
  1198.     ds = t_in[n-1] / (n_out-1);
  1199.     t_next = ds;
  1200.     t_out[0] = 0.0;
  1201.     flag[0] = TRUE;
  1202.     for (i = j = 1; i < n_out; i++) {
  1203.         if (j < n && t_in[j] < t_next) {
  1204.             t_out[i] = t_in[j];
  1205.             flag[i] = TRUE;
  1206.             j++;
  1207.             n_out++;
  1208.         }
  1209.         else {
  1210.             t_out[i] = t_next;
  1211.             t_next += ds;
  1212.         }
  1213.     }
  1214.     t_out[n_out-1] = t_in[n-1];
  1215.     if (t_out[n_out-1] == t_out[n_out-2]) n_out--;
  1216.     flag[n_out-1] = TRUE;
  1217.     
  1218.     intpol (t_in, x, n, n_out, t_out, x_tmp, stype);
  1219.     intpol (t_in, y, n, n_out, t_out, y_tmp, stype);
  1220.     
  1221.     /* Make sure interpolated function is bounded on each segment interval */
  1222.     
  1223.     i = 0;
  1224.     while (i < (n_out - 1)) {
  1225.         j = i + 1;
  1226.         while (j < n_out && !flag[j]) j++;
  1227.         x0 = x_tmp[i];    x1 = x_tmp[j];
  1228.         if (x0 > x1) d_swap (x0, x1);
  1229.         y0 = y_tmp[i];    y1 = y_tmp[j];
  1230.         if (y0 > y1) d_swap (y0, y1);
  1231.         for (k = i + 1; k < j; k++) {
  1232.             if (x_tmp[k] < x0)
  1233.                 x_tmp[k] = x0 + 1.0e-10;
  1234.             else if (x_tmp[k] > x1)
  1235.                 x_tmp[k] = x1 - 1.0e-10;
  1236.             if (y_tmp[k] < y0)
  1237.                 y_tmp[k] = y0 + 1.0e-10;
  1238.             else if (y_tmp[k] > y1)
  1239.                 y_tmp[k] = y1 - 1.0e-10;
  1240.         }
  1241.         i = j;
  1242.     }
  1243.     
  1244.     free ((char *)x);
  1245.     free ((char *)y);
  1246.     
  1247.     *x_in = x_tmp;
  1248.     *y_in = y_tmp;
  1249.     
  1250.     free ((char *) t_in);
  1251.     free ((char *) t_out);
  1252.     free (flag);
  1253.     
  1254.     return (n_out);
  1255. }
  1256.  
  1257. int get_plot_array () {      /* Allocate more space for plot arrays */
  1258.  
  1259.     gmt_n_alloc += GMT_CHUNK;
  1260.     gmt_x_plot = (double *) memory ((char *)gmt_x_plot, gmt_n_alloc, sizeof (double), "gmt");
  1261.     gmt_y_plot = (double *) memory ((char *)gmt_y_plot, gmt_n_alloc, sizeof (double), "gmt");
  1262.     gmt_pen = (int *) memory ((char *)gmt_pen, gmt_n_alloc, sizeof (int), "gmt");
  1263. }
  1264.  
  1265. int free_plot_array () {
  1266.     if (gmt_n_alloc) {
  1267.         free ((char *)gmt_x_plot);
  1268.         free ((char *)gmt_y_plot);
  1269.         free ((char *)gmt_pen);
  1270.     }
  1271. }
  1272.  
  1273. int comp_double_asc (point_1, point_2)
  1274. double    *point_1, *point_2;
  1275. {
  1276.     /* Returns -1 if point_1 is < that point_2,
  1277.        +1 if point_2 > point_1, and 0 if they are equal
  1278.     */
  1279.     int bad_1, bad_2;
  1280.     
  1281.     bad_1 = bad_float ((*point_1));
  1282.     bad_2 = bad_float ((*point_2));
  1283.     
  1284.     if (bad_1 && bad_2) return (0);
  1285.     if (bad_1) return (1);
  1286.     if (bad_2) return (-1);
  1287.     
  1288.     if ( (*point_1) < (*point_2) )
  1289.         return (-1);
  1290.     else if ( (*point_1) > (*point_2) )
  1291.         return (1);
  1292.     else
  1293.         return (0);
  1294. }
  1295.  
  1296. int comp_float_asc (point_1, point_2)
  1297. float    *point_1, *point_2;
  1298. {
  1299.     /* Returns -1 if point_1 is < that point_2,
  1300.        +1 if point_2 > point_1, and 0 if they are equal
  1301.     */
  1302.     int bad_1, bad_2;
  1303.     
  1304.     bad_1 = bad_float ((double)(*point_1));
  1305.     bad_2 = bad_float ((double)(*point_2));
  1306.     
  1307.     if (bad_1 && bad_2) return (0);
  1308.     if (bad_1) return (1);
  1309.     if (bad_2) return (-1);
  1310.     
  1311.     if ( (*point_1) < (*point_2) )
  1312.         return (-1);
  1313.     else if ( (*point_1) > (*point_2) )
  1314.         return (1);
  1315.     else
  1316.         return (0);
  1317. }
  1318.  
  1319. int comp_int_asc (point_1, point_2)
  1320. int    *point_1, *point_2;
  1321. {
  1322.     /* Returns -1 if point_1 is < that point_2,
  1323.        +1 if point_2 > point_1, and 0 if they are equal
  1324.     */
  1325.     
  1326.     if ( (*point_1) < (*point_2) )
  1327.         return (-1);
  1328.     else if ( (*point_1) > (*point_2) )
  1329.         return (1);
  1330.     else
  1331.         return (0);
  1332. }
  1333.  
  1334. struct EPS *gmt_epsinfo (program)
  1335. char *program;
  1336. {
  1337.     /* Supply info about the EPS file that will be created */
  1338.      
  1339.     int fno[5], id, i, n, n_fonts, last, move_up = FALSE;
  1340.     
  1341.     double y, dy;
  1342.     
  1343.     char title[512];
  1344.      
  1345.     struct passwd *pw;
  1346.     struct EPS *new;
  1347. #ifdef NeXT
  1348.     uid_t getuid();
  1349. #endif
  1350.     new = (struct EPS *) memory (CNULL, 1, sizeof (struct EPS), "gmt_epsinfo");
  1351.      
  1352.      /* First estimate the boundingbox coordinates */
  1353.      
  1354.     new->x0 = new->y0 = 0;
  1355.     new->x1 = rint (gmt_ppu[gmtdefs.measure_unit] * (z_project.xmax - z_project.xmin + 2.0 * ((gmtdefs.overlay) ? 1.0 : gmtdefs.x_origin)));
  1356.     y = z_project.ymax - z_project.ymin + 2.0 * ((gmtdefs.overlay) ? 1.0 : gmtdefs.y_origin);
  1357.      
  1358.     if (frame_info.header[0]) {    /* Make space for header text */
  1359.         move_up = (MAPPING || frame_info.side[2] == 2);
  1360.         dy = ((gmtdefs.tick_length > 0.0) ? gmtdefs.tick_length : 0.0) + 2.5 * gmtdefs.anot_offset;
  1361.         dy += ((move_up) ? (gmtdefs.anot_font_size + gmtdefs.label_font_size) / gmt_ppu[gmtdefs.measure_unit] : 0.0) + 2.5 * gmtdefs.anot_offset;
  1362.         y += dy;
  1363.     }
  1364.     new->y1 = rint (gmt_ppu[gmtdefs.measure_unit] * y);
  1365.     
  1366.     /* Get font names used */
  1367.     
  1368.     id = 0;
  1369.     if (gmtdefs.unix_time) {
  1370.         fno[0] = 0;
  1371.         fno[1] = 1;
  1372.         id = 2;
  1373.     }
  1374.     
  1375.     if (frame_info.header[0]) fno[id++] = gmtdefs.header_font;
  1376.     
  1377.     if (frame_info.label[0][0] || frame_info.label[1][0] || frame_info.label[2][0]) fno[id++] = gmtdefs.label_font;
  1378.     
  1379.     fno[id++] = gmtdefs.anot_font;
  1380.     
  1381.     qsort ((char *)fno, id, sizeof (int), comp_int_asc);
  1382.     
  1383.     last = -1;
  1384.     for (i = n_fonts = 0; i < id; i++) {
  1385.         if (fno[i] != last) {
  1386.             new->font[n_fonts++] = font_name[fno[i]];
  1387.             last = fno[i];
  1388.         }
  1389.     }
  1390.     new->font[n_fonts] = CNULL;
  1391.     
  1392.     /* Get user name and date */
  1393. #ifdef NeXT
  1394.     pw = getpwuid ( ( (int)getuid () ) );
  1395. #elif sony
  1396.     pw = getpwuid ( ( (int)getuid () ) );
  1397. #else
  1398.     pw = getpwnam ( cuserid ( CNULL));
  1399. #endif
  1400.  
  1401.     n = strlen (pw->pw_gecos) + 1;
  1402.     new->name = memory (CNULL, n, 1, "gmt_epsinfo");
  1403.     strcpy (new->name, pw->pw_gecos);
  1404.     
  1405.     sprintf (title, "GMT v%s Document from %s\0", GMT_VERSION, program);
  1406.     new->title = memory (CNULL, strlen (title) + 1, 1, "gmt_epsinfo");
  1407.     strcpy (new->title, title);
  1408.     
  1409.     return (new);
  1410. }
  1411.  
  1412. int median(x, n, xmin, xmax, m_initial, med)
  1413.  
  1414. double    *x, xmin, xmax, m_initial, *med;
  1415. int    n;
  1416.  
  1417. {
  1418.     double    lower_bound, upper_bound, m_guess, t_0, t_1, t_middle;
  1419.     double    lub, glb, xx, temp;
  1420.     int    i, n_above, n_below, n_equal, n_lub, n_glb;
  1421.     int    finished = FALSE;
  1422.     int    iteration = 0;
  1423.     
  1424.     m_guess = m_initial;
  1425.     lower_bound = xmin;
  1426.     upper_bound = xmax;
  1427.     t_0 = 0;
  1428.     t_1 = n - 1;
  1429.     t_middle = 0.5 * (n - 1);
  1430.     
  1431.     do {
  1432.  
  1433.         n_above = n_below = n_equal = n_lub = n_glb = 0;
  1434.         lub = xmax;
  1435.         glb = xmin;
  1436.         
  1437.         for (i = 0; i < n; i++) {
  1438.         
  1439.             xx = x[i];
  1440.             if (xx == m_guess) {
  1441.                 n_equal++;
  1442.             }
  1443.             else if (xx > m_guess) {
  1444.                 n_above++;
  1445.                 if (xx < lub) {
  1446.                     lub = xx;
  1447.                     n_lub = 1;
  1448.                 }
  1449.                 else if (xx == lub) {
  1450.                     n_lub++;
  1451.                 }
  1452.             }
  1453.             else {
  1454.                 n_below++;
  1455.                 if (xx > glb) {
  1456.                     glb = xx;
  1457.                     n_glb = 1;
  1458.                 }
  1459.                 else if (xx == glb) {
  1460.                     n_glb++;
  1461.                 }
  1462.             }
  1463.         }
  1464.         
  1465.         iteration++;
  1466.  
  1467.         /* Now test counts, watch multiple roots, think even/odd:  */
  1468.         
  1469.         if ( (abs(n_above - n_below)) <= n_equal) {
  1470.             
  1471.             if (n_equal) {
  1472.                 *med = m_guess;
  1473.             }
  1474.             else {
  1475.                 *med = 0.5 * (lub + glb);
  1476.             }
  1477.             finished = TRUE;
  1478.         }
  1479.         else if ( (abs( (n_above - n_lub) - (n_below + n_equal) ) ) < n_lub) {
  1480.         
  1481.             *med = lub;
  1482.             finished = TRUE;
  1483.         }
  1484.         else if ( (abs( (n_below - n_glb) - (n_above + n_equal) ) ) < n_glb) {
  1485.         
  1486.             *med = glb;
  1487.             finished = TRUE;
  1488.         }
  1489.         /* Those cases found the median; the next two will forecast a new guess:  */
  1490.         
  1491.         else if ( n_above > (n_below + n_equal) ) {  /* Guess is too low  */
  1492.         
  1493.             lower_bound = m_guess;
  1494.             t_0 = n_below + n_equal - 1;
  1495.             temp = lower_bound + (upper_bound - lower_bound) * (t_middle - t_0) / (t_1 - t_0);
  1496.             m_guess = (temp > lub) ? temp : lub;    /* Move guess at least to lub  */
  1497.         }
  1498.         else if ( n_below > (n_above + n_equal) ) {  /* Guess is too high  */
  1499.         
  1500.             upper_bound = m_guess;
  1501.             t_1 = n_below + n_equal - 1;
  1502.             temp = lower_bound + (upper_bound - lower_bound) * (t_middle - t_0) / (t_1 - t_0);
  1503.             m_guess = (temp < glb) ? temp : glb;    /* Move guess at least to glb  */
  1504.         }
  1505.         else {    /* If we get here, I made a mistake!  */
  1506.             fprintf (stderr,"GMT Fatal Error: Internal goof - please report to developers!\n");
  1507.             exit (-1);
  1508.         }
  1509.             
  1510.     } while (!finished);
  1511.     
  1512.     /* That's all, folks!  */
  1513.     return (iteration);
  1514. }
  1515.  
  1516. int mode(x, n, j, sort, mode_est)
  1517.  
  1518. double    *x, *mode_est;
  1519. int    n, j, sort;
  1520.  
  1521. {
  1522.     int    i, istop, multiplicity = 0;
  1523.     double    mid_point_sum = 0.0, length, short_length = 1.0e+30;
  1524.  
  1525.     if (sort) qsort((char *)x, n, sizeof(double), comp_double_asc);
  1526.  
  1527.     istop = n - j;
  1528.     
  1529.     for (i = 0; i < istop; i++) {
  1530.         length = x[i + j] - x[i];
  1531.         if (length < 0.0) {
  1532.             fprintf(stderr,"mode: Array not sorted in non-decreasing order.\n");
  1533.             return (-1);
  1534.         }
  1535.         else if (length == short_length) {
  1536.             multiplicity++;
  1537.             mid_point_sum += (0.5 * (x[i + j] + x[i]));
  1538.         }
  1539.         else if (length < short_length) {
  1540.             multiplicity = 1;
  1541.             mid_point_sum = (0.5 * (x[i + j] + x[i]));
  1542.             short_length = length;
  1543.         }
  1544.     }
  1545.  
  1546.     if (multiplicity - 1) mid_point_sum /= multiplicity;
  1547.     
  1548.     *mode_est = mid_point_sum;
  1549.     return (0);
  1550. }
  1551.  
  1552. int get_format (interval, format)
  1553. double interval;
  1554. char *format; {
  1555.     int i, j, ndec = 0;
  1556.     char text[80];
  1557.     
  1558.     /* Find number of decimals needed, if any, and create suitable format statement */
  1559.     
  1560.     sprintf (text, "%lg\0", interval);
  1561.     for (i = 0; text[i] && text[i] != '.'; i++);
  1562.     if (text[i]) {    /* Found a decimal point */
  1563.         for (j = i + 1; text[j]; j++);
  1564.         ndec = j - i - 1;
  1565.     }
  1566.     if (ndec > 0)
  1567.         sprintf (format, "%%.%dlf\0", ndec);
  1568.     else
  1569.         strcpy (format, "%lg");
  1570.     return (ndec);
  1571. }
  1572.  
  1573. int get_label_parameters (side, line_angle, text_angle, justify)
  1574. int side, *justify;
  1575. double line_angle, *text_angle; {
  1576.     int ok;
  1577.     
  1578.     *text_angle = line_angle;
  1579.     if (*text_angle < -90.0) *text_angle += 360.0;
  1580.     if (frame_info.horizontal && !(side%2)) *text_angle += 90.0;
  1581.     if (*text_angle >= 270.0 ) *text_angle -= 360.0;
  1582.     else if (*text_angle >= 90.0) *text_angle -= 180.0;
  1583.     
  1584.     switch (side) {
  1585.         case 0:        /* S */
  1586.             if (frame_info.horizontal)
  1587.                 *justify = 10;
  1588.             else
  1589.                 *justify = ((*text_angle) < 0.0) ? 5 : 7;
  1590.             break;
  1591.         case 1:        /* E */
  1592.             *justify = 5;
  1593.             break;
  1594.         case 2:        /* N */
  1595.             if (frame_info.horizontal)
  1596.                 *justify = 2;
  1597.             else
  1598.                 *justify = ((*text_angle) < 0.0) ? 7 : 5;
  1599.             break;
  1600.         case 3:        /* W */
  1601.             *justify = 7;
  1602.             break;
  1603.     }
  1604.     
  1605.     if (frame_info.horizontal) return (TRUE);
  1606.     
  1607.     
  1608.     switch (side) {
  1609.         case 0:        /* S */
  1610.         case 2:        /* N */
  1611.             ok = (fabs ((*text_angle)) >= gmtdefs.anot_min_angle);
  1612.             break;
  1613.         case 1:        /* E */
  1614.         case 3:        /* W */
  1615.             ok = (fabs ((*text_angle)) <= (90.0 - gmtdefs.anot_min_angle));
  1616.             break;
  1617.     }
  1618.     return (ok);
  1619. }
  1620.  
  1621. int    non_zero_winding(xp, yp, x, y, n_path)
  1622. double    xp, yp, *x, *y;
  1623. int    n_path;
  1624. {
  1625.     /* Routine returns (2) if (xp,yp) is inside the
  1626.        polygon x[n_path], y[n_path], (0) if outside,
  1627.        and (1) if exactly on the path edge.
  1628.        Uses non-zero winding rule in Adobe PostScript
  1629.        Language reference manual, section 4.6 on Painting.
  1630.        Path should have been closed first, so that
  1631.        x[n_path-1] = x[0], and y[n_path-1] = y[0].
  1632.  
  1633.        This is version 2, trying to kill a bug
  1634.        in above routine:  If point is on X edge,
  1635.        fails to discover that it is on edge.
  1636.  
  1637.        We are imagining a ray extending "up" from the
  1638.        point (xp,yp); the ray has equation x = xp, y >= yp.
  1639.        Starting with crossing_count = 0, we add 1 each time
  1640.        the path crosses the ray in the +x direction, and
  1641.        subtract 1 each time the ray crosses in the -x direction.
  1642.        After traversing the entire polygon, if (crossing_count)
  1643.        then point is inside.  We also watch for edge conditions.
  1644.  
  1645.        If two or more points on the path have x[i] == xp, then
  1646.        we have an ambiguous case, and we have to find the points
  1647.        in the path before and after this section, and check them.
  1648.        */
  1649.     
  1650.     int    i, j, k, jend, crossing_count, above;
  1651.     double    y_sect;
  1652.  
  1653.     above = FALSE;
  1654.     crossing_count = 0;
  1655.  
  1656.     /* First make sure first point in path is not a special case:  */
  1657.     j = jend = n_path - 1;
  1658.     if (x[j] == xp) {
  1659.         /* Trouble already.  We might get lucky:  */
  1660.         if (y[j] == yp) return(1);
  1661.  
  1662.         /* Go backward down the polygon until x[i] != xp:  */
  1663.         if (y[j] > yp) above = TRUE;
  1664.         i = j - 1;
  1665.         while (x[i] == xp && i > 0) {
  1666.             if (y[i] == yp) return (1);
  1667.             if (!(above) && y[i] > yp) above = TRUE;
  1668.             i--;
  1669.         }
  1670.  
  1671.         /* Now if i == 0 polygon is degenerate line x=xp;
  1672.            since we know xp,yp is inside bounding box,
  1673.            it must be on edge:  */
  1674.         if (i == 0) return(1);
  1675.  
  1676.         /* Now we want to mark this as the end, for later:  */
  1677.         jend = i;
  1678.  
  1679.         /* Now if (j-i)>1 there are some segments the point could be exactly on:  */
  1680.         for (k = i+1; k < j; k++) {
  1681.             if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) ) return (1);
  1682.         }
  1683.  
  1684.  
  1685.         /* Now we have arrived where i is > 0 and < n_path-1, and x[i] != xp.
  1686.             We have been using j = n_path-1.  Now we need to move j forward 
  1687.             from the origin:  */
  1688.         j = 1;
  1689.         while (x[j] == xp) {
  1690.             if (y[j] == yp) return (1);
  1691.             if (!(above) && y[j] > yp) above = TRUE;
  1692.             j++;
  1693.         }
  1694.  
  1695.         /* Now at the worst, j == jstop, and we have a polygon with only 1 vertex
  1696.             not at x = xp.  But now it doesn't matter, that would end us at
  1697.             the main while below.  Again, if j>=2 there are some segments to check:  */
  1698.         for (k = 0; k < j-1; k++) {
  1699.             if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) ) return (1);
  1700.         }
  1701.  
  1702.  
  1703.         /* Finally, we have found an i and j with points != xp.  If (above) we may have crossed the ray:  */
  1704.         if (above && x[i] < xp && x[j] > xp) 
  1705.             crossing_count++;
  1706.         else if (above && x[i] > xp && x[j] < xp) 
  1707.             crossing_count--;
  1708.  
  1709.         /* End nightmare scenario for x[0] == xp.  */
  1710.     }
  1711.  
  1712.     else {
  1713.         /* Get here when x[0] != xp:  */
  1714.         i = 0;
  1715.         j = 1;
  1716.         while (x[j] == xp && j < jend) {
  1717.             if (y[j] == yp) return (1);
  1718.             if (!(above) && y[j] > yp) above = TRUE;
  1719.             j++;
  1720.         }
  1721.         /* Again, if j==jend, (i.e., 0) then we have a polygon with only 1 vertex
  1722.             not on xp and we will branch out below.  */
  1723.  
  1724.         /* if ((j-i)>2) the point could be on intermediate segments:  */
  1725.         for (k = i+1; k < j-1; k++) {
  1726.             if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) ) return (1);
  1727.         }
  1728.  
  1729.         /* Now we have x[i] != xp and x[j] != xp.
  1730.             If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
  1731.             If not (above) and (j-i)>1, then we have not crossed it.
  1732.             If not (above) and j-i == 1, then we have to check the intersection point.  */
  1733.  
  1734.         if (x[i] < xp && x[j] > xp) {
  1735.             if (above) 
  1736.                 crossing_count++;
  1737.             else if ( (j-i) == 1) {
  1738.                 y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  1739.                 if (y_sect == yp) return (1);
  1740.                 if (y_sect > yp) crossing_count++;
  1741.             }
  1742.         }
  1743.         if (x[i] > xp && x[j] < xp) {
  1744.             if (above) 
  1745.                 crossing_count--;
  1746.             else if ( (j-i) == 1) {
  1747.                 y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  1748.                 if (y_sect == yp) return (1);
  1749.                 if (y_sect > yp) crossing_count--;
  1750.             }
  1751.         }
  1752.                     
  1753.         /* End easier case for x[0] != xp  */
  1754.     }
  1755.  
  1756.     /* Now MAIN WHILE LOOP begins:
  1757.         Set i = j, and search for a new j, and do as before.  */
  1758.  
  1759.     i = j;
  1760.     while (i < jend) {
  1761.         above = FALSE;
  1762.         j = i+1;
  1763.         while (x[j] == xp) {
  1764.             if (y[j] == yp) return (1);
  1765.             if (!(above) && y[j] > yp) above = TRUE;
  1766.             j++;
  1767.         }
  1768.         /* if ((j-i)>2) the point could be on intermediate segments:  */
  1769.         for (k = i+1; k < j-1; k++) {
  1770.             if ( (y[k] <= yp && y[k+1] >= yp) || (y[k] >= yp && y[k+1] <= yp) ) return (1);
  1771.         }
  1772.  
  1773.         /* Now we have x[i] != xp and x[j] != xp.
  1774.             If (above) and x[i] and x[j] on opposite sides, we are certain to have crossed the ray.
  1775.             If not (above) and (j-i)>1, then we have not crossed it.
  1776.             If not (above) and j-i == 1, then we have to check the intersection point.  */
  1777.  
  1778.         if (x[i] < xp && x[j] > xp) {
  1779.             if (above) 
  1780.                 crossing_count++;
  1781.             else if ( (j-i) == 1) {
  1782.                 y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  1783.                 if (y_sect == yp) return (1);
  1784.                 if (y_sect > yp) crossing_count++;
  1785.             }
  1786.         }
  1787.         if (x[i] > xp && x[j] < xp) {
  1788.             if (above) 
  1789.                 crossing_count--;
  1790.             else if ( (j-i) == 1) {
  1791.                 y_sect = y[i] + (y[j] - y[i]) * ( (xp - x[i]) / (x[j] - x[i]) );
  1792.                 if (y_sect == yp) return (1);
  1793.                 if (y_sect > yp) crossing_count--;
  1794.             }
  1795.         }
  1796.  
  1797.         /* That's it for this piece.  Advance i:  */
  1798.  
  1799.         i = j;
  1800.     }
  1801.  
  1802.     /* End of MAIN WHILE.  Get here when we have gone all around without landing on edge.  */
  1803.  
  1804.     if (crossing_count)
  1805.         return(2);
  1806.     else
  1807.         return(0);
  1808. }
  1809.     
  1810. /*
  1811.  * delaunay performs a Delaunay triangulation on the input data
  1812.  * and returns a list of indeces of the points for each triangle
  1813.  * found.  Algorithm translated from
  1814.  * Watson, D. F., ACORD: Automatic contouring of raw data,
  1815.  *   Computers & Geosciences, 8, 97-101, 1982.
  1816.  */
  1817.  
  1818. int delaunay (x_in, y_in, n, link)
  1819. double x_in[];    /* Input point x coordinates */
  1820. double y_in[];    /* Input point y coordinates */
  1821. int n;        /* Number of input poins */
  1822. int *link[];    /* pointer to List of point ids per triangle.  Vertices for triangle no i
  1823.            is in link[i*3], link[i*3+1], link[i*3+2] */
  1824. {
  1825.     int ix[3], iy[3];
  1826.     int i, j, ij, nuc, jt, km, id, isp, l1, l2, k, k1, jz, i2, kmt, kt, done, size;
  1827.     int *index, *istack, *x_tmp, *y_tmp;
  1828.     double det[2][3], *x_circum, *y_circum, *r2_circum, *x, *y;
  1829.     double xmin, xmax, ymin, ymax, datax, dx, dy, dsq, dd;
  1830.     
  1831.     size = 2 * n + 1;
  1832.     n += 3;
  1833.     
  1834.     index = (int *) memory (CNULL, 3 * size, sizeof (int), "delaunay");
  1835.     istack = (int *) memory (CNULL, size, sizeof (int), "delaunay");
  1836.     x_tmp = (int *) memory (CNULL, size, sizeof (int), "delaunay");
  1837.     y_tmp = (int *) memory (CNULL, size, sizeof (int), "delaunay");
  1838.     x_circum = (double *) memory (CNULL, size, sizeof (double), "delaunay");
  1839.     y_circum = (double *) memory (CNULL, size, sizeof (double), "delaunay");
  1840.     r2_circum = (double *) memory (CNULL, size, sizeof (double), "delaunay");
  1841.     x = (double *) memory (CNULL, n, sizeof (double), "delaunay");
  1842.     y = (double *) memory (CNULL, n, sizeof (double), "delaunay");
  1843.     
  1844.     x[0] = x[1] = -1.0;    x[2] = 5.0;
  1845.     y[0] = y[2] = -1.0;    y[1] = 5.0;
  1846.     x_circum[0] = y_circum[0] = 2.0;    r2_circum[0] = 18.0;
  1847.     
  1848.     ix[0] = ix[1] = 0;    ix[2] = 1;
  1849.     iy[0] = 1;    iy[1] = iy[2] = 2;
  1850.     
  1851.     for (i = 0; i < 3; i++) index[i] = i;
  1852.     for (i = 0; i < size; i++) istack[i] = i;
  1853.     
  1854.     xmin = ymin = 1.0e100;    xmax = ymax = -1.0e100;
  1855.     
  1856.     for (i = 3, j = 0; i < n; i++, j++) {    /* Copy data and get extremas */
  1857.         x[i] = x_in[j];
  1858.         y[i] = y_in[j];
  1859.         if (x[i] > xmax) xmax = x[i];
  1860.         if (x[i] < xmin) xmin = x[i];
  1861.         if (y[i] > ymax) ymax = y[i];
  1862.         if (y[i] < ymin) ymin = y[i];
  1863.     }
  1864.     
  1865.     /* Normalize data */
  1866.     
  1867.     datax = 1.0 / MAX (xmax - xmin, ymax - ymin);
  1868.     
  1869.     for (i = 3; i < n; i++) {
  1870.         x[i] = (x[i] - xmin) * datax;
  1871.         y[i] = (y[i] - ymin) * datax;
  1872.     }
  1873.     
  1874.     isp = id = 1;
  1875.     for (nuc = 3; nuc < n; nuc++) {
  1876.     
  1877.         km = 0;
  1878.         
  1879.         for (jt = 0; jt < isp; jt++) {    /* loop through established 3-tuples */
  1880.         
  1881.             ij = 3 * jt;
  1882.             
  1883.             /* Test if new data is within the jt circumcircle */
  1884.             
  1885.             dx = x[nuc] - x_circum[jt];
  1886.             if ((dsq = r2_circum[jt] - dx * dx) < 0.0) continue;
  1887.             dy = y[nuc] - y_circum[jt];
  1888.             if ((dsq -= dy * dy) < 0.0) continue;
  1889.  
  1890.             /* Delete this 3-tuple but save its edges */
  1891.             
  1892.             id--;
  1893.             istack[id] = jt;
  1894.             
  1895.             /* Add edges to x/y_tmp but delete if already listed */
  1896.             
  1897.             for (i = 0; i < 3; i++) {
  1898.                 l1 = ix[i];
  1899.                 l2 = iy[i];
  1900.                 if (km > 0) {
  1901.                     kmt = km;
  1902.                     for (j = 0, done = FALSE; !done && j < kmt; j++) {
  1903.                         if (index[ij+l1] != x_tmp[j]) continue;
  1904.                         if (index[ij+l2] != y_tmp[j]) continue;
  1905.                         km--;
  1906.                         if (j >= km) {
  1907.                             done = TRUE;
  1908.                             continue;
  1909.                         }
  1910.                         for (k = j; k < km; k++) {
  1911.                             k1 = k + 1;
  1912.                             x_tmp[k] = x_tmp[k1];
  1913.                             y_tmp[k] = y_tmp[k1];
  1914.                             done = TRUE;
  1915.                         }
  1916.                     }
  1917.                 }
  1918.                 else
  1919.                     done = FALSE;
  1920.                 if (!done) {
  1921.                     x_tmp[km] = index[ij+l1];
  1922.                     y_tmp[km] = index[ij+l2];
  1923.                     km++;
  1924.                 }
  1925.             }
  1926.         }
  1927.             
  1928.         /* Form new 3-tuples */
  1929.             
  1930.         for (i = 0; i < km; i++) {
  1931.             kt = istack[id];
  1932.             ij = 3 * kt;
  1933.             id++;
  1934.                 
  1935.             /* Calculate the circumcircle center and radius
  1936.              * squared of points ktetr[i,*] and place in tetr[kt,*] */
  1937.                 
  1938.             for (jz = 0; jz < 2; jz++) {
  1939.                 i2 = (jz == 0) ? x_tmp[i] : y_tmp[i];
  1940.                 det[jz][0] = x[i2] - x[nuc];
  1941.                 det[jz][1] = y[i2] - y[nuc];
  1942.                 det[jz][2] = 0.5 * (det[jz][0] * (x[i2] + x[nuc]) + det[jz][1] * (y[i2] + y[nuc]));
  1943.             }
  1944.             dd = 1.0 / (det[0][0] * det[1][1] - det[0][1] * det[1][0]);
  1945.             x_circum[kt] = (det[0][2] * det[1][1] - det[1][2] * det[0][1]) * dd;
  1946.             y_circum[kt] = (det[0][0] * det[1][2] - det[1][0] * det[0][2]) * dd;
  1947.             dx = x[nuc] - x_circum[kt];
  1948.             dy = y[nuc] - y_circum[kt];
  1949.             r2_circum[kt] = dx * dx + dy * dy;
  1950.             index[ij++] = x_tmp[i];
  1951.             index[ij++] = y_tmp[i];
  1952.             index[ij] = nuc;
  1953.         }
  1954.         isp += 2;
  1955.     }
  1956.     
  1957.     for (jt = i = 0; jt < isp; jt++) {
  1958.         ij = 3 * jt;
  1959.         if (index[ij] < 3 || r2_circum[jt] > 1.0) continue;
  1960.         index[i++] = index[ij++] - 3;
  1961.         index[i++] = index[ij++] - 3;
  1962.         index[i++] = index[ij++] - 3;
  1963.     }
  1964.     
  1965.     index = (int *) memory ((char *)index, i, sizeof (int), "delaunay");
  1966.     *link = index;
  1967.     
  1968.     free ((char *)istack);
  1969.     free ((char *)x_tmp);
  1970.     free ((char *)y_tmp);
  1971.     free ((char *)x_circum);
  1972.     free ((char *)y_circum);
  1973.     free ((char *)r2_circum);
  1974.     free ((char *)x);
  1975.     free ((char *)y);
  1976.     
  1977.     return (i/3);
  1978. }
  1979.  
  1980. /*
  1981.  * grd_init initializes a grd header to default values and copies the
  1982.  * command line to the header variable command
  1983.  */
  1984.  
  1985. int grd_init (header, argc, argv, update)
  1986. struct GRD_HEADER *header; 
  1987. int argc;
  1988. char **argv;
  1989. BOOLEAN update; {    /* TRUE if we only want to update command line */
  1990.     int i, len;
  1991.     
  1992.     /* Always update command line history */
  1993.     
  1994.     memset (header->command, 0, 320);
  1995.     if (argc > 0) {
  1996.         strcpy (header->command, argv[0]);
  1997.         len = strlen (header->command);
  1998.         for (i = 1; len < 320 && i < argc; i++) {
  1999.             len += strlen (argv[i]) + 1;
  2000.             if (len > 320) continue;
  2001.             strcat (header->command, " ");
  2002.             strcat (header->command, argv[i]);
  2003.         }
  2004.         header->command[len] = 0;
  2005.     }
  2006.  
  2007.     if (update) return;    /* Leave other variables unchanged */
  2008.     
  2009.     /* Here we initialize the variables to default settings */
  2010.     
  2011.     header->x_min = header->x_max = 0.0;
  2012.     header->y_min = header->y_max = 0.0;
  2013.     header->z_min = header->z_max = 0.0;
  2014.     header->x_inc = header->y_inc    = 0.0;
  2015.     header->z_scale_factor        = 1.0;
  2016.     header->z_add_offset        = 0.0;
  2017.     header->nx = header->ny        = 0;
  2018.     header->node_offset        = FALSE;
  2019.     
  2020.     memset (header->x_units, 0, 80);
  2021.     memset (header->y_units, 0, 80);
  2022.     memset (header->z_units, 0, 80);
  2023.     strcpy (header->x_units, "user_x_unit");
  2024.     strcpy (header->y_units, "user_y_unit");
  2025.     strcpy (header->z_units, "user_z_unit");
  2026.     memset (header->title, 0, 80);
  2027.     memset (header->remark, 0, 160);
  2028. }
  2029.  
  2030. int grd_shift (header, grd, shift)
  2031. struct GRD_HEADER *header;
  2032. float grd[];
  2033. double shift;        /* Rotate geographical, global grid in e-w direction */
  2034. {
  2035.     /* This function will shift a grid by shift degrees */
  2036.     
  2037.     int i, j, k, ij, nc, n_shift, width;
  2038.     float *tmp;
  2039.     
  2040.     tmp = (float *) memory (CNULL, header->nx, sizeof (float), "grd_shift");
  2041.     
  2042.     n_shift = rint (shift / header->x_inc);
  2043.     width = (header->node_offset) ? header->nx : header->nx - 1;
  2044.     nc = header->nx * sizeof (float);
  2045.     
  2046.     for (j = ij = 0; j < header->ny; j++, ij += header->nx) {
  2047.         for (i = 0; i < width; i++) {
  2048.             k = (i - n_shift) % width;
  2049.             if (k < 0) k += header->nx;
  2050.             tmp[k] = grd[ij+i];
  2051.         }
  2052.         if (!header->node_offset) tmp[width] = tmp[0];
  2053.         memcpy ((char *)&grd[ij], (char *)tmp, nc);
  2054.     }
  2055.     
  2056.     /* Shift boundaries */
  2057.     
  2058.     header->x_min += shift;
  2059.     header->x_max += shift;
  2060.     if (header->x_max < 0.0) {
  2061.         header->x_min += 360.0;
  2062.         header->x_max += 360.0;
  2063.     }
  2064.     else if (header->x_max > 360.0) {
  2065.         header->x_min -= 360.0;
  2066.         header->x_max -= 360.0;
  2067.     }
  2068.     
  2069.     free ((char *) tmp);
  2070. }
  2071.  
  2072. /* grd_setregion determines what wesn should be passed to grd_read.  It
  2073.   does so by using project_info.w,e,s,n which have been set correctly
  2074.   by map_setup. */
  2075.  
  2076. int grd_setregion (h, xmin, xmax, ymin, ymax)
  2077. struct GRD_HEADER *h;
  2078. double *xmin, *xmax, *ymin, *ymax; {
  2079.  
  2080.     /* Round off to nearest multiple of the grid spacing.  This should
  2081.        only affect the numbers when oblique projections or -R...r has been used */
  2082.  
  2083.     *xmin = floor (project_info.w / h->x_inc) * h->x_inc;
  2084.     *xmax = ceil (project_info.e / h->x_inc) * h->x_inc;
  2085.     if (fabs (h->x_max - h->x_min - 360.0) > SMALL) {    /* Not a periodic grid */
  2086.         if (*xmin < h->x_min) *xmin = h->x_min;
  2087.         if (*xmax > h->x_max) *xmax = h->x_max;
  2088.     }
  2089.     *ymin = floor (project_info.s / h->y_inc) * h->y_inc;
  2090.     if (*ymin < h->y_min) *ymin = h->y_min;
  2091.     *ymax = ceil (project_info.n / h->y_inc) * h->y_inc;
  2092.     if (*ymax > h->y_max) *ymax = h->y_max;
  2093. }
  2094.  
  2095. /* code for bicubic rectangle and bilinear interpolation of grd files
  2096.  *
  2097.  * Author:    Walter H F Smith
  2098.  * Date:    23 September 1993
  2099. */
  2100.  
  2101.  
  2102. void    bcr_init(grd, pad, bilinear)
  2103. struct GRD_HEADER *grd;
  2104. int    pad[];    /* padding on grid, as given to read_grd2 */
  2105. int    bilinear;    /* T/F we only want bilinear  */
  2106. {
  2107.     /* Initialize i,j so that they cannot look like they have been used:  */
  2108.     bcr.i = -10;
  2109.     bcr.j = -10;
  2110.  
  2111.     /* Initialize bilinear:  */
  2112.     bcr.bilinear = bilinear;
  2113.  
  2114.     /* Initialize ioff, joff, mx, my according to grd and pad:  */
  2115.     bcr.ioff = pad[0];
  2116.     bcr.joff = pad[3];
  2117.     bcr.mx = grd->nx + pad[0] + pad[1];
  2118.     bcr.my = grd->ny + pad[2] + pad[3];
  2119.  
  2120.     /* Initialize rx_inc, ry_inc, and offset:  */
  2121.     bcr.rx_inc = 1.0 / grd->x_inc;
  2122.     bcr.ry_inc = 1.0 / grd->y_inc;
  2123.     bcr.offset = (grd->node_offset) ? 0.5 : 0.0;
  2124.  
  2125.     /* Initialize ij_move:  */
  2126.     bcr.ij_move[0] = 0;
  2127.     bcr.ij_move[1] = 1;
  2128.     bcr.ij_move[2] = -bcr.mx;
  2129.     bcr.ij_move[3] = 1 - bcr.mx;
  2130.  
  2131.     return;
  2132. }
  2133.  
  2134. int    set_bcr_boundaries (grd, f)
  2135. struct GRD_HEADER *grd;
  2136. float f[];
  2137. {
  2138.     int i, j, ne, nw, se, sw;
  2139.     
  2140.     /* Get inside corner indeces */
  2141.     
  2142.     nw = bcr.joff * bcr.mx + bcr.ioff;
  2143.     ne = nw + grd->nx - 1;
  2144.     sw = (grd->ny + bcr.joff - 1) * bcr.mx + bcr.ioff;
  2145.     se = sw + grd->nx - 1;
  2146.     
  2147.     /* Set Natural spline boundary conditions  */
  2148.     
  2149.     for (i = 0; i < grd->nx; i++) {
  2150.         j = nw - bcr.mx + i;        /* Set top row  */
  2151.         f[j] = f[j+bcr.mx] + (f[j+bcr.mx] - f[j+bcr.mx+bcr.mx]);
  2152.         j = sw + bcr.mx + i;    /* Bottom row  */
  2153.         f[j] = f[j-bcr.mx] + (f[j-bcr.mx] - f[j-bcr.mx-bcr.mx]);
  2154.     }
  2155.     for (i = 0; i < grd->ny; i++) {
  2156.         j = nw - 1 + i * bcr.mx;    /* Left column  */
  2157.         f[j] = f[j+1] + (f[j+1] - f[j+2]);
  2158.         j += (grd->nx + 1);        /* Right column  */
  2159.         f[j] = f[j-1] + (f[j-1] - f[j-2]);
  2160.     }
  2161.     
  2162.     f[nw-bcr.mx-1] = f[nw-bcr.mx] - (f[nw] - f[nw-1]);    /* NW corner  */
  2163.     f[ne-bcr.mx+1] = f[ne-bcr.mx] + (f[nw+1] - f[nw]);    /* NE corner  */
  2164.     f[sw+bcr.mx-1] = f[sw+bcr.mx] - (f[sw] - f[sw-1]);    /* SW corner  */
  2165.     f[se+bcr.mx+1] = f[se+1] - (f[se] - f[se+bcr.mx]);    /* SE corner:  */
  2166. }
  2167.  
  2168. void    get_bcr_cardinals(x, y)
  2169. double    x;
  2170. double    y;
  2171. {
  2172.     /* Given x,y compute the cardinal functions.  Note x,y should be in
  2173.      * normalized range, usually [0,1) but sometimes a little outside this.  */
  2174.  
  2175.     double    xcf[2][2], ycf[2][2], tsq, tm1, tm1sq;
  2176.     int    vertex, verx, very, value, valx, valy;
  2177.  
  2178.     if (bcr.bilinear) {
  2179.         bcr.bl_basis[0] = (1.0 - x)*(1.0 - y);
  2180.         bcr.bl_basis[1] = x*(1.0 - y);
  2181.         bcr.bl_basis[2] = y*(1.0 - x);
  2182.         bcr.bl_basis[3] = x*y;
  2183.  
  2184.         return;
  2185.     }
  2186.  
  2187.     /* Get here when we need to do bicubic functions:  */
  2188.     tsq = x*x;
  2189.     tm1 = x - 1.0;
  2190.     tm1sq = tm1 * tm1;
  2191.     xcf[0][0] = (2.0*x + 1.0) * tm1sq;
  2192.     xcf[0][1] = x * tm1sq;
  2193.     xcf[1][0] = tsq * (3.0 - 2.0*x);
  2194.     xcf[1][1] = tsq * tm1;
  2195.  
  2196.     tsq = y*y;
  2197.     tm1 = y - 1.0;
  2198.     tm1sq = tm1 * tm1;
  2199.     ycf[0][0] = (2.0*y + 1.0) * tm1sq;
  2200.     ycf[0][1] = y * tm1sq;
  2201.     ycf[1][0] = tsq * (3.0 - 2.0*y);
  2202.     ycf[1][1] = tsq * tm1;
  2203.  
  2204.     for (vertex = 0; vertex < 4; vertex++) {
  2205.         verx = vertex%2;
  2206.         very = vertex/2;
  2207.         for (value = 0; value < 4; value++) {
  2208.             valx = value%2;
  2209.             valy = value/2;
  2210.  
  2211.             bcr.bcr_basis[vertex][value] = xcf[verx][valx] * ycf[very][valy];
  2212.         }
  2213.     }
  2214.     return;
  2215. }
  2216.  
  2217. void    get_bcr_ij(grd, xx, yy, ii, jj)
  2218. struct GRD_HEADER *grd;
  2219. double    xx, yy;
  2220. int    *ii, *jj;
  2221. {
  2222.     /* Given xx, yy in user's grdfile x and y units (not normalized),
  2223.        set ii,jj to the point to be used for the bqr origin.  This 
  2224.        should have jj in the range 1 grd->ny-1 and ii in the range 0 to
  2225.        grd->nx-2, so that the north and east edges will not have the
  2226.        origin on them.
  2227.        This function should NOT be called unless xx,yy are within the
  2228.        valid range of the grid.  */
  2229.  
  2230.     int    i, j;
  2231.  
  2232.     i = floor((xx-grd->x_min)*bcr.rx_inc - bcr.offset);
  2233.     if (i < 0 ) i = 0;
  2234.     if (i > grd->nx-2) i = grd->nx-2;
  2235.     j = ceil ((grd->y_max-yy)*bcr.ry_inc - bcr.offset);
  2236.     if (j < 1) j = 1;
  2237.     if (j > grd->ny-1) j = grd->ny-1;
  2238.     
  2239.     *ii = i;
  2240.     *jj = j;
  2241.  
  2242.     return;
  2243. }
  2244.     
  2245. void    get_bcr_xy(grd, xx, yy, x, y)
  2246. struct GRD_HEADER *grd;
  2247. double    xx, yy;
  2248. double    *x, *y;
  2249. {
  2250.     /* Given xx, yy in user's grdfile x and y units (not normalized),
  2251.        use the bcr.i and bcr.j to find x,y (normalized) which are the
  2252.        coordinates of xx,yy in the bcr rectangle.  */
  2253.  
  2254.     double    xorigin, yorigin;
  2255.  
  2256.     xorigin = (bcr.i + bcr.offset)*grd->x_inc + grd->x_min;
  2257.     yorigin = grd->y_max - (bcr.j + bcr.offset)*grd->y_inc;
  2258.  
  2259.     *x = (xx - xorigin) * bcr.rx_inc;
  2260.     *y = (yy - yorigin) * bcr.ry_inc;
  2261.  
  2262.     return;
  2263. }
  2264.  
  2265. void    get_bcr_nodal_values(z, ii, jj)
  2266. float    z[];
  2267. int    ii, jj;
  2268. {
  2269.     /* ii, jj is the point we want to use now, which is different from
  2270.        the bcr.i, bcr.j point we used last time this function was called.
  2271.        If (nan_condition == FALSE) && abs(ii-bcr.i) < 2 && abs(jj-bcr.j)
  2272.        < 2 then we can reuse some previous results.  */
  2273.  
  2274.     int    i, valstop, vertex, ij, ij_origin, k0, k1, k2, k3, dontneed[4];
  2275.  
  2276.     /* whattodo[vertex] says which vertices are previously known.  */
  2277.     for (i = 0; i < 4; i++) dontneed[i] = FALSE;
  2278.  
  2279.     valstop = (bcr.bilinear) ? 1 : 4;
  2280.  
  2281.     if (!(bcr.nan_condition) && (abs(ii-bcr.i) < 2 && abs(jj-bcr.j) < 2) ) {
  2282.         /* There was no nan-condition last time and we can use some
  2283.             previously computed results.  */
  2284.         switch (ii-bcr.i) {
  2285.             case 1:
  2286.                 /* We have moved to the right ...  */
  2287.                 switch (jj-bcr.j) {
  2288.                     case -1:
  2289.                         /* ... and up.  New 0 == old 3  */
  2290.                         dontneed[0] = TRUE;
  2291.                         for (i = 0; i < valstop; i++)
  2292.                             bcr.nodal_value[0][i] = bcr.nodal_value[3][i];
  2293.                         break;
  2294.                     case 0:
  2295.                         /* ... horizontally.  New 0 == old 1; New 2 == old 3  */
  2296.                         dontneed[0] = dontneed[2] = TRUE;
  2297.                         for (i = 0; i < valstop; i++) {
  2298.                             bcr.nodal_value[0][i] = bcr.nodal_value[1][i];
  2299.                             bcr.nodal_value[2][i] = bcr.nodal_value[3][i];
  2300.                         }
  2301.                         break;
  2302.                     case 1:
  2303.                         /* ... and down.  New 2 == old 1  */
  2304.                         dontneed[2] = TRUE;
  2305.                         for (i = 0; i < valstop; i++)
  2306.                             bcr.nodal_value[2][i] = bcr.nodal_value[1][i];
  2307.                         break;
  2308.                 }
  2309.                 break;
  2310.             case 0:
  2311.                 /* We have moved only ...  */
  2312.                 switch (jj-bcr.j) {
  2313.                     case -1:
  2314.                         /* ... up.  New 0 == old 2; New 1 == old 3  */
  2315.                         dontneed[0] = dontneed[1] = TRUE;
  2316.                         for (i = 0; i < valstop; i++) {
  2317.                             bcr.nodal_value[0][i] = bcr.nodal_value[2][i];
  2318.                             bcr.nodal_value[1][i] = bcr.nodal_value[3][i];
  2319.                         }
  2320.                         break;
  2321.                     case 0:
  2322.                         /* ... not at all.  This shouldn't happen  */
  2323.                         return;
  2324.                         break;
  2325.                     case 1:
  2326.                         /* ... down.  New 2 == old 0; New 3 == old 1  */
  2327.                         dontneed[2] = dontneed[3] = TRUE;
  2328.                         for (i = 0; i < valstop; i++) {
  2329.                             bcr.nodal_value[2][i] = bcr.nodal_value[0][i];
  2330.                             bcr.nodal_value[3][i] = bcr.nodal_value[1][i];
  2331.                         }
  2332.                         break;
  2333.                 }
  2334.                 break;
  2335.             case -1:
  2336.                 /* We have moved to the left ...  */
  2337.                 switch (jj-bcr.j) {
  2338.                     case -1:
  2339.                         /* ... and up.  New 1 == old 2  */
  2340.                         dontneed[1] = TRUE;
  2341.                         for (i = 0; i < valstop; i++)
  2342.                             bcr.nodal_value[1][i] = bcr.nodal_value[2][i];
  2343.                         break;
  2344.                     case 0:
  2345.                         /* ... horizontally.  New 1 == old 0; New 3 == old 2  */
  2346.                         dontneed[1] = dontneed[3] = TRUE;
  2347.                         for (i = 0; i < valstop; i++) {
  2348.                             bcr.nodal_value[1][i] = bcr.nodal_value[0][i];
  2349.                             bcr.nodal_value[3][i] = bcr.nodal_value[2][i];
  2350.                         }
  2351.                         break;
  2352.                     case 1:
  2353.                         /* ... and down.  New 3 == old 0  */
  2354.                         dontneed[3] = TRUE;
  2355.                         for (i = 0; i < valstop; i++)
  2356.                             bcr.nodal_value[3][i] = bcr.nodal_value[0][i];
  2357.                         break;
  2358.                 }
  2359.                 break;
  2360.         }
  2361.     }
  2362.         
  2363.     /* When we get here, we are ready to look for new values (and possibly derivatives)  */
  2364.  
  2365.     ij_origin = (jj + bcr.joff) * bcr.mx + (ii + bcr.ioff);
  2366.     bcr.i = ii;
  2367.     bcr.j = jj;
  2368.  
  2369.     for (vertex = 0; vertex < 4; vertex++) {
  2370.  
  2371.         if (dontneed[vertex]) continue;
  2372.  
  2373.         ij = ij_origin + bcr.ij_move[vertex];
  2374.         if (bad_float ((double)z[ij])) {
  2375.             bcr.nan_condition = TRUE;
  2376.             return;
  2377.         }
  2378.         bcr.nodal_value[vertex][0] = (double)z[ij];
  2379.  
  2380.         if (bcr.bilinear) continue;
  2381.  
  2382.         /* Get dz/dx:  */
  2383.         if (bad_float ((double)z[ij+1]) || bad_float((double)z[ij-1]) ){
  2384.             bcr.nan_condition = TRUE;
  2385.             return;
  2386.         }
  2387.         bcr.nodal_value[vertex][1] = 0.5 * (z[ij+1] - z[ij-1]);
  2388.  
  2389.         /* Get dz/dy:  */
  2390.         if (bad_float ((double)z[ij+bcr.mx]) || bad_float((double)z[ij-bcr.mx]) ){
  2391.             bcr.nan_condition = TRUE;
  2392.             return;
  2393.         }
  2394.         bcr.nodal_value[vertex][2] = 0.5 * (z[ij-bcr.mx] - z[ij+bcr.mx]);
  2395.  
  2396.         /* Get d2z/dxdy:  */
  2397.         k0 = ij + bcr.mx - 1;
  2398.         if (bad_float ((double)z[k0])) {
  2399.             bcr.nan_condition = TRUE;
  2400.             return;
  2401.         }
  2402.         k1 = k0 + 2;
  2403.         if (bad_float ((double)z[k1])) {
  2404.             bcr.nan_condition = TRUE;
  2405.             return;
  2406.         }
  2407.         k2 = ij - bcr.mx - 1;
  2408.         if (bad_float ((double)z[k2])) {
  2409.             bcr.nan_condition = TRUE;
  2410.             return;
  2411.         }
  2412.         k3 = k2 + 2;
  2413.         if (bad_float ((double)z[k3])) {
  2414.             bcr.nan_condition = TRUE;
  2415.             return;
  2416.         }
  2417.         bcr.nodal_value[vertex][3] = 0.25 * ( (z[k3] - z[k2]) - (z[k1] - z[k0]) );
  2418.     }
  2419.  
  2420.     /* If we get here, we have not found any NaNs that would screw us up.  */
  2421.     bcr.nan_condition = FALSE;
  2422.     return;
  2423. }
  2424.  
  2425. double    get_bcr_z(grd, xx, yy, data)
  2426. struct GRD_HEADER *grd;
  2427. double    xx, yy;
  2428. float    data[];
  2429. {
  2430.     /* Given xx, yy in user's grdfile x and y units (not normalized),
  2431.        this routine returns the desired interpolated value (bilinear
  2432.        or bicubic) at xx, yy.  */
  2433.  
  2434.     int    i, j, vertex, value;
  2435.     double    x, y, retval;
  2436.     
  2437.     if (xx < grd->x_min || xx > grd->x_max) return(gmt_NaN);
  2438.     if (yy < grd->y_min || yy > grd->y_max) return(gmt_NaN);
  2439.  
  2440.     if (fmod (xx, grd->x_inc) == 0.0 && fmod (yy, grd->y_inc) == 0.0) {    /* Special case, simply return this node value */
  2441.             i = floor((xx-grd->x_min)*bcr.rx_inc - bcr.offset);
  2442.             j = ceil ((grd->y_max-yy)*bcr.ry_inc - bcr.offset);
  2443.         return ((double)data[(j+bcr.joff)*bcr.mx+i+bcr.ioff]);
  2444.     }
  2445.  
  2446.     get_bcr_ij(grd, xx, yy, &i, &j);
  2447.  
  2448.     if (i != bcr.i || j != bcr.j)
  2449.         get_bcr_nodal_values(data, i, j);
  2450.  
  2451.     if (bcr.nan_condition) return(gmt_NaN);
  2452.  
  2453.     get_bcr_xy(grd, xx, yy, &x, &y);
  2454.  
  2455.     if (x == 0.0) {
  2456.         if (y == 0.0)
  2457.             return(bcr.nodal_value[0][0]);
  2458.         if (y == 1.0)
  2459.             return(bcr.nodal_value[2][0]);
  2460.     }
  2461.     if (x == 1.0) {
  2462.         if (y == 0.0)
  2463.             return(bcr.nodal_value[1][0]);
  2464.         if (y == 1.0)
  2465.             return(bcr.nodal_value[3][0]);
  2466.     }
  2467.  
  2468.     get_bcr_cardinals(x, y);
  2469.  
  2470.     retval = 0.0;
  2471.     if (bcr.bilinear) {
  2472.         for (vertex = 0; vertex < 4; vertex++) {
  2473.             retval += (bcr.nodal_value[vertex][0] * bcr.bl_basis[vertex]);
  2474.         }
  2475.         return(retval);
  2476.     }
  2477.     for (vertex = 0; vertex < 4; vertex++) {
  2478.         for (value = 0; value < 4; value++) {
  2479.             retval += (bcr.nodal_value[vertex][value]*bcr.bcr_basis[vertex][value]);
  2480.         }
  2481.     }
  2482.     return(retval);
  2483. }
  2484.  
  2485. /* Table I/O routines for ascii and binary io */
  2486.  
  2487. int ascii_input (fp, n, ptr)
  2488. FILE *fp;
  2489. int n;
  2490. double **ptr;
  2491. {
  2492.     char line[BUFSIZ], *p;
  2493.     int i = 0;
  2494.     double val;
  2495.     
  2496.     p = fgets (line, BUFSIZ, fp);
  2497.     if (!p) return (-1);
  2498.     
  2499.     p = strtok(line, " \t,");
  2500.     while (p && i < n) {
  2501.         gmt_data[i] = (sscanf (p, "%lf", &val) == 1) ? val : gmt_NaN;
  2502.         p = strtok(CNULL, " \t,");
  2503.         i++;
  2504.     }
  2505.     
  2506.     *ptr = gmt_data;
  2507.     return (i);
  2508. }
  2509.  
  2510. int bin_double_input (fp, n, ptr)
  2511. FILE *fp;
  2512. int n;
  2513. double **ptr;
  2514.  
  2515. {
  2516.     if (fread((char *) gmt_data, sizeof (double), n, fp) != n) return (EOF);
  2517.     *ptr = gmt_data;
  2518.     return (n);
  2519. }
  2520.     
  2521.  
  2522. int bin_float_input (fp, n, ptr)
  2523. FILE *fp;
  2524. int n;
  2525. double **ptr;
  2526. {
  2527.     int i;
  2528.     static float gmt_f[BUFSIZ];
  2529.     
  2530.     if (fread((char *) gmt_f, sizeof (float), n, fp) != n) return (EOF);
  2531.     for (i = 0; i < n; i++) gmt_data[i] = (double)gmt_f[i];
  2532.     *ptr = gmt_data;
  2533.     return (n);
  2534. }
  2535.     
  2536. int ascii_output (fp, n, ptr)
  2537. FILE *fp;
  2538. int n;
  2539. double *ptr;
  2540. {
  2541.     int i;
  2542.     
  2543.     n--;
  2544.     for (i = 0; i < n; i++) {
  2545.         (bad_float (ptr[i])) ? fprintf (fp, "NaN\t") : (fprintf (fp, gmtdefs.d_format, ptr[i]), putc ('\t', fp));
  2546.     }
  2547.     (bad_float (ptr[n])) ? fprintf (fp, "NaN\n") : (fprintf (fp, gmtdefs.d_format, ptr[n]), putc ('\n', fp));
  2548. }
  2549.  
  2550. int bin_double_output (fp, n, ptr)
  2551. FILE *fp;
  2552. int n;
  2553. double *ptr;
  2554.  
  2555. {
  2556.     fwrite((char *) ptr, sizeof (double), n, fp);
  2557. }
  2558.     
  2559. int bin_float_output (fp, n, ptr)
  2560. FILE *fp;
  2561. int n;
  2562. double *ptr;
  2563. {
  2564.     int i;
  2565.     static float gmt_f[BUFSIZ];
  2566.     
  2567.     for (i = 0; i < n; i++) gmt_f[i] = (float)ptr[i];
  2568.     fwrite((char *) gmt_f, sizeof (float), n, fp);
  2569. }
  2570.  
  2571.