home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / splitxyz.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-11  |  15.9 KB  |  501 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)splitxyz.c    2.15  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.  * 
  9.  * read a file of lon, lat, zvalue[, distance, azimuth]
  10.  * and split it into profile segments.
  11.  * 
  12.  * Author:    W. H. F. Smith
  13.  * Date:    24 July, 1991-1995
  14.  */
  15.  
  16. #include "gmt.h"
  17.  
  18. #define F_RES 1000    /* Number of points in filter halfwidth  */
  19.  
  20. struct  DATA    {
  21.         double  d;
  22.         double  a;
  23.         double  x;
  24.         double  y;
  25.         double  z;
  26.         double  w;
  27. }       *data;
  28.  
  29. main(argc, argv)
  30. int     argc;
  31. char    **argv;
  32. {
  33.     int     i, ndata, n_alloc, nprofiles, begin, end, ix, iy, n_in, n_fields;
  34.     
  35.     BOOLEAN    error = FALSE, dh_supplied = FALSE, ok, hilow, map_units = FALSE;
  36.     BOOLEAN binary = FALSE, single_precision = FALSE, more, xy_only = FALSE;
  37.  
  38.         double  desired_azim = 90.0, tolerance_azim = 360.0;
  39.         double  min_dist = 100.0, xy_filter = 0.0, z_filter = 0.0;
  40.     double    d_gap = 10.0, course_change = 0.0, d_to_km;
  41.     double    dy, dx, last_d, last_c, last_s, csum, ssum, this_c, this_s, dotprod;
  42.     double    mean_azim, fwork[F_RES], *in;
  43.     
  44.      char    buffer[512], basename[120], filename[120], format[80];
  45.      
  46.     FILE    *fpin = NULL, *fpout = NULL;
  47.  
  48.     void    filterz(), filterxy_setup(), filterxy();
  49.  
  50.     argc = gmt_begin (argc, argv);
  51.     basename[0] = '\0';    /* Later, if (basename[0]) we have to write files.  */
  52.  
  53.         for (i = 1; i < argc; i++) {
  54.                 if (argv[i][0] == '-') {
  55.                         switch (argv[i][1]) {
  56.          
  57.                 /* Common parameters */
  58.             
  59.                 case 'H':
  60.                 case 'V':
  61.                 case ':':
  62.                 case '\0':
  63.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  64.                     break;
  65.                 
  66.                 /* Supplemental parameters */
  67.                 
  68.                 case 'b':
  69.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  70.                     binary = TRUE;
  71.                     break;
  72.                                case 'A':
  73.                                         if ( (sscanf(&argv[i][2], "%lf/%lf", &desired_azim, &tolerance_azim)) != 2) {
  74.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Can't dechiper values\n", gmt_program);
  75.                         error = TRUE;
  76.                     }
  77.                                         break;
  78.                                 case 'C':
  79.                                         if ( (sscanf(&argv[i][2], "%lf", &course_change)) != 1) {
  80.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Can't dechiper value\n", gmt_program);
  81.                         error = TRUE;
  82.                     }
  83.                                         break;
  84.                                 case 'D':
  85.                                         if ( (sscanf(&argv[i][2], "%lf", &min_dist)) != 1) {
  86.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option: Can't dechiper value\n", gmt_program);
  87.                         error = TRUE;
  88.                     }
  89.                                         break;
  90.                 case 'F':
  91.                                         if ( (sscanf(&argv[i][2], "%lf/%lf", &xy_filter, &z_filter)) != 2) {
  92.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Can't dechiper values\n", gmt_program);
  93.                         error = TRUE;
  94.                     }
  95.                                         break;
  96.                                 case 'G':
  97.                                         if ( (sscanf(&argv[i][2], "%lf", &d_gap)) != 1) {
  98.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Can't dechiper value\n", gmt_program);
  99.                         error = TRUE;
  100.                     }
  101.                                         break;
  102.                 case 'M':
  103.                                         map_units = TRUE;
  104.                                         break;
  105.                 case 'N':
  106.                                         strcpy(basename, &argv[i][2]);
  107.                                         break;
  108.                 case 'S':
  109.                                         dh_supplied = TRUE;
  110.                                         break;
  111.                 case 'Z':
  112.                                         xy_only = TRUE;
  113.                                         break;
  114.                                 default:
  115.                                         error = TRUE;
  116.                     gmt_default_error (argv[i][1]);
  117.                                         break;
  118.                         }
  119.                 }
  120.                 else {
  121.                         if ( (fpin = fopen(argv[i], "r")) == NULL) {
  122.                                 fprintf(stderr,"splitxyz:  Cannot open %s\n", argv[i]);
  123.                                 error = TRUE;
  124.                         }
  125.                 }
  126.         }
  127.  
  128.     if (argc == 1 || gmt_quick) {    /* Display usage */
  129.          fprintf (stderr, "splitxyz %s - Split xyz[dh] files into segments\n\n", GMT_VERSION);
  130.         fprintf(stderr,"usage:  splitxyz [<xyz[dh]file>] -C<course_change> [-A<azimuth>/<tolerance>]\n");
  131.         fprintf(stderr,"\t[-D<minimum_distance>] [-F<xy_filter>/<z_filter>] [-G<gap>] [-H] [-M]\n");
  132.         fprintf(stderr,"\t[-N<namestem>] [-S] [-V] [-Z] [-:] [-b[d]]\n\n");
  133.          
  134.         if (gmt_quick) exit (-1);
  135.         
  136.         fprintf(stderr,"\tGive xyz[dh]file name or read stdin.\n");
  137.         fprintf(stderr,"\t-C  Profile ends when change of heading exceeds <course_change>.\n");
  138.         fprintf (stderr, "\n\tOPTIONS:\n");
  139.         fprintf(stderr,"\t-A Only write profile if mean direction is w/in +/- <tolerance>\n");
  140.         fprintf(stderr,"\t     of <azimuth>. [Default = All].\n");
  141.         fprintf(stderr,"\t-D Only write profile if length is at least <minimum_distance>.\n");
  142.         fprintf(stderr,"\t     [Default = 100 dist units].\n");
  143.         fprintf(stderr,"\t-F Filter the data.  Give full widths of cosine arch filters for xy and z.\n");
  144.         fprintf(stderr,"\t   Defaults are both widths = 0, giving no filtering.\n");
  145.         fprintf(stderr,"\t   Use negative width to highpass.\n");
  146.         fprintf(stderr,"\t-G Do not let profiles have gaps exceeding <gap>. [Default = 10 dist units].\n");
  147.         explain_option ('H');
  148.         fprintf(stderr,"\t-M Map units TRUE; x,y in degrees, dist units in km.  [Default dist unit = x,y unit].\n");
  149.         fprintf(stderr,"\t-N Write output to separate files named <namestem>.profile#.\n");
  150.         fprintf(stderr,"\t     [Default all to stdout, separated by >].\n");
  151.         fprintf(stderr,"\t-S dh is supplied.  Input is 5 col x,y,z,d,h with d non-decreasing.\n");
  152.         fprintf(stderr,"\t     [Default input is 3 col x,y,z only and computes d,h from the data].\n");
  153.         fprintf(stderr,"\t-Z No z-values.  Input is 2 col x,y only\n");
  154.         explain_option ('V');
  155.         explain_option (':');
  156.         fprintf (stderr, "\t-b means binary input (output is always ascii)\n");
  157.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  158.         explain_option ('.');
  159.         exit(-1);
  160.         }
  161.  
  162.     if (course_change <= 0.0) {
  163.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option: Course change tolerance must be positive\n", gmt_program);
  164.         error++;
  165.     }
  166.     if (tolerance_azim < 0.0) {
  167.         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Azimuth tolerance must be positive\n", gmt_program);
  168.         error++;
  169.     }
  170.     if (d_gap < 0.0) {
  171.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option: Data gap distance must be positive\n", gmt_program);
  172.         error++;
  173.     }
  174.     if (desired_azim < 0.0 || desired_azim > 360) {
  175.         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: azimuth must be in 0-360 degree range\n", gmt_program);
  176.         error++;
  177.     }
  178.     if (xy_only && dh_supplied) {
  179.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option: Cannot be used with -S option\n", gmt_program);
  180.         error++;
  181.     }
  182.     if (xy_only && z_filter != 0.0) {
  183.         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: Cannot specify z-filter while using -Z option\n", gmt_program);
  184.         error++;
  185.     }
  186.  
  187.     if (error) exit (-1);
  188.     
  189.     tolerance_azim *= D2R;
  190.     if (desired_azim > 180.0) desired_azim -= 180.0;    /* Put in Easterly strikes  */
  191.     desired_azim = D2R * (90.0 - desired_azim);    /* Work in cartesian angle and radians  */
  192.     course_change *= D2R;
  193.     d_to_km = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[N_ELLIPSOIDS-1].eq_radius / 360.0;
  194.  
  195.     if (fpin == NULL) fpin = stdin;
  196.         n_alloc = GMT_CHUNK;
  197.     ndata = 0;
  198.     i = 0;
  199.     n_in = (dh_supplied) ? 5 : ((xy_only) ? 2 : 3);
  200.     if (xy_only)
  201.         sprintf (format, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);
  202.     else
  203.         sprintf (format, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  204.         data = (struct DATA *) memory (CNULL, n_alloc, sizeof(struct DATA), "splitxyz");
  205.     ix = (gmtdefs.xy_toggle);    iy = 1 - ix;
  206.  
  207.     if (binary) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  208.     if (!binary && gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (buffer, 512, fpin);
  209.  
  210.     more = ((n_fields = gmt_input (fpin,  n_in, &in)) == n_in);
  211.                 
  212.     while (more) {
  213.         i++;
  214.                 if (ndata == n_alloc) {
  215.                         n_alloc += GMT_CHUNK;
  216.                         data = (struct DATA *) memory ((char *)data, n_alloc, sizeof(struct DATA), "splitxyz");
  217.                 }
  218.         if (!dh_supplied) {
  219.             data[ndata].x = in[ix];
  220.             data[ndata].y = in[iy];
  221.             if (!xy_only) data[ndata].z = in[2];
  222.             if (ndata) {
  223.                 dy = (data[ndata].y - data[ndata-1].y);
  224.                 dx = (data[ndata].x - data[ndata-1].x);
  225.                 if (map_units) {
  226.                     dy *= d_to_km;
  227.                     dx *= (d_to_km * cos (0.5 * (data[ndata].y + data[ndata-1].y) * D2R) );
  228.                 }
  229.                 if (dy == 0.0 && dx == 0.0) {
  230.                     data[ndata].d = data[ndata-1].d;
  231.                     data[ndata].a = data[ndata-1].a;
  232.                 }
  233.                 else {
  234.                     data[ndata].d = data[ndata-1].d + hypot(dx,dy);
  235.                     data[ndata].a = d_atan2(dy,dx);
  236.                 }
  237.             }
  238.             else {
  239.                 data[ndata].d = data[ndata].a = 0.0;
  240.             }
  241.             ndata++;
  242.         }
  243.         else {
  244.             data[ndata].x = in[ix];
  245.             data[ndata].y = in[iy];
  246.             data[ndata].z = in[2];
  247.             data[ndata].d = in[3];
  248.             data[ndata].a = D2R * (90.0 - in[4]);
  249.             ndata++;
  250.         }
  251.         more = ((n_fields = gmt_input (fpin,  n_in, &in)) == n_in);
  252.     }
  253.     
  254.     data = (struct DATA *) memory ((char *)data, ndata, sizeof(struct DATA), "splitxyz");
  255.     if (!dh_supplied) data[0].a = data[1].a;
  256.     if (fpin != stdin) fclose(fpin);
  257.     if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  258.     if (n_fields%n_in) {    /* Got garbage or multiple segment subheader */
  259.         fprintf (stderr, "%s:  Cannot read line %d, aborts\n", gmt_program, ndata);
  260.         exit (-1);
  261.     }
  262.  
  263.     /* Now we have read the data and can filter z if necessary.  */
  264.  
  265.     if (z_filter < 0.0) {
  266.         hilow = TRUE;
  267.         z_filter = fabs(z_filter);
  268.         filterz(data, ndata, z_filter, fwork, hilow);
  269.     }
  270.     else if (z_filter > 0.0) {
  271.         hilow = FALSE;
  272.         filterz(data, ndata, z_filter, fwork, hilow);
  273.     }
  274.  
  275.     if (xy_filter < 0.0) {
  276.         hilow = TRUE;
  277.         xy_filter = fabs(xy_filter);
  278.         filterxy_setup(fwork);
  279.     }
  280.     else if (xy_filter > 0.0) {
  281.         hilow = FALSE;
  282.         filterxy_setup(fwork);
  283.     }
  284.  
  285.     /* Now we are ready to search for segments.  */
  286.  
  287.     nprofiles = 0;
  288.     begin = 0;
  289.     end = begin;
  290.     while (end < ndata-1) {
  291.         last_d = data[begin].d;
  292.         last_c = cos(data[begin].a);
  293.         last_s = sin(data[begin].a);
  294.         csum = last_c;
  295.         ssum = last_s;
  296.         ok = TRUE;
  297.         while (ok && end < ndata-1) {
  298.             end++;
  299.             if (data[end].d - last_d > d_gap) {
  300.                 /* Fails due to too much distance gap  */
  301.                 ok = FALSE;
  302.                 continue;
  303.             }
  304.             this_c = cos(data[end].a);
  305.             this_s = sin(data[end].a);
  306.             dotprod = this_c * last_c + this_s * last_s;
  307.             if (fabs(dotprod) > 1.0) dotprod = copysign(1.0, dotprod);
  308.             if (d_acos(dotprod) > course_change) {
  309.                 /* Fails due to too much change in azimuth  */
  310.                 ok = FALSE;
  311.                 continue;
  312.             }
  313.             /* Get here when this point belongs with last one:  */
  314.             csum += this_c;
  315.             ssum += this_s;
  316.             last_c = this_c;
  317.             last_s = this_s;
  318.             last_d = data[end].d;
  319.         }
  320.         
  321.         /* Get here when we have found a beginning and end  */
  322.  
  323.         if (end - begin - 1) {
  324.         
  325.             /* There are at least two points in the list.  */
  326.             
  327.             if ((data[end-1].d - data[begin].d) >= min_dist) {
  328.             
  329.                 /* List is long enough.  Check strike. Compute mean_azim in range [-pi/2, pi/2]:  */
  330.                 
  331.                 mean_azim = d_asin(ssum/sqrt(csum*csum + ssum*ssum));
  332.                 mean_azim = fabs(mean_azim - desired_azim);
  333.                 if (mean_azim <= tolerance_azim || mean_azim >= (M_PI-tolerance_azim) ) {
  334.                 
  335.                     /* List has acceptable strike.  */
  336.                     
  337.                     if (xy_filter != 0.0) filterxy(data, begin, end, xy_filter, fwork, hilow);
  338.                     nprofiles++;
  339.                     
  340.                     /* If (basename) we write separate files, else stdout with > marks:  */
  341.                     
  342.                     if (basename[0]) {
  343.                         sprintf(filename, "%s.profile%d\0", basename, nprofiles);
  344.                         if ( (fpout = fopen(filename, "w")) == NULL) {
  345.                             fprintf(stderr,"splitxyz:  Cannot create %s\n", argv[i]);
  346.                             exit(-1);
  347.                         }
  348.                         if (xy_only)
  349.                             for(i = begin; i < end; i++) fprintf(fpout, format, data[i].x, data[i].y);
  350.                         else
  351.                             for(i = begin; i < end; i++) fprintf(fpout, format, data[i].x, data[i].y, data[i].z);
  352.                         fclose(fpout);
  353.                         if (gmtdefs.verbose) fprintf(stderr,"Wrote %d points to file %s\n", end-begin, filename);
  354.                     }
  355.                     else {
  356.                         printf("> Start profile # %d\n", nprofiles);
  357.                         if (xy_only)
  358.                             for (i = begin; i < end; i++) printf(format, data[i].x, data[i].y);
  359.                         else
  360.                             for (i = begin; i < end; i++) printf(format, data[i].x, data[i].y, data[i].z);
  361.                     }
  362.                 }
  363.             }
  364.         }
  365.         begin = end;
  366.     }
  367.  
  368.     /* Get here when all profiles have been found and written.  */
  369.  
  370.     if (gmtdefs.verbose) fprintf(stderr,"splitxyz:  Split %d data into %d files.\n", ndata, nprofiles);
  371.     free((char *)data);
  372.     
  373.     gmt_end (argc, argv);
  374. }
  375.  
  376.  
  377. void    filterz(data, ndata, z_filter, fwork, hilow)
  378. struct DATA    *data;
  379. int    ndata, hilow;
  380. double    z_filter, *fwork;
  381. {
  382.     
  383.     double    half_width, sum, dt, tmp;
  384.     int    i, j, k, istart, istop;
  385.  
  386.     half_width = 0.5 * z_filter;
  387.     sum = 0.0;
  388.     dt = F_RES/half_width;
  389.     tmp = M_PI / F_RES;
  390.     for (i = 0; i < F_RES; i++) {
  391.         fwork[i] = 1.0 + cos(i*tmp);
  392.         sum += fwork[i];
  393.     }
  394.     for (i = 1; i < F_RES; i++) {
  395.         fwork[i] /= sum;
  396.     }
  397.  
  398.     j = 0;
  399.     istart = 0;
  400.     istop = 0;
  401.     while(j < ndata) {
  402.         while(istart < ndata && data[istart].d - data[j].d <= -half_width) istart++;
  403.         while(istop < ndata && data[istop].d - data[j].d < half_width) istop++;
  404.         sum = 0.0;
  405.         data[j].w = 0.0;
  406.         for(i = istart; i < istop; i++) {
  407.             k = floor(dt*fabs(data[i].d - data[j].d));
  408.             sum += fwork[k];
  409.             data[j].w += (data[i].z * fwork[k]);
  410.         }
  411.         data[j].w /= sum;
  412.         j++;
  413.     }
  414.     if (hilow) {
  415.         for (i = 0; i < ndata; i++) data[i].z = data[i].z - data[i].w;
  416.     }
  417.     else {
  418.         for (i = 0; i < ndata; i++) data[i].z = data[i].w;
  419.     }
  420.     return;
  421. }
  422.  
  423. void    filterxy_setup(fwork)
  424. double    *fwork;
  425. {
  426.     int    i;
  427.     double    tmp, sum = 0.0;
  428.  
  429.     tmp = M_PI / F_RES;
  430.     for (i = 0; i < F_RES; i++) {
  431.         fwork[i] = 1.0 + cos(i*tmp);
  432.         sum += fwork[i];
  433.     }
  434.     for (i = 1; i < F_RES; i++) {
  435.         fwork[i] /= sum;
  436.     }
  437.     return;
  438. }
  439.  
  440. void    filterxy(data, begin, end, xy_filter, fwork, hilow)
  441. struct DATA    *data;
  442. int    begin, end, hilow;
  443. double    xy_filter, *fwork;
  444. {
  445.     int    i, j, k, istart, istop;
  446.     double    half_width, dt, sum;
  447.  
  448.     half_width = 0.5 * fabs(xy_filter);
  449.     dt = F_RES/half_width;
  450.  
  451.     j = begin;
  452.     istart = begin;
  453.     istop = begin;
  454.     while(j < end) {
  455.         while(istart < end && data[istart].d - data[j].d <= -half_width) istart++;
  456.         while(istop < end && data[istop].d - data[j].d < half_width) istop++;
  457.         sum = 0.0;
  458.         data[j].w = 0.0;
  459.         for(i = istart; i < istop; i++) {
  460.             k = floor(dt*fabs(data[i].d - data[j].d));
  461.             sum += fwork[k];
  462.             data[j].w += (data[i].x * fwork[k]);
  463.         }
  464.         data[j].w /= sum;
  465.         j++;
  466.     }
  467.     if (hilow) {
  468.         for (i = begin; i < end; i++) data[i].x = data[i].x - data[i].w;
  469.     }
  470.     else {
  471.         for (i = begin; i < end; i++) data[i].x = data[i].w;
  472.     }
  473.  
  474.     j = begin;
  475.     istart = begin;
  476.     istop = begin;
  477.     while(j < end) {
  478.         while(istart < end && data[istart].d - data[j].d <= -half_width) istart++;
  479.         while(istop < end && data[istop].d - data[j].d < half_width) istop++;
  480.         sum = 0.0;
  481.         data[j].w = 0.0;
  482.         for(i = istart; i < istop; i++) {
  483.             k = floor(dt*fabs(data[i].d - data[j].d));
  484.             sum += fwork[k];
  485.             data[j].w += (data[i].y * fwork[k]);
  486.         }
  487.         data[j].w /= sum;
  488.         j++;
  489.     }
  490.     if (hilow) {
  491.         for (i = begin; i < end; i++) data[i].y = data[i].y - data[i].w;
  492.     }
  493.     else {
  494.         for (i = begin; i < end; i++) data[i].y = data[i].w;
  495.     }
  496.  
  497.     return;
  498. }
  499.  
  500.  
  501.