home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / sample1d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-11  |  9.2 KB  |  317 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)sample1d.c    2.20  4/11/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * sample1d reads a 1-D dataset, and resamples the values on (1) user
  9.  * supplied time values or (2) equidistant time-values based on <timestart> and
  10.  * <dt>, both supplied at the command line. Choose among linear, cubic
  11.  * spline, and Akima's spline.  sample1d will handle multiple column files,
  12.  * user must choose which column contains the independent, monotonically
  13.  * increasing variable.
  14.  *
  15.  * Author:    Paul Wessel
  16.  * Date:    5-JAN-1990
  17.  * Version:    v2.0
  18.  *
  19.  */
  20.  
  21. #include "gmt.h"
  22.  
  23. double *t_out, *tmp, *ttime, *data, **col, **out;
  24. BOOLEAN *nan_flag;
  25.  
  26. main (argc, argv)
  27. int argc;
  28. char **argv; {
  29.     int i, j, k, t = 0, n, m, n_alloc, m_alloc, result, n_col;
  30.     
  31.     BOOLEAN error = FALSE, equidist = FALSE, first = TRUE, t_supplied = FALSE;
  32.     BOOLEAN got_t0 = FALSE, got_dt = FALSE, multi_segments = FALSE, ok;
  33.     
  34.     double dt = 0.0, t0 = 0.0, value, tt;
  35.     
  36.     char line[512], test[512], *ptr, type[3], *more, EOL_flag = '>';
  37.     
  38.     FILE *fpi = NULL, *fpt = NULL;
  39.     
  40.     type[0] = 'L';    type[1] = 'A';    type[2] = 'C';
  41.     
  42.     argc = gmt_begin (argc, argv);
  43.  
  44.     for (i = 1; i < argc; i++) {
  45.         if (argv[i][0] == '-') {
  46.             switch (argv[i][1]) {
  47.         
  48.                 /* Common parameters */
  49.             
  50.                 case 'H':
  51.                 case '\0':
  52.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  53.                     break;
  54.  
  55.                 /* Supplemental parameters */
  56.                 
  57.                 case 'A':
  58.                     gmtdefs.interpolant = 1;
  59.                     break;
  60.                 case 'C':
  61.                     gmtdefs.interpolant = 2;
  62.                     break;
  63.                 case 'I':
  64.                     dt = atof (&argv[i][2]);
  65.                     equidist = TRUE;
  66.                     got_dt = TRUE;
  67.                     break;
  68.                 case 'L':
  69.                     gmtdefs.interpolant = 0;
  70.                     break;
  71.                 case 'T':
  72.                     t = atoi (&argv[i][2]);
  73.                     break;
  74.                 case 'M':               /* Multiple line segments */
  75.                     multi_segments = TRUE;
  76.                     if (argv[i][2]) EOL_flag = argv[i][2];
  77.                     break;
  78.                 case 'N':
  79.                     if ((fpt = fopen (&argv[i][2], "r")) == NULL) {
  80.                         fprintf (stderr, "sample1d: Cannot open file %s\n", &argv[i][2]);
  81.                         exit (-1);
  82.                     }
  83.                     t_supplied = TRUE;
  84.                     break;
  85.                 case 'S':
  86.                     t0 = atof (&argv[i][2]);
  87.                     equidist = TRUE;
  88.                     got_t0 = TRUE;
  89.                     break;
  90.                     
  91.                 default:
  92.                     error = TRUE;
  93.                     gmt_default_error (argv[i][1]);
  94.                     break;
  95.             }
  96.         }
  97.         else if ((fpi = fopen (argv[i], "r")) == NULL) {
  98.             fprintf (stderr, "sample1d: Cannot open file %s\n", argv[i]);
  99.             exit (-1);
  100.         }
  101.         
  102.     }
  103.  
  104.     if (argc == 1 || gmt_quick) {    /* Display usage */
  105.         fprintf (stderr,"sample1d %s - Resampling of 1-D data sets\n\n", GMT_VERSION);
  106.         fprintf(stderr, "usage: sample1d <infile> [-A] [-C] [-H] [-I<t_inc>] [-L]  [-M[<flag>]] [-N<knotfile>]\n");
  107.         fprintf(stderr, "    [-S<xstart>] [-T<time_col>]\n\n");
  108.         
  109.         if (gmt_quick) exit(-1);
  110.         
  111.         fprintf(stderr, "    <infile> is a multicolumn ASCII table. [Default is standard input]\n");
  112.         fprintf(stderr, "\n\tOPTIONS:\n");
  113.         fprintf(stderr, "    -A use Akima spline interpolation. [Default is -%c]\n", type[gmtdefs.interpolant]);
  114.         fprintf(stderr, "    -C use cubic spline interpolation. [Default is -%c]\n", type[gmtdefs.interpolant]);
  115.         explain_option ('H');
  116.         fprintf(stderr, "    -I <x_inc> sets equidistant grid interval\n");
  117.         fprintf(stderr, "    -L use linear interpolation. [Default is -%c]\n", type[gmtdefs.interpolant]);
  118.         fprintf(stderr, "    -M Input files each consist of multiple segments separated by one record.\n");
  119.         fprintf(stderr, "       whose first character is <flag> [>].\n");
  120.         fprintf(stderr, "    -N <knotfile> is an ASCII table with the desired time positions in column 0\n");
  121.         fprintf(stderr, "       Overrides the -I and -S settings.  If none of -I, -S, and -N is set\n");
  122.         fprintf(stderr, "       then <tstart> = first input point, <t_inc> = (t[1] - t[0])\n");
  123.         fprintf(stderr, "    -S <xstart> sets the first output point\n");
  124.         fprintf(stderr, "    -T gives column number of the independent variable (time) [Default is 0 (first)]\n");
  125.         explain_option ('.');
  126.         exit(-1);
  127.     }
  128.     
  129.     if (t < 0) {
  130.         fprintf (stderr, "%s: GMT SYNTAX ERROR -T option: Column number cannot be negative\n", gmt_program);
  131.         error++;
  132.     }
  133.     if (t_supplied && equidist) {
  134.         fprintf (stderr, "%s: GMT SYNTAX ERROR: Specify only one of -N and -S\n", gmt_program);
  135.         error++;
  136.     }
  137.     if (!t_supplied && !equidist) {
  138.         fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one of -N and -S\n", gmt_program);
  139.         error++;
  140.     }
  141.     if (got_dt && dt <= 0.0) {
  142.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Must specify positive increment\n", gmt_program);
  143.         error++;
  144.     }
  145.     
  146.     if (error) exit (-1);
  147.     
  148.     if (fpi == NULL) fpi = stdin;
  149.  
  150.     if (gmtdefs.io_header) {
  151.         for (i = 0; i < gmtdefs.n_header_recs; i++) {
  152.             fgets (line, 512, fpi);
  153.             printf ("%s", line);
  154.         }
  155.     }
  156.     if (multi_segments) {
  157.         fgets (line, 512, fpi);
  158.         printf ("%s", line);
  159.     }
  160.     
  161.     i = 0;
  162.     more = fgets (line, 512, fpi);
  163.     
  164.     while (more) {
  165.  
  166.         i++;
  167.         n = 0;
  168.     
  169.         if (first) {    /* Allocate memory first time around */
  170.             first = FALSE;
  171.             strcpy (test, line);
  172.             ptr = strtok (test, " \t");
  173.             n_col = 0;
  174.             while (ptr) {
  175.                 ptr = strtok (CNULL, " \t");
  176.                 n_col++;
  177.             }
  178.             if (t >= n_col) {
  179.                 fprintf (stderr, "sample1d: time_col exceeds range of columns!\n");
  180.                 exit (-1);
  181.             }
  182.             col = (double **) memory (CNULL, n_col, sizeof (double *), "sample1d");
  183.             out = (double **) memory (CNULL, n_col, sizeof (double *), "sample1d");
  184.             nan_flag = (BOOLEAN *) memory (CNULL, n_col, sizeof (BOOLEAN), "sample1d");
  185.             n_alloc = GMT_CHUNK;
  186.             for (j = 0; j < n_col; j++) col[j] = (double *) memory (CNULL, n_alloc, sizeof (double), "sample1d");
  187.         }
  188.         
  189.         while ((more && !multi_segments) || (more && multi_segments && line[0] != EOL_flag)) {
  190.         
  191.             ptr = strtok (line, " \t");
  192.             j = 0;
  193.             while (ptr && j < n_col) {
  194.                 ok = sscanf (ptr, "%lf", &value);
  195.                 if (ok)
  196.                     col[j][n] = value;
  197.                 else {
  198.                     nan_flag[j] = TRUE;
  199.                     col[j][n] = gmt_NaN;
  200.                 }
  201.                 ptr = strtok (CNULL, " \t");
  202.                 j++;
  203.             }
  204.             if (ptr || j != n_col)
  205.                 fprintf (stderr, "sample1d: Column mismatch in record # %d\n", i);
  206.             else
  207.                 n++;
  208.         
  209.             if (n == n_alloc) {    /* Get more memory */
  210.                 n_alloc += GMT_CHUNK;
  211.                 for (j = 0; j < n_col; j++) col[j] = (double *) memory ((char *)col[j], n_alloc, sizeof (double), "sample1d");
  212.             }
  213.             
  214.             more = fgets (line, 512, fpi);
  215.         }
  216.             
  217.     
  218.         if (nan_flag[t]) {
  219.             fprintf (stderr, "sample1d: Independent column has NaN's!\n");
  220.             exit (-1);
  221.         }
  222.     
  223.         if (t_supplied) {    /* read file with abscissae */
  224.             m_alloc = GMT_CHUNK;
  225.             t_out = (double *) memory ((char *)t_out, m_alloc, sizeof (double), "sample1d");
  226.             m = 0;
  227.             if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fpt);
  228.             while (fgets (line, 512, fpt)) {
  229.                 sscanf (line, "%lf", &t_out[m]);
  230.                 m++;
  231.                 if (m == m_alloc) {    /* Get more memory */
  232.                     m_alloc += GMT_CHUNK;
  233.                     t_out = (double *) memory ((char *)t_out, m_alloc, sizeof (double), "sample1d");
  234.                 }
  235.             }
  236.             fclose(fpt);
  237.         }
  238.         else {    /* Generate evenly spaced grid */
  239.             if (!got_dt) dt = col[t][1] - col[t][0];
  240.             if (!got_t0) {
  241.                 t0 = floor (col[t][0] / dt) * dt;
  242.                 if (t0 < col[t][0]) t0 += dt;
  243.             }
  244.             m = m_alloc = rint ((col[t][n-1] - t0) / dt) + 1;
  245.             if (m < 0) {
  246.                 fprintf (stderr, "sample1d: Abscissa must be monotonically increasing!\n");
  247.                 exit (-1);
  248.             }
  249.             t_out = (double *) memory (CNULL, m_alloc, sizeof (double), "sample1d");
  250.             t_out[0] = t0;
  251.             i = 1;
  252.             while (i < m && (tt = t0 + i * dt) <= col[t][n-1]) {
  253.                 t_out[i] = tt;
  254.                 i++;
  255.             }
  256.             m = i;
  257.             if (fabs (t_out[m-1]-col[t][n-1]) < SMALL) {    /* Fix roundoff */
  258.                 t_out[m-1] = col[t][n-1];
  259.             }
  260.         }
  261.     
  262.         for (j = 0; !error && j < n_col; j++) {
  263.             if (j == t) continue;
  264.  
  265.             tmp = (double *) memory (CNULL, m, sizeof (double), "sample1d");
  266.         
  267.             if (nan_flag[j]) {    /* NaN's present, need "clean" time and data columns */
  268.                 ttime = (double *) memory (CNULL, n, sizeof (double), "sample1d");
  269.                 data = (double *) memory (CNULL, n, sizeof (double), "sample1d");
  270.                 for (i = k = 0; i < n; i++) {
  271.                     if ( bad_float (col[j][i]) ) continue;
  272.                     ttime[k] = col[t][i];
  273.                     data[k] = col[j][i];
  274.                     k++;
  275.                 }
  276.                 result = intpol (ttime, data, k, m, t_out, tmp, gmtdefs.interpolant);
  277.                 free ((char *)ttime);
  278.                 free ((char *)data);
  279.             }
  280.             else
  281.                 result = intpol (col[t], col[j], n, m, t_out, tmp, gmtdefs.interpolant);
  282.                 
  283.             if (result != 0) {
  284.                 fprintf(stderr, "sample1d: Error from intpol (%d) for column %d!\n", result, j);
  285.                 error = TRUE;
  286.             }
  287.             else
  288.                 out[j] = tmp;
  289.         }
  290.     
  291.         out[t] = t_out;
  292.  
  293.         for (i = 0; !error && i < m; i++) {
  294.             printf (gmtdefs.d_format, out[0][i]);
  295.             for (j = 1; j < n_col; j++) {
  296.                 putchar ('\t');
  297.                 printf (gmtdefs.d_format, out[j][i]);
  298.             }
  299.             putchar ('\n');
  300.         }
  301.         
  302.         for (j = 1; j < n_col; j++) free ((char *)out[j]);
  303.         
  304.         if (more) {
  305.             if (multi_segments) printf ("%s", line);
  306.             more = fgets (line, 512, fpi);
  307.         }
  308.     }
  309.     
  310.     if (fpi != stdin) fclose(fpi);
  311.     
  312.     for (j = 0; j < n_col; j++) free ((char *)col[j]);
  313.     free ((char *)t_out);
  314.     
  315.     gmt_end (argc, argv);
  316. }
  317.