/*-------------------------------------------------------------------- * $Id: backtracker.c,v 1.2.4.3 2002/02/27 17:58:55 pwessel Exp $ * * Copyright (c) 1999-2002 by P. Wessel * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: www.soest.hawaii.edu/wessel *--------------------------------------------------------------------*/ /* * Program for moving points along small circles on a sphere given a * set of plate motion stage (Euler) poles. * backtracker can move a point forward or backward in time. * It can do so either along flowlines or hotspot tracks. * It can move point to final position or generate a track between * starting point and final position. The output track is a GMT * multisegment file and can be plotted with psxy -M. * * Author: Paul Wessel, SOEST, Univ. of Hawaii, Honolulu, HI, USA * Date: 29-DEC-1999 * Version: 1.1 * *------------------------------------------------------------------------- * The ASCII Euler file must have following format: * * 1. Any number of comment lines starting with # in first column * 2. Any number of blank lines (just carriage return, no spaces) * 2. Any number of stage pole records which each have the format: * lon(deg) lat(deg) tstart(Ma) tstop(Ma) ccw-angle(deg) * 3. stage records must go from oldest to youngest rotation * 4. Note tstart is larger (older) that tstop for each record * 5. No gaps allowed: tstart must equal the previous records tstop * * Example: Duncan & Clague [1985] Pacific-Hotspot rotations: * * # Time in Ma, angles in degrees * # lon lat tstart tend ccw-angle * 165 85 150 100 24.0 * 284 36 100 74 15.0 * 265 22 74 65 7.5 * 253 17 65 42 14.0 * 285 68 42 0 34.0 * * ASCII point location file(s) must have the following format: * * 1. Any number of comment lines starting with # in first column * 2. Any number of blank lines (just carriage return, no spaces) * 3. For special header records, see -H * 4. Any number of data recordswhich each have the format: * lon lat age (or lat lon age, see -: option). age in Ma. * * Binary files cannot have header records, and data fields must all be * either single or double precision (see -bi option) */ #include "spotter.h" main (int argc, char **argv) { struct EULER *p; /* Pointer to array of stage poles */ int n_points; /* Number of data points read */ int n_chunk; /* Total length or array returned by libeuler functions */ int n_track; /* Number of points in a track segment */ int n_stages; /* Number of stage poles */ int n_segments; /* Number of path segments written out */ int n_skipped = 0; /* Number of points skipped because t < 0 */ int track_limit = 0; /* Flag for limiting the output track extent [no limit] */ int n_args; int n_files = 0, fno; int n_fields, n_expected_fields; int n_read = 0; int n_out; int i, j, k; /* Misc. counters */ BOOLEAN error = FALSE; /* Set to TRUE if arguments are inconsistent */ BOOLEAN flowline = FALSE; /* TRUE means we want flowlines, FALSE we want hotspot tracks */ BOOLEAN forward = FALSE; /* TRUE we go FROM hotspot to seamount, FALSE is reverse */ BOOLEAN make_path = FALSE; /* TRUE means create continuous path, FALSE works on discrete points */ BOOLEAN seg_files = FALSE; /* TRUE will write individual files for each segment */ BOOLEAN finite = FALSE; /* TRUE if stage pole file contains finite rotation poles instead */ BOOLEAN stage_id = FALSE; BOOLEAN first = TRUE; BOOLEAN done = FALSE; BOOLEAN nofile = TRUE; double d_km = 0.0; /* Step interval along calculated tracks */ double t_zero = 0.0; /* Current age in Ma */ double upper_age = 0.0; /* Extend oldest age back to this time, in Ma */ double *c; /* Array of track chunks returned by libeuler rutines */ double lon, lat; /* Seamounts location in decimal degrees */ double age; /* Age of seamount, in Ma */ double *in, *out; /* i/o arrays used by GMT */ double t_low, t_high; /* upper/lower age/stage to ouput for track segments */ char *euler_file = CNULL; /* Name pointer for file with stage poles */ char *file_stem = CNULL; /* Name pointer for file stem for individual track files */ char buffer[BUFSIZ]; /* Input buffer for reading data */ char type[50]; /* What kind of line (flowline or hotspot track) */ char dir[8]; /* From or To */ FILE *fp = NULL; /* File pointer for input data */ FILE *fpo = NULL; /* File pointer for output data */ PFI spot_func; /* Pointer to the requiered forth/back track function */ out = (double *)NULL; argc = GMT_begin (argc, argv); #ifdef DEBUG fprintf (stderr, "DEBUG mode\n"); #endif for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'H': case 'V': case ':': case '\0': error += GMT_get_common_args (argv[i], NULL, NULL, NULL, NULL); break; /* Supplemental parameters */ case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'A': /* Output only an age-limited segment of the track */ if (argv[i][2]) { /* Gave specific limits for all input points */ sscanf (&argv[i][2], "%lf/%lf", &t_low, &t_high); track_limit = 1; } else { /* Limits for each input point given in columns 4 and 5 */ track_limit = 2; } break; case 'C': /* Use finite rotation poles */ finite = TRUE; break; case 'D': /* Specify in which direction we should project */ switch (argv[i][2]) { case 'B': /* Go FROM hotspot TO seamount */ case 'b': forward = FALSE; break; case 'F': /* Go FROM seamount TO hotspot */ case 'f': forward = TRUE; break; default: error++; fprintf (stderr, "%s ERROR Option -D: Append b or f\n", GMT_program); break; } break; case 'L': /* Specify what kind of track to project */ switch (argv[i][2]) { case 'F': /* Calculate flowlines */ stage_id = TRUE; case 'f': flowline = TRUE; break; case 'B': /* Calculate hotspot tracks */ stage_id = TRUE; case 'b': flowline = FALSE; break; default: error++; fprintf (stderr, "%s ERROR Option -L: Append f or b\n", GMT_program); break; } d_km = (argv[i][3]) ? atof (&argv[i][3]) : -1.0; break; case 'E': /* File with stage poles */ euler_file = &argv[i][2]; break; case 'S': /* Set file stem for individual output files */ if (argv[i][2]) { file_stem = &argv[i][2]; seg_files = TRUE; } else { fprintf (stderr, "%s ERROR Option -S: Append a file stem\n", GMT_program); error++; } break; case 'T': /* Current age [0 Ma] */ t_zero = atof (&argv[i][2]); break; case 'N': /* Extend oldest stage back to this time [no extension] */ upper_age = atof (&argv[i][2]); break; default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else n_files++; } if (argc == 1 || GMT_quick) { fprintf (stderr, "%s %s - Forward and backward flowlines and hotspot tracks\n\n", GMT_program, SPOTTER_VERSION); fprintf (stderr, "usage: %s [infile(s)] -E [-A[young/old]] [-C] [-Df|b] [-H[]] [-Lf|b]\n", GMT_program); fprintf (stderr, "\t [-N] [-S] [-T] [-V] [-:] [-bi[s][]] [-bo[s][]]\n\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\tinfiles (in ASCII or binary) has 3 or more columns. If no file(s) is given, standard input is read.\n"); fprintf (stderr, "\tFirst 3 columns must have lon, lat (or lat, lon, see -:) and age (Ma)\n"); fprintf (stderr, "\t-E specifies the stage pole file to be used\n\n"); fprintf (stderr, "\tOPTIONS:\n\n"); fprintf (stderr, "\t-A Output tracks for ages (or stages, see -L) between young and old [Default is entire track]\n"); fprintf (stderr, "\t If no limit is given, then each seamount should have their limits in columns 4 and 5 instead.\n"); fprintf (stderr, "\t Only applicable in conjunction with the -L option.\n"); fprintf (stderr, "\t-C The file given with -E contains finite rotation poles [Default is stage poles]\n"); fprintf (stderr, "\t-Db move backward in time [Default]\n"); fprintf (stderr, "\t-Df move foreward in time\n"); GMT_explain_option ('H'); fprintf (stderr, "\t-Lb Compute hotspot tracks sampled every interval [Default projects single points].\n"); fprintf (stderr, "\t-Lf Compute flowline for seamounts of unknown but maximum age [Default projects single points]\n"); fprintf (stderr, "\t If no is given, the start/stop points for each stage are returned\n"); fprintf (stderr, "\t If B and F is used instead, stage id is returned as z-value [Default is predicted ages]\n"); fprintf (stderr, "\t-N extends earliest stage pole back to [no extension]\n"); fprintf (stderr, "\t-S writes tracks to individual files .# instead of to stdout (requires -L)\n"); fprintf (stderr, "\t-T sets the current age in Ma [0]\n"); GMT_explain_option ('V'); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf (stderr, "\t Default is 3 input columns\n"); GMT_explain_option ('o'); GMT_explain_option ('.'); fprintf (stderr, "\t Output produced by -Lf|t can be plotted with psxy using the -M option\n"); exit (EXIT_FAILURE); } if (GMT_io.binary[0] && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 3 + ((track_limit == 2) ? 2 : 0); if (GMT_io.binary[0] && GMT_io.ncol[0] < 3) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data (-bi) must have at least 3 columns\n", GMT_program); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ /* Load in the stage poles */ n_stages = spotter_init (euler_file, &p, flowline, finite, &upper_age); if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } spot_func = ((flowline + forward) == 1) ? spotter_forthtrack : spotter_backtrack; if (fabs (d_km) > SMALL) { /* User wants to interpolate tracks rather than project individual points */ make_path = TRUE; (flowline) ? sprintf (type, "Flowline") : sprintf (type, "Hotspot track"); (forward) ? sprintf (dir, "from") : sprintf (dir, "to"); } else if (track_limit > 0) { /* Limits on track set but track option not selected */ fprintf (stderr, "%s: GMT SYNTAX ERROR. -A requires -L\n", GMT_program); exit (EXIT_FAILURE); } fpo = GMT_stdout; n_out = (seg_files) ? 4 : 3; /* Append smt id number as 4th column when individual files are requested */ /* Read the seamount data from file or stdin */ n_points = n_segments = 0; if (n_files > 0) nofile = FALSE; else n_files = 1; n_args = (argc > 1) ? argc : 2; n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 3 + ((track_limit == 2) ? 2 : 0); for (fno = 1; !done && fno < n_args; fno++) { /* Loop over input files, if any */ if (!nofile && argv[fno][0] == '-') continue; if (nofile) { /* Just read standard input */ fp = GMT_stdin; done = TRUE; if (gmtdefs.verbose) fprintf (stderr, "%s: Reading from standard input\n", GMT_program); } else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]); continue; } if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]); if (gmtdefs.io_header) { for (i = 0; i < gmtdefs.n_header_recs; i++) { GMT_fgets (buffer, BUFSIZ, fp); if (first && !GMT_io.binary[1]) fprintf (GMT_stdout, "%s", buffer); } first = FALSE; } while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ n_read++; while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) { if (!make_path) GMT_write_segmentheader (GMT_stdout, n_expected_fields); n_fields = GMT_input (fp, &n_expected_fields, &in); n_read++; } if (GMT_io.status & GMT_IO_EOF) continue; if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d (skipped)\n", GMT_program, n_fields, n_expected_fields, n_read); continue; } if (in[2] < 0.0) { /* Negative ages are flags for points to be skipped */ n_skipped++; continue; } if (!out) out = (double *) GMT_memory (VNULL, (size_t)MAX(n_out, n_expected_fields), sizeof (double), GMT_program); lon = in[0] * D2R; lat = in[1] * D2R; if (track_limit) { if (track_limit == 2) t_low = in[3], t_high = in[4]; age = t_high; /* No point working more than necessary */ } else age = in[2]; if (age > upper_age) { /* Points older than oldest stage cannot be used */ fprintf (stderr, "%s: Seamount at line %d has age (%lg) > oldest stage (%lg) (skipped)\n", GMT_program, n_read, in[2], upper_age); n_skipped++; continue; } n_chunk = (*spot_func) (&lon, &lat, &age, 1, p, n_stages, d_km, t_zero, TRUE + stage_id, &c); if (make_path) { /* Asked for paths, now write out several multiple segment tracks */ if (seg_files) { sprintf (buffer, "%s.%d", file_stem, n_points); if ((fpo = GMT_fopen (buffer, GMT_io.w_mode)) == NULL) { fprintf (stderr, "%s: Error - cannot create file %s\n", GMT_program, buffer); exit (EXIT_FAILURE); } out[3] = (double)n_points; /* Put the seamount id number in 4th column */ } i = 0; n_track = irint (c[i++]); sprintf (GMT_io.segment_header, "> %s %s %lg %lg\n", type, dir, in[0], in[1]); GMT_write_segmentheader (fpo, n_out); for (j = 0; j < n_track; j++, i += 3) { out[2] = c[i+2]; if (track_limit && (out[2] < t_low || out[2] > t_high)) continue; out[0] = c[i] * R2D; out[1] = c[i+1] * R2D; GMT_output (fpo, n_out, out); } if (seg_files) GMT_fclose (fpo); GMT_free ((void *)c); } else { /* Just return the projected locations */ out[0] = lon * R2D; out[1] = lat * R2D; for (k = 2; k < n_expected_fields; k++) out[k] = in[k]; GMT_output (fpo, n_expected_fields, out); } n_points++; } if (fp != stdin) GMT_fclose (fp); } if (gmtdefs.verbose) { if (make_path) fprintf (stderr, "%s: %d segments written\n", GMT_program, n_points); else fprintf (stderr, "%s: %d points projected\n", GMT_program, n_points); } if (gmtdefs.verbose && n_skipped) fprintf (stderr, "%s: %d points skipped because age < 0\n", GMT_program, n_skipped); /* Clean up and exit */ GMT_end (argc, argv); }