home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / surface.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-28  |  60.8 KB  |  1,873 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)surface.c    2.27  4/28/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * surface.c:  a gridding program.
  9.  * reads xyz triples and fits a surface to the data.
  10.  * surface satisfies (1 - T) D4 z - T D2 z = 0,
  11.  * where D4 is the biharmonic operator,
  12.  * D2 is the Laplacian,
  13.  * and T is a "tension factor" between 0 and 1.
  14.  * End member T = 0 is the classical minimum curvature
  15.  * surface.  T = 1 gives a harmonic surface.  Use T = 0.25
  16.  * or so for potential data; something more for topography.
  17.  *
  18.  * Program includes overrelaxation for fast convergence and
  19.  * automatic optimal grid factorization.
  20.  *
  21.  * See Smith & Wessel (Geophysics, 3, 293-305, 1990) for details.
  22.  *
  23.  * Authors: Walter H. F. Smith and Paul Wessel
  24.  * Date: April, 1988.
  25.  *
  26.  * Version 3.0 1 Nov 1988:  New argument switches standardized.
  27.  *    W. H. F. Smith
  28.  *
  29.  * Version 3.2 5 Dec 1988:    Search radius default = 0.0
  30.  *    to skip this step - unnecessary now that planes are
  31.  *    removed, unless grid starts at 1 for diabolical prime lattice.
  32.  *    Ver 3.1 removed best fit plane, undocumented changes.
  33.  *    W. H. F. Smith
  34.  *
  35.  * Version 4.0 -- 5 Dec 1988:    The scratch file has been eliminated,
  36.  *    resulting in a small gain in speed.   It is replaced by the 
  37.  *    in core struct BRIGGS.  Division by xinc or grid_xinc has been
  38.  *    replaced with multiplies by reciprocals, which also made a small
  39.  *    improvement in speed.  The computation of i, j to index a point
  40.  *    was changed from rint(x) to floor(x + 0.5) to make it follow the
  41.  *    convention of blockmean.  This is actually significant at times.
  42.  *    The significance was discovered when the routine "throw_away_unusables"
  43.  *    was written.  This routine and struct BRIGGS are the major differences
  44.  *    between V3.2 and V4.0.  Modifications by w.h.f. smith
  45.  *
  46.  * V 4.1 -- 26 July, 1989:    W.H.F. Smith fixed up the arg loop and usage
  47.  *    message, and added automatic scaling of the range of z values as an
  48.  *    experiment.  Often this improves the accuracy/convergence of numerical
  49.  *    work; simple testing today suggests it makes no difference, but I will
  50.  *    leave it in.
  51.  *
  52.  * V 4.2 -- 4 June, 1991-1995:    P. Wessel upgraded to GMT v2.0 grdfile i/o.
  53.  *    Added feature to constrain solution to be within lower and/or
  54.  *    upper bounds.  Bounds can be min/max input data value, another
  55.  *    value outside the input data range, or be provided by a grdfile
  56.  *    whose values all are outside the input data range.
  57.  *
  58.  * V 4.3 -- 26 February, 1992:    W. H. F. Smith added option to suggest better
  59.  *    dimensions, and fixed bug in -L option to -Lu -Ll so that full paths
  60.  *    to lower/upper bound files may be specified.
  61.  *
  62.  *
  63.  */
  64.  
  65. #include "gmt.h"
  66.  
  67. #define OUTSIDE 2000000000    /* Index number indicating data is outside useable area */
  68.  
  69. int n_alloc = GMT_CHUNK;
  70. int npoints=0;            /* Number of data points */
  71. int nx=0;            /* Number of nodes in x-dir. */
  72. int ny=0;            /* Number of nodes in y-dir. (Final grid) */
  73. int    mx = 0;
  74. int    my = 0;
  75. int    ij_sw_corner, ij_se_corner, ij_nw_corner, ij_ne_corner;
  76. int block_nx;            /* Number of nodes in x-dir for a given grid factor */
  77. int block_ny;            /* Number of nodes in y-dir for a given grid factor */
  78. int max_iterations=250;        /* Max iter per call to iterate */
  79. int total_iterations = 0;
  80. int grid, old_grid;        /* Node spacings  */
  81. int    grid_east;
  82. int n_fact = 0;            /* Number of factors in common (ny-1, nx-1) */
  83. int factors[32];        /* Array of common factors */
  84. int long_verbose = FALSE;
  85. int n_empty;            /* No of unconstrained nodes at initialization  */
  86. int set_low = 0;        /* 0 unconstrained,1 = by min data value, 2 = by user value */
  87. int set_high = 0;        /* 0 unconstrained,1 = by max data value, 2 = by user value */
  88. int constrained = FALSE;    /* TRUE if set_low or set_high is TRUE */
  89. double low_limit, high_limit;    /* Constrains on range of solution */
  90. double xmin, xmax, ymin, ymax;    /* minmax coordinates */
  91. float *lower, *upper;        /* arrays for minmax values, if set */
  92. double xinc, yinc;        /* Size of each grid cell (final size) */
  93. double grid_xinc, grid_yinc;    /* size of each grid cell for a given grid factor */
  94. double r_xinc, r_yinc, r_grid_xinc, r_grid_yinc;    /* Reciprocals  */
  95. double converge_limit = 0.0;    /* Convergence limit */
  96. double radius = 0.0;            /* Search radius for initializing grid  */
  97. double    tension = 0.0;        /* Tension parameter on the surface  */
  98. double    boundary_tension = 0.0;
  99. double    interior_tension = 0.0;
  100. double    a0_const_1, a0_const_2;    /* Constants for off grid point equation  */
  101. double    e_2, e_m2, one_plus_e2;
  102. double eps_p2, eps_m2, two_plus_ep2, two_plus_em2;
  103. double    x_edge_const, y_edge_const;
  104. double    epsilon = 1.0;
  105. double    z_mean;
  106. double    z_scale = 1.0;        /* Root mean square range of z after removing planar trend  */
  107. double    r_z_scale = 1.0;    /* reciprocal of z_scale  */
  108. double    plane_c0, plane_c1, plane_c2;    /* Coefficients of best fitting plane to data  */
  109. double small;            /* Let data point coincide with node if distance < small */
  110. float *u;            /* Pointer to grid array */
  111. float *v2;            /* Pointer to v.2.0 grid array */
  112. char *iu;            /* Pointer to grid info array */
  113. char mode_type[2] = {'I','D'};    /* D means include data points when iterating
  114.                  * I means just interpolate from larger grid */
  115.  
  116. int    offset[25][12];        /* Indices of 12 nearby points in 25 cases of edge conditions  */
  117. double        coeff[2][12];    /* Coefficients for 12 nearby points, constrained and unconstrained  */
  118.  
  119. double relax_old, relax_new = 1.4;    /* Coefficients for relaxation factor to speed up convergence */
  120.  
  121. int compare_points();
  122.  
  123. struct DATA {
  124.     float x;
  125.     float y;
  126.     float z;
  127.     int index;
  128. } *data;        /* Data point and index to node it currently constrains  */
  129.  
  130. struct BRIGGS {
  131.     double b[6];
  132. } *briggs;        /* Coefficients in Taylor series for Laplacian(z) a la I. C. Briggs (1974)  */
  133.  
  134. struct SUGGESTION {    /* Used to find top ten list of faster grid dimensions  */
  135.     int    nx;
  136.     int    ny;
  137.     double    factor;    /* Speed up by a factor of factor  */
  138. };
  139.  
  140.  
  141. FILE    *fp_in = NULL;    /* File pointer  */
  142.  
  143. main(argc, argv)
  144. int argc;
  145. char **argv; {
  146.  
  147.     void    suggest_sizes_for_surface();
  148.     int    i, error = FALSE, size_query = FALSE, binary = FALSE, single_precision = FALSE;
  149.     char    modifier, *grdfile, low[100], high[100];
  150.     struct GRD_HEADER h;
  151.     grdfile = 0;
  152.     
  153.     argc = gmt_begin (argc, argv);
  154.  
  155.     grd_init (&h, argc, argv, FALSE);
  156.  
  157.     xmin = ymin = 0.0;
  158.     xmax = ymax = 0.0;
  159.     xinc = yinc = 0.0;
  160.  
  161.     /* New in v4.3:  Default to unconstrained:  */
  162.     set_low = set_high = 0; 
  163.  
  164.     for (i = 1; i < argc; i++) {
  165.         if (argv[i][0] == '-') {
  166.             switch (argv[i][1]) {
  167.               
  168.                 /* Common parameters */
  169.                       
  170.                 case 'V':
  171.                     if (argv[i][2] == 'L' || argv[i][2] == 'l') long_verbose = TRUE;
  172.                 case 'H':
  173.                 case 'R':
  174.                 case ':':
  175.                 case '\0':
  176.                                       error += get_common_args (argv[i], &xmin, &xmax, &ymin, &ymax);
  177.                                       break;
  178.                               
  179.                 /* Supplemental parameters */
  180.                               
  181.                 case 'b':    /* Input triplets are binary, not ascii */
  182.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  183.                     binary = TRUE;
  184.                     break;
  185.                 case 'A':
  186.                     epsilon = atof (&argv[i][2]);
  187.                     break;
  188.                 case 'C':
  189.                     converge_limit = atof (&argv[i][2]);
  190.                     break;
  191.                 case 'G':
  192.                     grdfile = &argv[i][2];
  193.                     break;
  194.                 case 'I':
  195.                     gmt_getinc (&argv[i][2], &xinc, &yinc);
  196.                     break;
  197.                 case 'L':    /* Set limits */
  198.                         /* This is new, to use -Ll and -Lu:  */
  199.                     switch (argv[i][2]) {
  200.                         case 'l':
  201.                             /* Lower limit  */
  202.                             if (argv[i][3] == 0) {
  203.                                 fprintf(stderr,"%s: GMT SYNTAX ERROR -Ll option: No argument given\n", gmt_program);
  204.                                 error++;
  205.                             }
  206.                             strcpy (low, &argv[i][3]);
  207.                             if (!access (low, R_OK))    /* file exists */
  208.                                 set_low = 3;
  209.                             else if (low[0] == 'd')
  210.                                 set_low = 1;
  211.                             else {
  212.                                 set_low = 2;
  213.                                 low_limit = atof (&argv[i][3]);
  214.                             }
  215.                             break;
  216.                         case 'u':
  217.                             /* Upper limit  */
  218.                             if (argv[i][3] == 0) {
  219.                                 fprintf(stderr,"%s: GMT SYNTAX ERROR -Lu option: No argument given\n", gmt_program);
  220.                                 error++;
  221.                             }
  222.                             strcpy (high, &argv[i][3]);
  223.                             if (!access (high, R_OK))    /* file exists */
  224.                                 set_high = 3;
  225.                             else if (high[0] == 'd')
  226.                                 set_high = 1;
  227.                             else {
  228.                                 set_high = 2;
  229.                                 high_limit = atof (&argv[i][3]);
  230.                             }
  231.                             break;
  232.                         default:
  233.                             fprintf(stderr,"%s: GMT SYNTAX ERROR -L option: Unrecongized modifier %c\n", gmt_program, argv[i][2]);
  234.                             error = TRUE;
  235.                             break;
  236.                     }
  237.                     break;
  238.                 case 'N':
  239.                     max_iterations = atoi (&argv[i][2]);
  240.                     break;
  241.                 case 'S':
  242.                     radius = atof (&argv[i][2]);
  243.                     modifier = argv[i][strlen(argv[i])-1];
  244.                     if (modifier == 'm' || modifier == 'M') radius /= 60.0;
  245.                     break;
  246.                 case 'T':
  247.                     modifier = argv[i][strlen(argv[i])-1];
  248.                     if (modifier == 'b' || modifier == 'B') {
  249.                         boundary_tension = atof (&argv[i][2]);
  250.                     }
  251.                     else if (modifier == 'i' || modifier == 'I') {
  252.                         interior_tension = atof (&argv[i][2]);
  253.                     }
  254.                     else if (modifier >= '0' && modifier <= '9') {
  255.                         tension = atof (&argv[i][2]);
  256.                     }
  257.                     else {
  258.                         fprintf(stderr,"%s: GMT SYNTAX ERROR -T option: Unrecongized modifier %c\n", gmt_program, modifier);
  259.                         error = TRUE;
  260.                     }
  261.                     break;
  262.                 case 'Q':
  263.                     size_query = TRUE;
  264.                     break;
  265.                 case 'Z':
  266.                     relax_new = atof (&argv[i][2]);
  267.                     break;
  268.                 default:
  269.                     error = TRUE;
  270.                     gmt_default_error (argv[i][1]);
  271.                     break;
  272.             }
  273.         }
  274.         else {
  275.             if ((fp_in = fopen(argv[i],"r")) == NULL) {
  276.                 fprintf (stderr, "surface: cannot open input data file %s\n", argv[i]);
  277.                 exit(-1);
  278.             }
  279.         }
  280.     }
  281.     
  282.     if (argc == 1 || gmt_quick) {    /* Display usage */
  283.         fprintf (stderr, "surface %s - Adjustable tension continuous curvature surface gridding\n\n", GMT_VERSION);
  284.         fprintf (stderr, "usage: surface [xyz-file] -G<output_grdfile_name> -I<xinc>[m|c][/<yinc>[m|c]]\n");
  285.         fprintf (stderr, "\t-R<west>/<east>/<south>/<north> [-A<aspect_ratio>] [-C<convergence_limit>]\n");
  286.         fprintf (stderr, "\t[-H] [-Ll<limit>] [-Lu<limit>] [-N<n_iterations>] ] [-S<search_radius>[m]]\n");
  287.         fprintf (stderr, "\t[-T<tension>[i][b] ] [-Q] [-V[l]] [-Z<over_relaxation_parameter>] [-:] [-b[d]]\n\n");
  288.         
  289.         if (gmt_quick) exit (-1);
  290.         
  291.         fprintf (stderr, "\tsurface will read from standard input or [xyz-file].\n\n");
  292.         fprintf (stderr, "\tRequired arguments to surface:\n");
  293.         fprintf (stderr, "\t-G sets output grd File name\n");
  294.         fprintf (stderr, "\t-I sets the Increment of the grid; enter xinc, optionally xinc/yinc.\n");
  295.         fprintf (stderr, "\t\tDefault is yinc = xinc.  Append an m [or c] to xinc or yinc to indicate minutes [or seconds]\n");
  296.         fprintf (stderr, "\t\te.g.  -I10m/5m grids longitude every 10 minutes, latitude every 5 minutes.\n\n");
  297.         explain_option ('R');
  298.         fprintf (stderr, "\n\tOPTIONS:\n");
  299.         fprintf (stderr, "\t-A<aspect_ratio>  = 1.0  by default which gives an isotropic solution.\n");
  300.         fprintf (stderr, "\t\ti.e. xinc and yinc assumed to give derivatives of equal weight; if not, specify\n");
  301.         fprintf (stderr, "\t\t<aspect_ratio> such that yinc = xinc / <aspect_ratio>.\n");
  302.         fprintf (stderr, "\t\te.g. if gridding lon,lat use <aspect_ratio> = cosine(middle of lat range).\n");
  303.         fprintf (stderr, "\t-C<convergence_limit> iteration stops when max abs change is less than <c.l.>\n");
  304.         fprintf (stderr, "\t\tdefault will choose 0.001 of the range of your z data (1 ppt precision).\n");
  305.         fprintf (stderr, "\t\tEnter your own convergence limit in same units as z data.\n");
  306.         explain_option ('H');
  307.         fprintf (stderr, "\t-L constrain the range of output values:\n");
  308.         fprintf (stderr, "\t\t-Ll<limit> specifies lower limit; forces solution to be >= <limit>.\n");
  309.         fprintf (stderr, "\t\t-Lu<limit> specifies upper limit; forces solution to be <= <limit>.\n");
  310.         fprintf (stderr, "\t\t<limit> can be any number, or the letter d for min (or max) input data value,\n");
  311.         fprintf (stderr, "\t\tor the filename of a grdfile with bounding values.  [Default solution unconstrained].\n");
  312.         fprintf (stderr, "\t\tExample:  -Ll0 gives a non-negative solution.\n");
  313.         fprintf (stderr, "\t-N sets max <n_iterations> in each cycle; default = 250.\n");
  314.         fprintf (stderr, "\t-S sets <search_radius> to initialize grid; default = 0 will skip this step.\n");
  315.         fprintf (stderr, "\t\tThis step is slow and not needed unless grid dimensions are pathological;\n");
  316.         fprintf (stderr, "\t\ti.e., have few or no common factors.\n");
  317.         fprintf (stderr, "\t\tAppend m to give <search_radius> in minutes.\n");
  318.         fprintf (stderr, "\t-T adds Tension to the gridding equation; use a value between 0 and 1.\n");
  319.         fprintf (stderr, "\t\tdefault = 0 gives minimum curvature (smoothest; bicubic) solution.\n");
  320.         fprintf (stderr, "\t\t1 gives a harmonic spline solution (local max/min occur only at data points).\n");
  321.         fprintf (stderr, "\t\ttypically 0.25 or more is good for potential field (smooth) data;\n");
  322.         fprintf (stderr, "\t\t0.75 or so for topography.  Experiment.\n");
  323.         fprintf (stderr, "\t\tAppend B or b to set tension in boundary conditions only;\n");
  324.         fprintf (stderr, "\t\tAppend I or i to set tension in interior equations only;\n");
  325.         fprintf (stderr, "\t\tNo appended letter sets tension for both to same value.\n");
  326.         fprintf (stderr, "\t-Q Query for grid sizes that might run faster than your -R -I give.\n");
  327.         explain_option ('V');
  328.         fprintf (stderr, "\t\tAppend l for long verbose\n");
  329.         fprintf (stderr, "\t-Z sets <over_relaxation parameter>.  Default = 1.4\n");
  330.         fprintf (stderr, "\t\tUse a value between 1 and 2.  Larger number accelerates convergence but can be unstable.\n");
  331.         fprintf (stderr, "\t\tUse 1 if you want to be sure to have (slow) stable convergence.\n\n");
  332.         explain_option (':');
  333.         fprintf (stderr, "\t-b means input triplets are binary.  Default is ascii.\n");
  334.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");
  335.         explain_option ('.');
  336.         fprintf (stderr, "\t(For additional details, see Smith & Wessel, Geophysics, 55, 293-305, 1990.)\n");
  337.         exit (-1);
  338.     }
  339.  
  340.     if (!project_info.region_supplied) {
  341.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  342.         error++;
  343.     }
  344.     if (xinc <= 0.0 || yinc <= 0.0) {
  345.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", gmt_program);
  346.         error++;
  347.     }
  348.     if (max_iterations < 1) {
  349.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  Max iterations must be nonzero\n", gmt_program);
  350.         error++;
  351.     }
  352.     if (relax_new < 1.0 || relax_new > 2.0) {
  353.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option.  Relaxation value must be 1 <= z <= 2\n", gmt_program);
  354.         error++;
  355.     }
  356.     if (!grdfile) {
  357.         fprintf (stderr, "%s: GMT SYNTAX ERROR option -G:  Must specify output file\n", gmt_program);
  358.         error++;
  359.     }    
  360.     if (binary && gmtdefs.io_header) {
  361.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  362.         error++;
  363.     }
  364.     
  365.     if (error) exit (-1);
  366.  
  367.     if (binary) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  368.     if (tension != 0.0) {
  369.         boundary_tension = tension;
  370.         interior_tension = tension;
  371.     }
  372.     relax_old = 1.0 - relax_new;
  373.  
  374.     nx = rint((xmax - xmin)/xinc) + 1;
  375.     ny = rint((ymax - ymin)/yinc) + 1;
  376.     mx = nx + 4;
  377.     my = ny + 4;
  378.     r_xinc = 1.0 / xinc;
  379.     r_yinc = 1.0 / yinc;
  380.  
  381.     /* New stuff here for v4.3:  Check out the grid dimensions:  */
  382.     grid = gcd_euclid(nx-1, ny-1);
  383.  
  384.     if (gmtdefs.verbose || size_query || grid == 1) fprintf (stderr, "W: %.3lf E: %.3lf S: %.3lf N: %.3lf nx: %d ny: %d\n",
  385.         xmin, xmax, ymin, ymax, nx-1, ny-1);
  386.     if (grid == 1) fprintf(stderr,"surface:  WARNING:  Your grid dimensions are mutually prime.\n");
  387.     if (grid == 1 || size_query) suggest_sizes_for_surface(nx-1, ny-1);
  388.     if (size_query) exit(0);
  389.  
  390.     /* New idea: set grid = 1, read data, setting index.  Then throw
  391.         away data that can't be used in end game, constraining
  392.         size of briggs->b[6] structure.  */
  393.     
  394.     grid = 1;
  395.     set_grid_parameters();
  396.     if (fp_in == NULL) fp_in = stdin;
  397.     read_data();
  398.     throw_away_unusables();
  399.     remove_planar_trend();
  400.     rescale_z_values();
  401.     load_constraints(low, high);
  402.     
  403.     /* Set up factors and reset grid to first value  */
  404.     
  405.     grid = gcd_euclid(nx-1, ny-1);
  406.     n_fact = get_prime_factors(grid, factors);
  407.     set_grid_parameters();
  408.     while ( block_nx < 4 || block_ny < 4 ) {
  409.         smart_divide();
  410.         set_grid_parameters();
  411.     }
  412.     set_offset();
  413.     set_index();
  414.     /* Now the data are ready to go for the first iteration.  */
  415.  
  416.     /* Allocate more space  */
  417.     
  418.     briggs = (struct BRIGGS *) memory (CNULL, npoints, sizeof(struct BRIGGS), "surface");
  419.     iu = (char *) memory (CNULL, mx * my, sizeof(char), "surface");
  420.     u = (float *) memory (CNULL, mx * my, sizeof(float), "surface");
  421.  
  422.     if (radius > 0) initialize_grid(); /* Fill in nodes with a weighted avg in a search radius  */
  423.  
  424.     if (gmtdefs.verbose) fprintf(stderr,"Grid\tMode\tIteration\tMax Change\tConv Limit\tTotal Iterations\n");
  425.     
  426.     set_coefficients();
  427.     
  428.     old_grid = grid;
  429.     find_nearest_point ();
  430.     iterate (1);
  431.      
  432.     while (grid > 1) {
  433.         smart_divide ();
  434.         set_grid_parameters();
  435.         set_offset();
  436.         set_index ();
  437.         fill_in_forecast ();
  438.         iterate(0);
  439.         old_grid = grid;
  440.         find_nearest_point ();
  441.         iterate (1);
  442.     }
  443.     
  444.     if (gmtdefs.verbose) check_errors ();
  445.  
  446.     replace_planar_trend();
  447.     
  448.     free ((char *) data);
  449.     free ((char *) briggs);
  450.     free (iu);
  451.     if (set_low) free ((char *)lower);
  452.     if (set_high) free ((char *)upper);
  453.  
  454.     write_output(&h, grdfile);
  455.  
  456.     free ((char *) u);
  457.  
  458.     
  459.     gmt_end (argc, argv);
  460. }
  461.  
  462. int    set_coefficients()
  463. {
  464.     double    e_4, loose, a0;
  465.     
  466.     loose = 1.0 - interior_tension;
  467.     e_2 = epsilon * epsilon;
  468.     e_4 = e_2 * e_2;
  469.     eps_p2 = e_2;
  470.     eps_m2 = 1.0/e_2;
  471.     one_plus_e2 = 1.0 + e_2;
  472.     two_plus_ep2 = 2.0 + 2.0*eps_p2;
  473.     two_plus_em2 = 2.0 + 2.0*eps_m2;
  474.     
  475.     x_edge_const = 4 * one_plus_e2 - 2 * (interior_tension / loose);
  476.     e_m2 = 1.0 / e_2;
  477.     y_edge_const = 4 * (1.0 + e_m2) - 2 * (interior_tension * e_m2 / loose);
  478.  
  479.     
  480.     a0 = 1.0 / ( (6 * e_4 * loose + 10 * e_2 * loose + 8 * loose - 2 * one_plus_e2) + 4*interior_tension*one_plus_e2);
  481.     a0_const_1 = 2 * loose * (1.0 + e_4);
  482.     a0_const_2 = 2.0 - interior_tension + 2 * loose * e_2;
  483.     
  484.     coeff[1][4] = coeff[1][7] = -loose;
  485.     coeff[1][0] = coeff[1][11] = -loose * e_4;
  486.     coeff[0][4] = coeff[0][7] = -loose * a0;
  487.     coeff[0][0] = coeff[0][11] = -loose * e_4 * a0;
  488.     coeff[1][5] = coeff[1][6] = 2 * loose * one_plus_e2;
  489.     coeff[0][5] = coeff[0][6] = (2 * coeff[1][5] + interior_tension) * a0;
  490.     coeff[1][2] = coeff[1][9] = coeff[1][5] * e_2;
  491.     coeff[0][2] = coeff[0][9] = coeff[0][5] * e_2;
  492.     coeff[1][1] = coeff[1][3] = coeff[1][8] = coeff[1][10] = -2 * loose * e_2;
  493.     coeff[0][1] = coeff[0][3] = coeff[0][8] = coeff[0][10] = coeff[1][1] * a0;
  494.     
  495.     e_2 *= 2;        /* We will need these in boundary conditions  */
  496.     e_m2 *= 2;
  497.     
  498.     ij_sw_corner = 2 * my + 2;            /*  Corners of array of actual data  */
  499.     ij_se_corner = ij_sw_corner + (nx - 1) * my;
  500.     ij_nw_corner = ij_sw_corner + (ny - 1);
  501.     ij_ne_corner = ij_se_corner + (ny - 1);
  502.  
  503. }
  504.  
  505. int    set_offset()
  506. {
  507.     int    add_w[5], add_e[5], add_s[5], add_n[5], add_w2[5], add_e2[5], add_s2[5], add_n2[5];
  508.     int    i, j, kase;
  509.     
  510.     add_w[0] = -my; add_w[1] = add_w[2] = add_w[3] = add_w[4] = -grid_east;
  511.     add_w2[0] = -2 * my;  add_w2[1] = -my - grid_east;  add_w2[2] = add_w2[3] = add_w2[4] = -2 * grid_east;
  512.     add_e[4] = my; add_e[0] = add_e[1] = add_e[2] = add_e[3] = grid_east;
  513.     add_e2[4] = 2 * my;  add_e2[3] = my + grid_east;  add_e2[2] = add_e2[1] = add_e2[0] = 2 * grid_east;
  514.  
  515.     add_n[4] = 1; add_n[3] = add_n[2] = add_n[1] = add_n[0] = grid;
  516.     add_n2[4] = 2;  add_n2[3] = grid + 1;  add_n2[2] = add_n2[1] = add_n2[0] = 2 * grid;
  517.     add_s[0] = -1; add_s[1] = add_s[2] = add_s[3] = add_s[4] = -grid;
  518.     add_s2[0] = -2;  add_s2[1] = -grid - 1;  add_s2[2] = add_s2[3] = add_s2[4] = -2 * grid;
  519.  
  520.     for (i = 0, kase = 0; i < 5; i++) {
  521.         for (j = 0; j < 5; j++, kase++) {
  522.             offset[kase][0] = add_n2[j];
  523.             offset[kase][1] = add_n[j] + add_w[i];
  524.             offset[kase][2] = add_n[j];
  525.             offset[kase][3] = add_n[j] + add_e[i];
  526.             offset[kase][4] = add_w2[i];
  527.             offset[kase][5] = add_w[i];
  528.             offset[kase][6] = add_e[i];
  529.             offset[kase][7] = add_e2[i];
  530.             offset[kase][8] = add_s[j] + add_w[i];
  531.             offset[kase][9] = add_s[j];
  532.             offset[kase][10] = add_s[j] + add_e[i];
  533.             offset[kase][11] = add_s2[j];
  534.         }
  535.     }
  536. }
  537.  
  538.  
  539.  
  540. int fill_in_forecast () {
  541.  
  542.     /* Fills in bilinear estimates into new node locations
  543.        after grid is divided.   
  544.      */
  545.  
  546.     int i, j, ii, jj, index_0, index_1, index_2, index_3;
  547.     int index_new;
  548.     double delta_x, delta_y, a0, a1, a2, a3;
  549.     double old_size;
  550.     
  551.         
  552.     old_size = 1.0 / (double)old_grid;
  553.  
  554.     /* first do from southwest corner */
  555.     
  556.     for (i = 0; i < nx-1; i += old_grid) {
  557.         
  558.         for (j = 0; j < ny-1; j += old_grid) {
  559.             
  560.             /* get indices of bilinear square */
  561.             index_0 = ij_sw_corner + i * my + j;
  562.             index_1 = index_0 + old_grid * my;
  563.             index_2 = index_1 + old_grid;
  564.             index_3 = index_0 + old_grid;
  565.             
  566.             /* get coefficients */
  567.             a0 = u[index_0];
  568.             a1 = u[index_1] - a0;
  569.             a2 = u[index_3] - a0;
  570.             a3 = u[index_2] - a0 - a1 - a2;
  571.             
  572.             /* find all possible new fill ins */
  573.             
  574.             for (ii = i;  ii < i + old_grid; ii += grid) {
  575.                 delta_x = (ii - i) * old_size;
  576.                 for (jj = j;  jj < j + old_grid; jj += grid) {
  577.                     index_new = ij_sw_corner + ii * my + jj;
  578.                     if (index_new == index_0) continue;
  579.                     delta_y = (jj - j) * old_size;
  580.                     u[index_new] = a0 + a1 * delta_x + delta_y * ( a2 + a3 * delta_x);    
  581.                     iu[index_new] = 0;
  582.                 }
  583.             }
  584.             iu[index_0] = 5;
  585.         }
  586.     }
  587.     
  588.     /* now do linear guess along east edge */
  589.     
  590.     for (j = 0; j < (ny-1); j += old_grid) {
  591.         index_0 = ij_se_corner + j;
  592.         index_3 = index_0 + old_grid;
  593.         for (jj = j;  jj < j + old_grid; jj += grid) {
  594.             index_new = ij_se_corner + jj;
  595.             delta_y = (jj - j) * old_size;
  596.             u[index_new] = u[index_0] + delta_y * (u[index_3] - u[index_0]);
  597.             iu[index_new] = 0;
  598.         }
  599.         iu[index_0] = 5;
  600.     }
  601.     /* now do linear guess along north edge */
  602.     for (i = 0; i < (nx-1); i += old_grid) {
  603.         index_0 = ij_nw_corner + i * my;
  604.         index_1 = index_0 + old_grid * my;
  605.         for (ii = i;  ii < i + old_grid; ii += grid) {
  606.             index_new = ij_nw_corner + ii * my;
  607.             delta_x = (ii - i) * old_size;
  608.             u[index_new] = u[index_0] + delta_x * (u[index_1] - u[index_0]);
  609.             iu[index_new] = 0;
  610.         }
  611.         iu[index_0] = 5;
  612.     }
  613.     /* now set northeast corner to fixed and we're done */
  614.     iu[ij_ne_corner] = 5;
  615. }
  616.  
  617. int compare_points (point_1, point_2)
  618. struct DATA *point_1, *point_2; {
  619.         /*  Routine for qsort to sort data structure for fast access to data by node location.
  620.             Sorts on index first, then on radius to node corresponding to index, so that index
  621.             goes from low to high, and so does radius.
  622.         */
  623.     int block_i, block_j, index_1, index_2;
  624.     double x0, y0, dist_1, dist_2;
  625.     
  626.     index_1 = point_1->index;
  627.     index_2 = point_2->index;
  628.     if (index_1 < index_2)
  629.         return (-1);
  630.     else if (index_1 > index_2)
  631.         return (1);
  632.     else if (index_1 == OUTSIDE)
  633.         return (0);
  634.     else {    /* Points are in same grid cell, find the one who is nearest to grid point */
  635.         block_i = point_1->index/block_ny;
  636.         block_j = point_1->index%block_ny;
  637.         x0 = xmin + block_i * grid_xinc;
  638.         y0 = ymin + block_j * grid_yinc;
  639.         dist_1 = (point_1->x - x0) * (point_1->x - x0) + (point_1->y - y0) * (point_1->y - y0);
  640.         dist_2 = (point_2->x - x0) * (point_2->x - x0) + (point_2->y - y0) * (point_2->y - y0);
  641.         if (dist_1 < dist_2)
  642.             return (-1);
  643.         if (dist_1 > dist_2)
  644.             return (1);
  645.         else
  646.             return (0);
  647.     }
  648. }
  649.  
  650. int smart_divide () {
  651.         /* Divide grid by its largest prime factor */
  652.     grid /= factors[n_fact - 1];
  653.     n_fact--;
  654. }
  655.  
  656. int set_index () {
  657.         /* recomputes data[k].index for new value of grid,
  658.            sorts data on index and radii, and throws away
  659.            data which are now outside the useable limits. */
  660.     int i, j, k, k_skipped = 0;
  661.  
  662.     for (k = 0; k < npoints; k++) {
  663.         i = floor(((data[k].x-xmin)*r_grid_xinc) + 0.5);
  664.         j = floor(((data[k].y-ymin)*r_grid_yinc) + 0.5);
  665.         if (i < 0 || i >= block_nx || j < 0 || j >= block_ny) {
  666.             data[k].index = OUTSIDE;
  667.             k_skipped++;
  668.         }
  669.         else
  670.             data[k].index = i * block_ny + j;
  671.     }
  672.     
  673.     qsort ((char *)data, npoints, sizeof (struct DATA), compare_points);
  674.     
  675.     npoints -= k_skipped;
  676.     
  677. }
  678.  
  679. int find_nearest_point() {
  680.     int i, j, k, last_index, block_i, block_j, iu_index, briggs_index;
  681.     double x0, y0, dx, dy, xys, xy1, btemp;
  682.     double b0, b1, b2, b3, b4, b5;
  683.     
  684.     last_index = -1;
  685.     small = 0.05 * ((grid_xinc < grid_yinc) ? grid_xinc : grid_yinc);
  686.  
  687.     for (i = 0; i < nx; i += grid)    /* Reset grid info */
  688.         for (j = 0; j < ny; j += grid)
  689.             iu[ij_sw_corner + i*my + j] = 0;
  690.     
  691.     briggs_index = 0;
  692.     for (k = 0; k < npoints; k++) {    /* Find constraining value  */
  693.         if (data[k].index != last_index) {
  694.             block_i = data[k].index/block_ny;
  695.             block_j = data[k].index%block_ny;
  696.             last_index = data[k].index;
  697.              iu_index = ij_sw_corner + (block_i * my + block_j) * grid;
  698.              x0 = xmin + block_i*grid_xinc;
  699.              y0 = ymin + block_j*grid_yinc;
  700.              dx = (data[k].x - x0)*r_grid_xinc;
  701.              dy = (data[k].y - y0)*r_grid_yinc;
  702.              if (fabs(dx) < small && fabs(dy) < small) {
  703.                  iu[iu_index] = 5;
  704.                  u[iu_index] = data[k].z;
  705.              }
  706.              else {
  707.                  if (dx >= 0.0) {
  708.                      if (dy >= 0.0)
  709.                          iu[iu_index] = 1;
  710.                      else
  711.                          iu[iu_index] = 4;
  712.                  }
  713.                  else {
  714.                      if (dy >= 0.0)
  715.                          iu[iu_index] = 2;
  716.                      else
  717.                          iu[iu_index] = 3;
  718.                  }
  719.                  dx = fabs(dx);
  720.                  dy = fabs(dy);
  721.                  btemp = 2 * one_plus_e2 / ( (dx + dy) * (1.0 + dx + dy) );
  722.                  b0 = 1.0 - 0.5 * (dx + (dx * dx)) * btemp;
  723.                  b3 = 0.5 * (e_2 - (dy + (dy * dy)) * btemp);
  724.                  xys = 1.0 + dx + dy;
  725.                  xy1 = 1.0 / xys;
  726.                  b1 = (e_2 * xys - 4 * dy) * xy1;
  727.                  b2 = 2 * (dy - dx + 1.0) * xy1;
  728.                  b4 = b0 + b1 + b2 + b3 + btemp;
  729.                  b5 = btemp * data[k].z;
  730.                  briggs[briggs_index].b[0] = b0;
  731.                  briggs[briggs_index].b[1] = b1;
  732.                  briggs[briggs_index].b[2] = b2;
  733.                  briggs[briggs_index].b[3] = b3;
  734.                  briggs[briggs_index].b[4] = b4;
  735.                  briggs[briggs_index].b[5] = b5;
  736.                  briggs_index++;
  737.              }
  738.          }
  739.      }
  740. }
  741.  
  742.                         
  743. int set_grid_parameters()
  744. {            
  745.     block_ny = (ny - 1) / grid + 1;
  746.     block_nx = (nx - 1) / grid + 1;
  747.     grid_xinc = grid * xinc;
  748.     grid_yinc = grid * yinc;
  749.     grid_east = grid * my;
  750.     r_grid_xinc = 1.0 / grid_xinc;
  751.     r_grid_yinc = 1.0 / grid_yinc;
  752. }
  753.  
  754. int initialize_grid()
  755. {    /*
  756.      * For the initial gridsize, compute weighted averages of data inside the search radius
  757.      * and assign the values to u[i,j] where i,j are multiples of gridsize.
  758.      */
  759.      int    irad, jrad, i, j, imin, imax, jmin, jmax, index_1, index_2, k, ki, kj, k_index;
  760.      double    r, rfact, sum_w, sum_zw, weight, x0, y0;
  761.  
  762.      irad = ceil(radius/grid_xinc);
  763.      jrad = ceil(radius/grid_yinc);
  764.      rfact = -4.5/(radius*radius);
  765.      
  766.      for (i = 0; i < block_nx; i ++ ) {
  767.          x0 = xmin + i*grid_xinc;
  768.          for (j = 0; j < block_ny; j ++ ) {
  769.              y0 = ymin + j*grid_yinc;
  770.              imin = i - irad;
  771.              if (imin < 0) imin = 0;
  772.              imax = i + irad;
  773.              if (imax >= block_nx) imax = block_nx - 1;
  774.              jmin = j - jrad;
  775.              if (jmin < 0) jmin = 0;
  776.              jmax = j + jrad;
  777.              if (jmax >= block_ny) jmax = block_ny - 1;
  778.              index_1 = imin*block_ny + jmin;
  779.              index_2 = imax*block_ny + jmax + 1;
  780.              sum_w = sum_zw = 0.0;
  781.              k = 0;
  782.              while (k < npoints && data[k].index < index_1) k++;
  783.              for (ki = imin; k < npoints && ki <= imax && data[k].index < index_2; ki++) {
  784.                  for (kj = jmin; k < npoints && kj <= jmax && data[k].index < index_2; kj++) {
  785.                      k_index = ki*block_ny + kj;
  786.                      while (k < npoints && data[k].index < k_index) k++;
  787.                      while (k < npoints && data[k].index == k_index) {
  788.                          r = (data[k].x-x0)*(data[k].x-x0) + (data[k].y-y0)*(data[k].y-y0);
  789.                          weight = exp (rfact*r);
  790.                          sum_w += weight;
  791.                          sum_zw += weight*data[k].z;
  792.                          k++;
  793.                      }
  794.                  }
  795.              }
  796.              if (sum_w == 0.0) {
  797.                  fprintf (stderr, "surface: Warning: no data inside search radius at: %.8lg %.8lg\n", x0, y0);
  798.                  u[ij_sw_corner + (i * my + j) * grid] = z_mean;
  799.              }
  800.              else {
  801.                  u[ij_sw_corner + (i*my+j)*grid] = sum_zw/sum_w;
  802.              }
  803.         }
  804.     }
  805. }
  806.  
  807.  
  808. int new_initialize_grid()
  809. {    /*
  810.      * For the initial gridsize, load constrained nodes with weighted avg of their data;
  811.      * and then do something with the unconstrained ones.
  812.      */
  813.      int    k, k_index, u_index, block_i, block_j;
  814.      double    sum_w, sum_zw, weight, x0, y0, dx, dy, dx_scale, dy_scale;
  815.  
  816.     dx_scale = 4.0 / grid_xinc;
  817.     dy_scale = 4.0 / grid_yinc;
  818.     n_empty = block_ny * block_nx;
  819.     k = 0;
  820.     while (k < npoints) {
  821.         block_i = data[k].index / block_ny;
  822.         block_j = data[k].index % block_ny;
  823.         x0 = xmin + block_i*grid_xinc;
  824.         y0 = ymin + block_j*grid_yinc;
  825.         u_index = ij_sw_corner + (block_i*my + block_j) * grid;
  826.         k_index = data[k].index;
  827.         
  828.         dy = (data[k].y - y0) * dy_scale;
  829.         dx = (data[k].x - x0) * dx_scale;
  830.         sum_w = 1.0 / (1.0 + dx*dx + dy*dy);
  831.         sum_zw = data[k].z * sum_w;
  832.         k++;
  833.  
  834.         while (k < npoints && data[k].index == k_index) {
  835.             
  836.             dy = (data[k].y - y0) * dy_scale;
  837.             dx = (data[k].x - x0) * dx_scale;
  838.             weight = 1.0 / (1.0 + dx*dx + dy*dy);
  839.             sum_zw += data[k].z * weight;
  840.             sum_w += weight;
  841.             sum_zw += weight*data[k].z;
  842.             k++;
  843.          }
  844.          u[u_index] = sum_zw/sum_w;
  845.          iu[u_index] = 5;
  846.          n_empty--;
  847.      }
  848. }
  849.  
  850. int read_data()
  851. {
  852.     int    i, j, k, kmax, kmin, ix, iy, n_fields, more, skip;
  853.     double    *in, zmin = 1.0e38, zmax = -1.0e38;
  854.     char    buffer[512];
  855.  
  856.     data = (struct DATA *) memory (CNULL, n_alloc, sizeof(struct DATA), "surface");
  857.     
  858.     /* Read in xyz data and computes index no and store it in a structure */
  859.     
  860.     ix = (gmtdefs.xy_toggle) ? 1 : 0;        iy = 1 - ix;              /* Set up which columns have x and y */
  861.     k = 0;
  862.     z_mean = 0;
  863.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (buffer, 512, fp_in);
  864.     
  865.     more = ((n_fields = gmt_input (fp_in,  3, &in)) == 3);
  866.         
  867.     while (more) {
  868.         skip = FALSE;
  869.         
  870.         if (bad_float ((double)in[2])) skip = TRUE;
  871.         
  872.         if (!skip) {
  873.             i = floor(((in[ix]-xmin)*r_grid_xinc) + 0.5);
  874.             if (i < 0 || i >= block_nx) skip = TRUE;
  875.             j = floor(((in[iy]-ymin)*r_grid_yinc) + 0.5);
  876.             if (j < 0 || j >= block_ny) skip = TRUE;
  877.         }
  878.         if (!skip) {
  879.             data[k].index = i * block_ny + j;
  880.             data[k].x = in[ix];
  881.             data[k].y = in[iy];
  882.             data[k].z = in[2];
  883.             if (zmin > in[2]) {
  884.                 zmin = in[2];
  885.                 kmin = k;
  886.             }
  887.             if (zmax < in[2]) {
  888.                 zmax = in[2];
  889.                 kmax = k;
  890.             }
  891.             k++;
  892.             z_mean += in[2];
  893.             if (k == n_alloc) {
  894.                 n_alloc += GMT_CHUNK;
  895.                 data = (struct DATA *) memory ((char *)data, n_alloc, sizeof(struct DATA), "surface");
  896.             }
  897.         }
  898.         
  899.         more = ((n_fields = gmt_input (fp_in,  3, &in)) == 3);
  900.     }
  901.     
  902.     if (fp_in != stdin) fclose (fp_in);
  903.     if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  904.     if (n_fields%3) {    /* Got garbage or multiple segment subheader */
  905.         fprintf (stderr, "%s:  Cannot read X Y Z at line %d, aborts\n", gmt_program, k);
  906.         exit (-1);
  907.     }
  908.  
  909.     npoints = k;
  910.     
  911.     if (npoints == 0) {
  912.         fprintf (stderr, "%s:  No datapoints inside region, aborts\n", gmt_program);
  913.         exit (-1);
  914.     }
  915.     
  916.     z_mean /= k;
  917.     if( converge_limit == 0.0 ) {
  918.         converge_limit = 0.001 * z_scale; /* c_l = 1 ppt of L2 scale */
  919.     }
  920.     if (gmtdefs.verbose) {
  921.         fprintf(stderr, "surface: Minimum value of your dataset x,y,z at: %g %g %g\n",
  922.             data[kmin].x, data[kmin].y, data[kmin].z);
  923.         fprintf(stderr, "surface: Maximum value of your dataset x,y,z at: %g %g %g\n",
  924.             data[kmax].x, data[kmax].y, data[kmax].z);
  925.     }
  926.     data = (struct DATA *) memory ((char *)data, npoints, sizeof(struct DATA), "surface");
  927.     
  928.     if (set_low == 1)
  929.         low_limit = data[kmin].z;
  930.     else if (set_low == 2 && low_limit > data[kmin].z) {
  931.     /*    low_limit = data[kmin].z;    */
  932.         fprintf (stderr, "surface: Warning:  Your lower value is > than min data value.\n");
  933.     }
  934.     if (set_high == 1)
  935.         high_limit = data[kmax].z;
  936.     else if (set_high == 2 && high_limit < data[kmax].z) {
  937.     /*    high_limit = data[kmax].z;    */
  938.         fprintf (stderr, "surface: Warning:  Your upper value is < than max data value.\n");
  939.     }
  940.  
  941. }
  942.  
  943. int write_output(h, grdfile)
  944. struct GRD_HEADER *h;
  945. char *grdfile;
  946. {    /* Uses v.2.0 netCDF grd format */
  947.     int    index, i, j, dummy[4];
  948.     
  949.     dummy[3] = dummy[2] = dummy[1] = dummy[0] = 0;
  950.     
  951.     h->x_min = xmin;
  952.     h->x_max = xmax;
  953.     h->y_min = ymin;
  954.     h->y_max = ymax;
  955.     h->nx = nx;
  956.     h->ny = ny;
  957.     h->x_inc = xinc;
  958.     h->y_inc = yinc;
  959.     strcpy (h->title, "surface");
  960.     
  961.     v2 = (float *) memory (CNULL, nx * ny, sizeof (float), "surface");
  962.     index = ij_sw_corner;
  963.     for(i = 0; i < nx; i++, index += my) for (j = 0; j < ny; j++) v2[j*nx+i] = u[index + ny - j - 1];
  964.     if (write_grd (grdfile, h, v2, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  965.         fprintf (stderr, "surface: Error writing file %s\n", grdfile);
  966.         exit (-1);
  967.     }
  968.     free ((char *)v2);
  969. }
  970.     
  971. int    iterate(mode)
  972.     int    mode;
  973. {
  974.  
  975.     int    i, j, k, ij, kase, briggs_index, ij_v2;
  976.     int    x_case, y_case, x_w_case, x_e_case, y_s_case, y_n_case;
  977.     int    iteration_count = 0;
  978.     
  979.     double    current_limit = converge_limit / grid;
  980.     double    change, max_change = 0.0, busum, sum_ij;
  981.     double    b0, b1, b2, b3, b4, b5;
  982.     
  983.     double    x_0_const = 4.0 * (1.0 - boundary_tension) / (2.0 - boundary_tension);
  984.     double    x_1_const = (3 * boundary_tension - 2.0) / (2.0 - boundary_tension);
  985.     double    y_denom = 2 * epsilon * (1.0 - boundary_tension) + boundary_tension;
  986.     double    y_0_const = 4 * epsilon * (1.0 - boundary_tension) / y_denom;
  987.     double    y_1_const = (boundary_tension - 2 * epsilon * (1.0 - boundary_tension) ) / y_denom;
  988.  
  989.     do {
  990.         briggs_index = 0;    /* Reset the constraint table stack pointer  */
  991.         
  992.         max_change = -1.0;
  993.         
  994.         /* Fill in auxiliary boundary values (in new way) */
  995.         
  996.         /* First set d2[]/dn2 = 0 along edges:  */
  997.         /* New experiment : (1-T)d2[]/dn2 + Td[]/dn = 0  */
  998.         
  999.         
  1000.         
  1001.         for (i = 0; i < nx; i += grid) {
  1002.             /* set d2[]/dy2 = 0 on south side:  */
  1003.             ij = ij_sw_corner + i * my;
  1004.             /* u[ij - 1] = 2 * u[ij] - u[ij + grid];  */
  1005.             u[ij - 1] = y_0_const * u[ij] + y_1_const * u[ij + grid];
  1006.             /* set d2[]/dy2 = 0 on north side:  */
  1007.             ij = ij_nw_corner + i * my;
  1008.             /* u[ij + 1] = 2 * u[ij] - u[ij - grid];  */
  1009.             u[ij + 1] = y_0_const * u[ij] + y_1_const * u[ij - grid];
  1010.             
  1011.         }
  1012.         
  1013.         for (j = 0; j < ny; j += grid) {
  1014.             /* set d2[]/dx2 = 0 on west side:  */
  1015.             ij = ij_sw_corner + j;
  1016.             /* u[ij - my] = 2 * u[ij] - u[ij + grid_east];  */
  1017.             u[ij - my] = x_1_const * u[ij + grid_east] + x_0_const * u[ij];
  1018.             /* set d2[]/dx2 = 0 on east side:  */
  1019.             ij = ij_se_corner + j;
  1020.             /* u[ij + my] = 2 * u[ij] - u[ij - grid_east];  */
  1021.             u[ij + my] = x_1_const * u[ij - grid_east] + x_0_const * u[ij];
  1022.         }
  1023.             
  1024.         /* Now set d2[]/dxdy = 0 at each corner:  */
  1025.         
  1026.         ij = ij_sw_corner;
  1027.         u[ij - my - 1] = u[ij + grid_east - 1] + u[ij - my + grid] - u[ij + grid_east + grid];
  1028.                 
  1029.         ij = ij_nw_corner;
  1030.         u[ij - my + 1] = u[ij + grid_east + 1] + u[ij - my - grid] - u[ij + grid_east - grid];
  1031.                 
  1032.         ij = ij_se_corner;
  1033.         u[ij + my - 1] = u[ij - grid_east - 1] + u[ij + my + grid] - u[ij - grid_east + grid];
  1034.                 
  1035.         ij = ij_ne_corner;
  1036.         u[ij + my + 1] = u[ij - grid_east + 1] + u[ij + my - grid] - u[ij - grid_east - grid];
  1037.         
  1038.         /* Now set (1-T)dC/dn + Tdu/dn = 0 at each edge :  */
  1039.         /* New experiment:  only dC/dn = 0  */
  1040.         
  1041.         x_w_case = 0;
  1042.         x_e_case = block_nx - 1;
  1043.         for (i = 0; i < nx; i += grid, x_w_case++, x_e_case--) {
  1044.         
  1045.             if(x_w_case < 2)
  1046.                 x_case = x_w_case;
  1047.             else if(x_e_case < 2)
  1048.                 x_case = 4 - x_e_case;
  1049.             else
  1050.                 x_case = 2;
  1051.                 
  1052.             /* South side :  */
  1053.             kase = x_case * 5;
  1054.             ij = ij_sw_corner + i * my;
  1055.             u[ij + offset[kase][11]] = 
  1056.                 (u[ij + offset[kase][0]] + eps_m2*(u[ij + offset[kase][1]] + u[ij + offset[kase][3]]
  1057.                     - u[ij + offset[kase][8]] - u[ij + offset[kase][10]])
  1058.                     + two_plus_em2 * (u[ij + offset[kase][9]] - u[ij + offset[kase][2]]) );
  1059.                 /*  + tense * eps_m2 * (u[ij + offset[kase][2]] - u[ij + offset[kase][9]]) / (1.0 - tense);  */
  1060.             /* North side :  */
  1061.             kase = x_case * 5 + 4;
  1062.             ij = ij_nw_corner + i * my;
  1063.             u[ij + offset[kase][0]] = 
  1064.                 -(-u[ij + offset[kase][11]] + eps_m2 * (u[ij + offset[kase][1]] + u[ij + offset[kase][3]]
  1065.                     - u[ij + offset[kase][8]] - u[ij + offset[kase][10]])
  1066.                     + two_plus_em2 * (u[ij + offset[kase][9]] - u[ij + offset[kase][2]]) );
  1067.                 /*  - tense * eps_m2 * (u[ij + offset[kase][2]] - u[ij + offset[kase][9]]) / (1.0 - tense);  */
  1068.         }
  1069.         
  1070.         y_s_case = 0;
  1071.         y_n_case = block_ny - 1;
  1072.         for (j = 0; j < ny; j += grid, y_s_case++, y_n_case--) {
  1073.                 
  1074.             if(y_s_case < 2)
  1075.                 y_case = y_s_case;
  1076.             else if(y_n_case < 2)
  1077.                 y_case = 4 - y_n_case;
  1078.             else
  1079.                 y_case = 2;
  1080.             
  1081.             /* West side :  */
  1082.             kase = y_case;
  1083.             ij = ij_sw_corner + j;
  1084.             u[ij+offset[kase][4]] = 
  1085.                 u[ij + offset[kase][7]] + eps_p2 * (u[ij + offset[kase][3]] + u[ij + offset[kase][10]]
  1086.                 -u[ij + offset[kase][1]] - u[ij + offset[kase][8]])
  1087.                 + two_plus_ep2 * (u[ij + offset[kase][5]] - u[ij + offset[kase][6]]);
  1088.                 /*  + tense * (u[ij + offset[kase][6]] - u[ij + offset[kase][5]]) / (1.0 - tense);  */
  1089.             /* East side :  */
  1090.             kase = 20 + y_case;
  1091.             ij = ij_se_corner + j;
  1092.             u[ij + offset[kase][7]] = 
  1093.                 - (-u[ij + offset[kase][4]] + eps_p2 * (u[ij + offset[kase][3]] + u[ij + offset[kase][10]]
  1094.                 - u[ij + offset[kase][1]] - u[ij + offset[kase][8]])
  1095.                 + two_plus_ep2 * (u[ij + offset[kase][5]] - u[ij + offset[kase][6]]) );
  1096.                 /*  - tense * (u[ij + offset[kase][6]] - u[ij + offset[kase][5]]) / (1.0 - tense);  */
  1097.         }
  1098.  
  1099.             
  1100.             
  1101.         /* That's it for the boundary points.  Now loop over all data  */
  1102.         
  1103.         x_w_case = 0;
  1104.         x_e_case = block_nx - 1;
  1105.         for (i = 0; i < nx; i += grid, x_w_case++, x_e_case--) {
  1106.         
  1107.             if(x_w_case < 2)
  1108.                 x_case = x_w_case;
  1109.             else if(x_e_case < 2)
  1110.                 x_case = 4 - x_e_case;
  1111.             else
  1112.                 x_case = 2;
  1113.             
  1114.             y_s_case = 0;
  1115.             y_n_case = block_ny - 1;
  1116.             
  1117.             ij = ij_sw_corner + i * my;
  1118.             
  1119.             for (j = 0; j < ny; j += grid, ij += grid, y_s_case++, y_n_case--) {
  1120.     
  1121.                 if (iu[ij] == 5) continue;    /* Point is fixed  */
  1122.                 
  1123.                 if(y_s_case < 2)
  1124.                     y_case = y_s_case;
  1125.                 else if(y_n_case < 2)
  1126.                     y_case = 4 - y_n_case;
  1127.                 else
  1128.                     y_case = 2;
  1129.                 
  1130.                 kase = x_case * 5 + y_case;
  1131.                 sum_ij = 0.0;
  1132.  
  1133.                 if (iu[ij] == 0) {        /* Point is unconstrained  */
  1134.                     for (k = 0; k < 12; k++) {
  1135.                         sum_ij += (u[ij + offset[kase][k]] * coeff[0][k]);
  1136.                     }
  1137.                 }
  1138.                 else {                /* Point is constrained  */
  1139.                 
  1140.                     b0 = briggs[briggs_index].b[0];
  1141.                     b1 = briggs[briggs_index].b[1];
  1142.                     b2 = briggs[briggs_index].b[2];
  1143.                     b3 = briggs[briggs_index].b[3];
  1144.                     b4 = briggs[briggs_index].b[4];
  1145.                     b5 = briggs[briggs_index].b[5];
  1146.                     briggs_index++;
  1147.                     if (iu[ij] < 3) {
  1148.                         if (iu[ij] == 1) {    /* Point is in quadrant 1  */
  1149.                             busum = b0 * u[ij + offset[kase][10]]
  1150.                                 + b1 * u[ij + offset[kase][9]]
  1151.                                 + b2 * u[ij + offset[kase][5]]
  1152.                                 + b3 * u[ij + offset[kase][1]];
  1153.                         }
  1154.                         else {            /* Point is in quadrant 2  */
  1155.                             busum = b0 * u[ij + offset[kase][8]]
  1156.                                 + b1 * u[ij + offset[kase][9]]
  1157.                                 + b2 * u[ij + offset[kase][6]]
  1158.                                 + b3 * u[ij + offset[kase][3]];
  1159.                         }
  1160.                     }
  1161.                     else {
  1162.                         if (iu[ij] == 3) {    /* Point is in quadrant 3  */
  1163.                             busum = b0 * u[ij + offset[kase][1]]
  1164.                                 + b1 * u[ij + offset[kase][2]]
  1165.                                 + b2 * u[ij + offset[kase][6]]
  1166.                                 + b3 * u[ij + offset[kase][10]];
  1167.                         }
  1168.                         else {        /* Point is in quadrant 4  */
  1169.                             busum = b0 * u[ij + offset[kase][3]]
  1170.                                 + b1 * u[ij + offset[kase][2]]
  1171.                                 + b2 * u[ij + offset[kase][5]]
  1172.                                 + b3 * u[ij + offset[kase][8]];
  1173.                         }
  1174.                     }
  1175.                     for (k = 0; k < 12; k++) {
  1176.                         sum_ij += (u[ij + offset[kase][k]] * coeff[1][k]);
  1177.                     }
  1178.                     sum_ij = (sum_ij + a0_const_2 * (busum + b5))
  1179.                         / (a0_const_1 + a0_const_2 * b4);
  1180.                 }
  1181.                 
  1182.                 /* New relaxation here  */
  1183.                 sum_ij = u[ij] * relax_old + sum_ij * relax_new;
  1184.                 
  1185.                 if (constrained) {    /* Must check limits.  Note lower/upper is v2 format and need ij_v2! */
  1186.                     ij_v2 = (ny - j - 1) * nx + i;
  1187.                     if (set_low && !bad_float((double)lower[ij_v2]) && sum_ij < lower[ij_v2])
  1188.                         sum_ij = lower[ij_v2];
  1189.                     else if (set_high && !bad_float((double)upper[ij_v2]) && sum_ij > upper[ij_v2])
  1190.                         sum_ij = upper[ij_v2];
  1191.                 }
  1192.                     
  1193.                 change = fabs(sum_ij - u[ij]);
  1194.                 u[ij] = sum_ij;
  1195.                 if (change > max_change) max_change = change;
  1196.             }
  1197.         }
  1198.         iteration_count++;
  1199.         total_iterations++;
  1200.         max_change *= z_scale;    /* Put max_change into z units  */
  1201.         if (long_verbose) fprintf(stderr,"%4d\t%c\t%8d\t%10lg\t%10lg\t%10d\n",
  1202.             grid, mode_type[mode], iteration_count, max_change, current_limit, total_iterations);
  1203.  
  1204.     } while (max_change > current_limit && iteration_count < max_iterations);
  1205.     
  1206.     if (gmtdefs.verbose && !long_verbose) fprintf(stderr,"%4d\t%c\t%8d\t%10lg\t%10lg\t%10d\n",
  1207.         grid, mode_type[mode], iteration_count, max_change, current_limit, total_iterations);
  1208.  
  1209.     return(iteration_count);
  1210. }
  1211.  
  1212. int check_errors () {
  1213.  
  1214.     int    i, j, k, ij, n_nodes, move_over[12];    /* move_over = offset[kase][12], but grid = 1 so move_over is easy  */
  1215.     
  1216.     double    x0, y0, dx, dy, mean_error, mean_squared_error, z_est, z_err, curvature, c;
  1217.     double    du_dx, du_dy, d2u_dx2, d2u_dxdy, d2u_dy2, d3u_dx3, d3u_dx2dy, d3u_dxdy2, d3u_dy3;
  1218.     
  1219.     double    x_0_const = 4.0 * (1.0 - boundary_tension) / (2.0 - boundary_tension);
  1220.     double    x_1_const = (3 * boundary_tension - 2.0) / (2.0 - boundary_tension);
  1221.     double    y_denom = 2 * epsilon * (1.0 - boundary_tension) + boundary_tension;
  1222.     double    y_0_const = 4 * epsilon * (1.0 - boundary_tension) / y_denom;
  1223.     double    y_1_const = (boundary_tension - 2 * epsilon * (1.0 - boundary_tension) ) / y_denom;
  1224.     
  1225.     
  1226.     move_over[0] = 2;
  1227.     move_over[1] = 1 - my;
  1228.     move_over[2] = 1;
  1229.     move_over[3] = 1 + my;
  1230.     move_over[4] = -2 * my;
  1231.     move_over[5] = -my;
  1232.     move_over[6] = my;
  1233.     move_over[7] = 2 * my;
  1234.     move_over[8] = -1 - my;
  1235.     move_over[9] = -1;
  1236.     move_over[10] = -1 + my;
  1237.     move_over[11] = -2;
  1238.  
  1239.     mean_error = 0;
  1240.     mean_squared_error = 0;
  1241.     
  1242.     /* First update the boundary values  */
  1243.  
  1244.     for (i = 0; i < nx; i ++) {
  1245.         ij = ij_sw_corner + i * my;
  1246.         u[ij - 1] = y_0_const * u[ij] + y_1_const * u[ij + 1];
  1247.         ij = ij_nw_corner + i * my;
  1248.         u[ij + 1] = y_0_const * u[ij] + y_1_const * u[ij - 1];
  1249.     }
  1250.  
  1251.     for (j = 0; j < ny; j ++) {
  1252.         ij = ij_sw_corner + j;
  1253.         u[ij - my] = x_1_const * u[ij + my] + x_0_const * u[ij];
  1254.         ij = ij_se_corner + j;
  1255.         u[ij + my] = x_1_const * u[ij - my] + x_0_const * u[ij];
  1256.     }
  1257.  
  1258.     ij = ij_sw_corner;
  1259.     u[ij - my - 1] = u[ij + my - 1] + u[ij - my + 1] - u[ij + my + 1];
  1260.     ij = ij_nw_corner;
  1261.     u[ij - my + 1] = u[ij + my + 1] + u[ij - my - 1] - u[ij + my - 1];
  1262.     ij = ij_se_corner;
  1263.     u[ij + my - 1] = u[ij - my - 1] + u[ij + my + 1] - u[ij - my + 1];
  1264.     ij = ij_ne_corner;
  1265.     u[ij + my + 1] = u[ij - my + 1] + u[ij + my - 1] - u[ij - my - 1];
  1266.  
  1267.     for (i = 0; i < nx; i ++) {
  1268.                 
  1269.         ij = ij_sw_corner + i * my;
  1270.         u[ij + move_over[11]] = 
  1271.             (u[ij + move_over[0]] + eps_m2*(u[ij + move_over[1]] + u[ij + move_over[3]]
  1272.                 - u[ij + move_over[8]] - u[ij + move_over[10]])
  1273.                 + two_plus_em2 * (u[ij + move_over[9]] - u[ij + move_over[2]]) );
  1274.                     
  1275.         ij = ij_nw_corner + i * my;
  1276.         u[ij + move_over[0]] = 
  1277.             -(-u[ij + move_over[11]] + eps_m2 * (u[ij + move_over[1]] + u[ij + move_over[3]]
  1278.                 - u[ij + move_over[8]] - u[ij + move_over[10]])
  1279.                 + two_plus_em2 * (u[ij + move_over[9]] - u[ij + move_over[2]]) );
  1280.     }
  1281.         
  1282.     for (j = 0; j < ny; j ++) {
  1283.             
  1284.         ij = ij_sw_corner + j;
  1285.         u[ij+move_over[4]] = 
  1286.             u[ij + move_over[7]] + eps_p2 * (u[ij + move_over[3]] + u[ij + move_over[10]]
  1287.             -u[ij + move_over[1]] - u[ij + move_over[8]])
  1288.             + two_plus_ep2 * (u[ij + move_over[5]] - u[ij + move_over[6]]);
  1289.                 
  1290.         ij = ij_se_corner + j;
  1291.         u[ij + move_over[7]] = 
  1292.             - (-u[ij + move_over[4]] + eps_p2 * (u[ij + move_over[3]] + u[ij + move_over[10]]
  1293.             - u[ij + move_over[1]] - u[ij + move_over[8]])
  1294.             + two_plus_ep2 * (u[ij + move_over[5]] - u[ij + move_over[6]]) );
  1295.     }
  1296.  
  1297.     /* That resets the boundary values.  Now we can test all data.  
  1298.         Note that this loop checks all values, even though only nearest were used.  */
  1299.     
  1300.     for (k = 0; k < npoints; k++) {
  1301.         i = data[k].index/ny;
  1302.         j = data[k].index%ny;
  1303.          ij = ij_sw_corner + i * my + j;
  1304.          if ( iu[ij] == 5 ) continue;
  1305.          x0 = xmin + i*xinc;
  1306.          y0 = ymin + j*yinc;
  1307.          dx = (data[k].x - x0)*r_xinc;
  1308.          dy = (data[k].y - y0)*r_yinc;
  1309.  
  1310.          du_dx = 0.5 * (u[ij + move_over[6]] - u[ij + move_over[5]]);
  1311.          du_dy = 0.5 * (u[ij + move_over[2]] - u[ij + move_over[9]]);
  1312.          d2u_dx2 = u[ij + move_over[6]] + u[ij + move_over[5]] - 2 * u[ij];
  1313.          d2u_dy2 = u[ij + move_over[2]] + u[ij + move_over[9]] - 2 * u[ij];
  1314.          d2u_dxdy = 0.25 * (u[ij + move_over[3]] - u[ij + move_over[1]]
  1315.                  - u[ij + move_over[10]] + u[ij + move_over[8]]);
  1316.          d3u_dx3 = 0.5 * ( u[ij + move_over[7]] - 2 * u[ij + move_over[6]]
  1317.                      + 2 * u[ij + move_over[5]] - u[ij + move_over[4]]);
  1318.          d3u_dy3 = 0.5 * ( u[ij + move_over[0]] - 2 * u[ij + move_over[2]]
  1319.                      + 2 * u[ij + move_over[9]] - u[ij + move_over[11]]);
  1320.          d3u_dx2dy = 0.5 * ( ( u[ij + move_over[3]] + u[ij + move_over[1]] - 2 * u[ij + move_over[2]] )
  1321.                      - ( u[ij + move_over[10]] + u[ij + move_over[8]] - 2 * u[ij + move_over[9]] ) );
  1322.          d3u_dxdy2 = 0.5 * ( ( u[ij + move_over[3]] + u[ij + move_over[10]] - 2 * u[ij + move_over[6]] )
  1323.                      - ( u[ij + move_over[1]] + u[ij + move_over[8]] - 2 * u[ij + move_over[5]] ) );
  1324.  
  1325.          /* 3rd order Taylor approx:  */
  1326.              
  1327.          z_est = u[ij] + dx * (du_dx +  dx * ( (0.5 * d2u_dx2) + dx * (d3u_dx3 / 6.0) ) )
  1328.                 + dy * (du_dy +  dy * ( (0.5 * d2u_dy2) + dy * (d3u_dy3 / 6.0) ) )
  1329.                  + dx * dy * (d2u_dxdy) + (0.5 * dx * d3u_dx2dy) + (0.5 * dy * d3u_dxdy2);
  1330.              
  1331.          z_err = z_est - data[k].z;
  1332.          mean_error += z_err;
  1333.          mean_squared_error += (z_err * z_err);
  1334.      }
  1335.      mean_error /= npoints;
  1336.      mean_squared_error = sqrt( mean_squared_error / npoints);
  1337.      
  1338.      curvature = 0.0;
  1339.      n_nodes = nx * ny;
  1340.      
  1341.      for (i = 0; i < nx; i++) {
  1342.          for (j = 0; j < ny; j++) {
  1343.              ij = ij_sw_corner + i * my + j;
  1344.              c = u[ij + move_over[6]] + u[ij + move_over[5]]
  1345.                  + u[ij + move_over[2]] + u[ij + move_over[9]] - 4.0 * u[ij + move_over[6]];
  1346.             curvature += (c * c);
  1347.         }
  1348.     }
  1349.  
  1350.      fprintf(stderr, "Fit info: N data points  N nodes\tmean error\trms error\tcurvature\n");
  1351.      fprintf (stderr,"\t%8d\t%8d\t%.8lg\t%.8lg\t%.8lg\n", npoints, n_nodes, mean_error, mean_squared_error,
  1352.          curvature);
  1353.  }
  1354.  
  1355. int    remove_planar_trend()
  1356. {
  1357.  
  1358.     int    i;
  1359.     double    a, b, c, d, xx, yy, zz;
  1360.     double    sx, sy, sz, sxx, sxy, sxz, syy, syz;
  1361.     
  1362.     sx = sy = sz = sxx = sxy = sxz = syy = syz = 0.0;
  1363.     
  1364.     for (i = 0; i < npoints; i++) {
  1365.  
  1366.         xx = (data[i].x - xmin) * r_xinc;
  1367.         yy = (data[i].y - ymin) * r_yinc;
  1368.         zz = data[i].z;
  1369.         
  1370.         sx += xx;
  1371.         sy += yy;
  1372.         sz += zz;
  1373.         sxx +=(xx * xx);
  1374.         sxy +=(xx * yy);
  1375.         sxz +=(xx * zz);
  1376.         syy +=(yy * yy);
  1377.         syz +=(yy * zz);
  1378.     }
  1379.     
  1380.     d = npoints*sxx*syy + 2*sx*sy*sxy - npoints*sxy*sxy - sx*sx*syy - sy*sy*sxx;
  1381.     
  1382.     if (d == 0.0) {
  1383.         plane_c0 = plane_c1 = plane_c2 = 0.0;
  1384.         return(0);
  1385.     }
  1386.     
  1387.     a = sz*sxx*syy + sx*sxy*syz + sy*sxy*sxz - sz*sxy*sxy - sx*sxz*syy - sy*syz*sxx;
  1388.     b = npoints*sxz*syy + sz*sy*sxy + sy*sx*syz - npoints*sxy*syz - sz*sx*syy - sy*sy*sxz;
  1389.     c = npoints*sxx*syz + sx*sy*sxz + sz*sx*sxy - npoints*sxy*sxz - sx*sx*syz - sz*sy*sxx;
  1390.  
  1391.     plane_c0 = a / d;
  1392.     plane_c1 = b / d;
  1393.     plane_c2 = c / d;
  1394.  
  1395.     for (i = 0; i < npoints; i++) {
  1396.  
  1397.         xx = (data[i].x - xmin) * r_xinc;
  1398.         yy = (data[i].y - ymin) * r_yinc;
  1399.         
  1400.         data[i].z -=(plane_c0 + plane_c1 * xx + plane_c2 * yy);
  1401.     }
  1402.  
  1403.     return(0);
  1404. }
  1405.  
  1406. int    replace_planar_trend()
  1407. {
  1408.     int    i, j, ij;
  1409.  
  1410.      for (i = 0; i < nx; i++) {
  1411.          for (j = 0; j < ny; j++) {
  1412.              ij = ij_sw_corner + i * my + j;
  1413.              u[ij] = (u[ij] * z_scale) + (plane_c0 + plane_c1 * i + plane_c2 * j);
  1414.         }
  1415.     }
  1416.     return(0);
  1417. }
  1418.  
  1419. int    throw_away_unusables()
  1420. {
  1421.     /* This is a new routine to eliminate data which will become
  1422.         unusable on the final iteration, when grid = 1.
  1423.         It assumes grid = 1 and set_grid_parameters has been
  1424.         called.  We sort, mark redundant data as OUTSIDE, and
  1425.         sort again, chopping off the excess.
  1426.         
  1427.         Experimental modification 5 Dec 1988 by Smith, as part
  1428.         of a new implementation using core memory for b[6]
  1429.         coefficients, eliminating calls to temp file.
  1430.     */
  1431.     
  1432.     int    last_index, n_outside, k;
  1433.     
  1434.     /* Sort the data  */
  1435.     
  1436.     qsort ((char *)data, npoints, sizeof (struct DATA), compare_points);
  1437.     
  1438.     /* If more than one datum is indexed to same node, only the first should be kept.
  1439.         Mark the additional ones as OUTSIDE
  1440.     */
  1441.     last_index = -1;
  1442.     n_outside = 0;
  1443.     for (k = 0; k < npoints; k++) {
  1444.         if (data[k].index == last_index) {
  1445.             data[k].index = OUTSIDE;
  1446.             n_outside++;
  1447.         }
  1448.         else {
  1449.             last_index = data[k].index;
  1450.         }
  1451.     }
  1452.     /* Sort again; this time the OUTSIDE points will be thrown away  */
  1453.     
  1454.     qsort ((char *)data, npoints, sizeof (struct DATA), compare_points);
  1455.     npoints -= n_outside;
  1456.     data = (struct DATA *) memory ((char *)data, npoints, sizeof(struct DATA), "surface");
  1457.     if (gmtdefs.verbose && (n_outside)) {
  1458.         fprintf(stderr,"surface: %d unusable points were supplied; these will be ignored.\n", n_outside);
  1459.         fprintf(stderr,"\tYou should have pre-processed the data with blockmean or blockmedian.\n");
  1460.     }
  1461.  
  1462.     return(0);
  1463. }
  1464.  
  1465. int    rescale_z_values()
  1466. {
  1467.     int    i;
  1468.     double    ssz = 0.0;
  1469.  
  1470.     for (i = 0; i < npoints; i++) {
  1471.         ssz += (data[i].z * data[i].z);
  1472.     }
  1473.     
  1474.     /* Set z_scale = rms(z):  */
  1475.     
  1476.     z_scale = sqrt(ssz / npoints);
  1477.     r_z_scale = 1.0 / z_scale;
  1478.  
  1479.     for (i = 0; i < npoints; i++) {
  1480.         data[i].z *= r_z_scale;
  1481.     }
  1482.     return (0);
  1483. }
  1484.  
  1485. int load_constraints (low, high)
  1486. char *low, *high; {
  1487.     int i, j, ij, n_trimmed, dummy[4];
  1488.     double yy;
  1489.     struct GRD_HEADER hdr;
  1490.     
  1491.     dummy[3] = dummy[2] = dummy[1] = dummy[0] = 0;
  1492.     
  1493.     /* Load lower/upper limits, verify range, deplane, and rescale */
  1494.     
  1495.     if (set_low > 0) {
  1496.         lower = (float *) memory (CNULL, nx * ny, sizeof (float), "surface");
  1497.         if (set_low < 3)
  1498.             for (i = 0; i < nx * ny; i++) lower[i] = low_limit;
  1499.         else {
  1500.             if (read_grd_info (low, &hdr)) {
  1501.                 fprintf (stderr, "surface: Error opening file %s\n", low);
  1502.                 exit (-1);
  1503.             }
  1504.             if (hdr.nx != nx || hdr.ny != ny) {
  1505.                 fprintf (stderr, "surface: lower limit file not of proper dimension!\n");
  1506.                 exit (-1);
  1507.             }
  1508.             if (read_grd (low, &hdr, lower, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  1509.                 fprintf (stderr, "surface: Error reading file %s\n", low);
  1510.                 exit (-1);
  1511.             }
  1512. /* Comment this out:    n_trimmed = 0;
  1513.             for (i = 0; i < nx * ny; i++) if (lower[i] > low_limit) {
  1514.                 lower[i] = low_limit;
  1515.                 n_trimmed++;
  1516.             }
  1517.             if (n_trimmed) fprintf (stderr, "surface: %d lower limit values > min data, reset to min data!\n");
  1518. */
  1519.         }
  1520.             
  1521.         for (j = ij = 0; j < ny; j++) {
  1522.             yy = ny - j - 1;
  1523.             for (i = 0; i < nx; i++, ij++) {
  1524.                 if (bad_float ((double)lower[ij])) continue;
  1525.                 lower[ij] -= (plane_c0 + plane_c1 * i + plane_c2 * yy);
  1526.                 lower[ij] *= r_z_scale;
  1527.             }
  1528.         }
  1529.         constrained = TRUE;
  1530.     }
  1531.     if (set_high > 0) {
  1532.         upper = (float *) memory (CNULL, nx * ny, sizeof (float), "surface");
  1533.         if (set_high < 3)
  1534.             for (i = 0; i < nx * ny; i++) upper[i] = high_limit;
  1535.         else {
  1536.             if (read_grd_info (high, &hdr)) {
  1537.                 fprintf (stderr, "surface: Error opening file %s\n", high);
  1538.                 exit (-1);
  1539.             }
  1540.             if (hdr.nx != nx || hdr.ny != ny) {
  1541.                 fprintf (stderr, "surface: upper limit file not of proper dimension!\n");
  1542.                 exit (-1);
  1543.             }
  1544.             if (read_grd (high, &hdr, upper, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  1545.                 fprintf (stderr, "surface: Error reading file %s\n", high);
  1546.                 exit (-1);
  1547.             }
  1548. /* Comment this out:    n_trimmed = 0;
  1549.             for (i = 0; i < nx * ny; i++) if (upper[i] < high_limit) {
  1550.                 upper[i] = high_limit;
  1551.                 n_trimmed++;
  1552.             }
  1553.             if (n_trimmed) fprintf (stderr, "surface: %d upper limit values < max data, reset to max data!\n");
  1554. */
  1555.         }
  1556.         for (j = ij = 0; j < ny; j++) {
  1557.             yy = ny - j - 1;
  1558.             for (i = 0; i < nx; i++, ij++) {
  1559.                 if (bad_float ((double)upper[ij])) continue;
  1560.                 upper[ij] -= (plane_c0 + plane_c1 * i + plane_c2 * yy);
  1561.                 upper[ij] *= r_z_scale;
  1562.             }
  1563.         }
  1564.         constrained = TRUE;
  1565.     }
  1566. }
  1567.  
  1568. double    guess_surface_time(nx, ny)
  1569. int    nx, ny;
  1570. {
  1571.     /* Routine to guess a number proportional to the operations
  1572.      * required by surface working on a user-desired grid of
  1573.      * size nx by ny, where nx = (xmax - xmin)/dx, and same for
  1574.      * ny.  (That is, one less than actually used in routine.)
  1575.      *
  1576.      * This is based on the following untested conjecture:
  1577.      *     The operations are proportional to T = nxg*nyg*L,
  1578.      *    where L is a measure of the distance that data
  1579.      *    constraints must propagate, and nxg, nyg are the
  1580.      *     current size of the grid.
  1581.      *    For nx,ny relatively prime, we will go through only
  1582.      *     one grid cycle, L = max(nx,ny), and T = nx*ny*L.
  1583.      *    But for nx,ny whose greatest common divisor is a highly
  1584.      *     composite number, we will have L equal to the division
  1585.      *     step made at each new grid cycle, and nxg,nyg will
  1586.      *     also be smaller than nx,ny.  Thus we can hope to find
  1587.      *    some nx,ny for which the total value of T is small.
  1588.      *
  1589.      * The above is pure speculation and has not been derived
  1590.      * empirically.  In actual practice, the distribution of the
  1591.      * data, both spatially and in terms of their values, will
  1592.      * have a strong effect on convergence.
  1593.      *
  1594.      * W. H. F. Smith, 26 Feb 1992.  */
  1595.  
  1596.     int    gcd_euclid();    /* Finds the greatest common divisor  */
  1597.     int    get_prime_factors();
  1598.     int    gcd;        /* Current value of the gcd  */
  1599.     int    nxg, nyg;    /* Current value of the grid dimensions  */
  1600.     int    nfactors;    /* Number of prime factors of current gcd  */
  1601.     int    factor;        /* Currently used factor  */
  1602.     /* Doubles are used below, even though the values will be integers,
  1603.         because the multiplications might reach sizes of O(n**3)  */
  1604.     double    t_sum;        /* Sum of values of T at each grid cycle  */
  1605.     double    length;        /* Current propagation distance.  */
  1606.  
  1607.  
  1608.     gcd = gcd_euclid(nx, ny);
  1609.     if (gcd > 1) {
  1610.         nfactors = get_prime_factors(gcd, factors);
  1611.         nxg = nx/gcd;
  1612.         nyg = ny/gcd;
  1613.         if (nxg < 3 || nyg < 3) {
  1614.             factor = factors[nfactors - 1];
  1615.             nfactors--;
  1616.             gcd /= factor;
  1617.             nxg *= factor;
  1618.             nyg *= factor;
  1619.         }
  1620.     }
  1621.     else {
  1622.         nxg = nx;
  1623.         nyg = ny;
  1624.     }
  1625.     length = MAX(nxg, nyg);
  1626.     t_sum = nxg * (nyg * length);    /* Make it double at each multiply  */
  1627.  
  1628.     /* Are there more grid cycles ?  */
  1629.     while (gcd > 1) {
  1630.         factor = factors[nfactors - 1];
  1631.         nfactors--;
  1632.         gcd /= factor;
  1633.         nxg *= factor;
  1634.         nyg *= factor;
  1635.         length = factor;
  1636.         t_sum += nxg * (nyg * length);
  1637.     }
  1638.     return(t_sum);
  1639. }
  1640.  
  1641.  
  1642. void    suggest_sizes_for_surface(nx, ny)
  1643. int    nx, ny;
  1644. {
  1645.     /* Calls guess_surface_time for a variety of trial grid
  1646.      * sizes, where the trials are highly composite numbers
  1647.      * with lots of factors of 2, 3, and 5.  The sizes are
  1648.      * within the range (nx,ny) - (2*nx, 2*ny).  Prints to
  1649.      * stderr the values which are an improvement over the
  1650.      * user's original nx,ny.
  1651.      * Should be called with nx=(xmax-xmin)/dx, and ditto
  1652.      * for ny; that is, one smaller than the lattice used
  1653.      * in surface.c
  1654.      *
  1655.      * W. H. F. Smith, 26 Feb 1992.  */
  1656.  
  1657.     double    guess_surface_time();
  1658.     double    users_time;    /* Time for user's nx, ny  */
  1659.     double    current_time;    /* Time for current nxg, nyg  */
  1660.     int    i;
  1661.     int    nxg, nyg;    /* Guessed by this routine  */
  1662.     int    nx2, ny2, nx3, ny3, nx5, ny5;    /* For powers  */
  1663.     int    xstop, ystop;    /* Set to 2*nx, 2*ny  */
  1664.     int    n_sug = 0;    /* N of suggestions found  */
  1665.     int    compare_sugs();    /* Sort suggestions decreasing  */
  1666.     struct SUGGESTION *sug = NULL;
  1667.     
  1668.     users_time = guess_surface_time(nx, ny);
  1669.     xstop = 2*nx;
  1670.     ystop = 2*ny;
  1671.  
  1672.     for (nx2 = 2; nx2 <= xstop; nx2 *= 2) {
  1673.       for (nx3 = 1; nx3 <= xstop; nx3 *= 3) {
  1674.         for (nx5 = 1; nx5 <= xstop; nx5 *= 5) {
  1675.         nxg = nx2 * nx3 * nx5;
  1676.         if (nxg < nx || nxg > xstop) continue;
  1677.  
  1678.         for (ny2 = 2; ny2 <= ystop; ny2 *= 2) {
  1679.           for (ny3 = 1; ny3 <= ystop; ny3 *= 3) {
  1680.             for (ny5 = 1; ny5 <= ystop; ny5 *= 5) {
  1681.             nyg = ny2 * ny3 * ny5;
  1682.             if (nyg < ny || nyg > ystop) continue;
  1683.  
  1684.             current_time = guess_surface_time(nxg, nyg);
  1685.             if (current_time < users_time) {
  1686.                 n_sug++;
  1687.                 sug = (struct SUGGESTION *)memory((char *)sug, n_sug, sizeof(struct SUGGESTION), "suggest_sizes_for_surface");
  1688.                 sug[n_sug-1].nx = nxg;
  1689.                 sug[n_sug-1].ny = nyg;
  1690.                 sug[n_sug-1].factor = users_time/current_time;
  1691.             }
  1692.  
  1693.             }
  1694.           }
  1695.         }
  1696.  
  1697.         }
  1698.       }
  1699.     }
  1700.  
  1701.     if (n_sug) {
  1702.         qsort((char *)sug, n_sug, sizeof(struct SUGGESTION), compare_sugs);
  1703.         for (i = 0; i < n_sug && i < 10; i++) {
  1704.             fprintf(stderr,"surface:  HINT:  Choosing nx = %d, ny = %d might cut run time by a factor of %.8lg\n",
  1705.                 sug[i].nx, sug[i].ny, sug[i].factor);
  1706.         }
  1707.         free((char *)sug);
  1708.     }
  1709.     else {
  1710.         fprintf(stderr,"surface: Cannot suggest any nx,ny better than your -R -I define.\n");
  1711.     }
  1712.     return;
  1713. }
  1714.  
  1715. int    compare_sugs(point_1, point_2)
  1716. struct SUGGESTION    *point_1, *point_2;
  1717. {
  1718.     /* Sorts sugs into DESCENDING order!  */
  1719.     if (point_1->factor < point_2->factor)
  1720.         return(1);
  1721.     else if (point_1->factor > point_2->factor)
  1722.         return(-1);
  1723.     else
  1724.         return(0);
  1725. }
  1726.  
  1727. int    get_prime_factors(n, f)
  1728. int    n, f[];
  1729. {
  1730.     /* Fills the integer array f with the prime factors of n.
  1731.      * Returns the number of locations filled in f, which is
  1732.      * one if n is prime.
  1733.      *
  1734.      * f[] should have been malloc'ed to enough space before
  1735.      * calling prime_factors().  We can be certain that f[32]
  1736.      * is enough space, for if n fits in a long, then n < 2**32,
  1737.      * and so it must have fewer than 32 prime factors.  I think
  1738.      * that in general, ceil(log2((double)n)) is enough storage
  1739.      * space for f[].
  1740.      *
  1741.      * Tries 2,3,5 explicitly; then alternately adds 2 or 4
  1742.      * to the previously tried factor to obtain the next trial
  1743.      * factor.  This is done with the variable two_four_toggle.
  1744.      * With this method we try 7,11,13,17,19,23,25,29,31,35,...
  1745.      * up to a maximum of sqrt(n).  This shortened list results
  1746.      * in 1/3 fewer divisions than if we simply tried all integers
  1747.      * between 5 and sqrt(n).  We can reduce the size of the list
  1748.      * of trials by an additional 20% by removing the multiples
  1749.      * of 5, which are equal to 30m +/- 5, where m >= 1.  Starting
  1750.      * from 25, these are found by alternately adding 10 or 20.
  1751.      * To do this, we use the variable ten_twenty_toggle.
  1752.      *
  1753.      * W. H. F. Smith, 26 Feb 1992, after D.E. Knuth, vol. II  */
  1754.  
  1755.     int    current_factor;    /* The factor currently being tried  */
  1756.     int    max_factor;    /* Don't try any factors bigger than this  */
  1757.     int    n_factors = 0;    /* Returned; one if n is prime  */
  1758.     int    two_four_toggle = 0;    /* Used to add 2 or 4 to get next trial factor  */
  1759.     int    ten_twenty_toggle = 0;    /* Used to add 10 or 20 to skip_five  */
  1760.     int    skip_five = 25;    /* Used to skip multiples of 5 in the list  */
  1761.     int    m;    /* Used to keep a working copy of n  */
  1762.  
  1763.  
  1764.     /* Initialize m and max_factor  */
  1765.     m = abs(n);
  1766.     if (m < 2) return(0);
  1767.     max_factor = floor(sqrt((double)m));
  1768.  
  1769.     /* First find the 2s  */
  1770.     current_factor = 2;
  1771.     while(!(m%current_factor)) {
  1772.         m /= current_factor;
  1773.         f[n_factors] = current_factor;
  1774.         n_factors++;
  1775.     }
  1776.     if (m == 1) return(n_factors);
  1777.  
  1778.     /* Next find the 3s  */
  1779.     current_factor = 3;
  1780.     while(!(m%current_factor)) {
  1781.         m /= current_factor;
  1782.         f[n_factors] = current_factor;
  1783.         n_factors++;
  1784.     }
  1785.     if (m == 1) return(n_factors);
  1786.  
  1787.     /* Next find the 5s  */
  1788.     current_factor = 5;
  1789.     while(!(m%current_factor)) {
  1790.         m /= current_factor;
  1791.         f[n_factors] = current_factor;
  1792.         n_factors++;
  1793.     }
  1794.     if (m == 1) return(n_factors);
  1795.  
  1796.     /* Now try all the rest  */
  1797.  
  1798.     while (m > 1 && current_factor <= max_factor) {
  1799.  
  1800.         /* Current factor is either 2 or 4 more than previous value  */
  1801.  
  1802.         if (two_four_toggle) {
  1803.             current_factor += 4;
  1804.             two_four_toggle = 0;
  1805.         }
  1806.         else {
  1807.             current_factor += 2;
  1808.             two_four_toggle = 1;
  1809.         }
  1810.  
  1811.         /* If current factor is a multiple of 5, skip it.  But first,
  1812.             set next value of skip_five according to 10/20 toggle:  */
  1813.  
  1814.         if (current_factor == skip_five) {
  1815.             if (ten_twenty_toggle) {
  1816.                 skip_five += 20;
  1817.                 ten_twenty_toggle = 0;
  1818.             }
  1819.             else {
  1820.                 skip_five += 10;
  1821.                 ten_twenty_toggle = 1;
  1822.             }
  1823.             continue;
  1824.         }
  1825.  
  1826.         /* Get here when current_factor is not a multiple of 2,3 or 5:  */
  1827.  
  1828.         while(!(m%current_factor)) {
  1829.             m /= current_factor;
  1830.             f[n_factors] = current_factor;
  1831.             n_factors++;
  1832.         }
  1833.     }
  1834.  
  1835.     /* Get here when all factors up to floor(sqrt(n)) have been tried.  */
  1836.  
  1837.     if (m > 1) {
  1838.         /* m is an additional prime factor of n  */
  1839.         f[n_factors] = m;
  1840.         n_factors++;
  1841.     }
  1842.     return (n_factors);
  1843. }
  1844. /* gcd_euclid.c  Greatest common divisor routine  */
  1845.  
  1846. #define IABS(i)    (((i) < 0) ? -(i) : (i))
  1847.  
  1848. int    gcd_euclid(a,b)
  1849. int    a, b;
  1850. {
  1851.     /* Returns the greatest common divisor of u and v by Euclid's method.
  1852.      * I have experimented also with Stein's method, which involves only
  1853.      * subtraction and left/right shifting; Euclid is faster, both for
  1854.      * integers of size 0 - 1024 and also for random integers of a size
  1855.      * which fits in a long integer.  Stein's algorithm might be better
  1856.      * when the integers are HUGE, but for our purposes, Euclid is fine.
  1857.      *
  1858.      * Walter H. F. Smith, 25 Feb 1992, after D. E. Knuth, vol. II  */
  1859.  
  1860.     int    u,v,r;
  1861.  
  1862.     u = MAX(IABS(a), IABS(b));
  1863.     v = MIN(IABS(a), IABS(b));
  1864.  
  1865.     while (v > 0) {
  1866.         r = u%v;    /* Knuth notes that u < 2v 40% of the time;  */
  1867.         u = v;        /* thus we could have tried a subtraction  */
  1868.         v = r;        /* followed by an if test to do r = u%v  */
  1869.     }
  1870.     return(u);
  1871. }
  1872.  
  1873.