JSOC_Documentation: proj/timed/apps/travel

00001 /*
00002  *  travel_times.c                       ~rick/hmi/timed
00003  *
00004  *  Calculate travel times for selected annuli in a tracked data cube using
00005  *    Gabor wavelet fitting and the method of Gizon & Birch
00006  *  The module expects a 3-dimensional real dataset and a pixel locator data
00007  *    file of a special format
00008  *
00009  *  Responsible:  Junwei Zhao & Rick Bogart       RBogart@solar.Stanford.EDU
00010  *
00011  *  Usage:
00012  *    travel_times [-v] in= pxloc= out= ...
00013  *
00014  *
00015  *  Arguments: (type        default         description)
00016  *  in  DataSet
00017  *      hmi_test.tdVtrack_synop[2010.07.08_04:00:00_TAI][][0.0][0.0]
00018  *                          Input dataset
00019  *          A series of records containing 3-d data cubes as one
00020  *          of their segments. It is assumed that all records in
00021  *          the dataset are from the same dataseries, or at least
00022  *          share a common segment structure
00023  *  pxloc   DataSet     hmi.tdpixlist[]
00024  *                  Annulus pixel locator data records
00025  *          A set of records specifying the different annuli to
00026  *          be used for analysis
00027  *  out DataSeries  hmi_test.tdVtimes_synop
00028  *                  Output data series
00029  *      copy    String      "+" Comma separated list of keys to
00030  *          propagate forward; if the list includes "+", all
00031  *          keys in the default list of propagation keys, plus
00032  *          all prime keys in the input records will be copied
00033  *
00034  *  Flags
00035  *
00036  *  Notes:
00037  *    All the work is done by the FORTRAN subroutines mainsub1_() & mainsub2_(),
00038  *  whose arguments are:
00039  *    (float) velo
00040  *    (float) wid
00041  *    (int) num_annuli
00042  *    (float) time[max(num_annuli)]
00043  *    (float) coef[5]
00044  *    (int) n_start
00045  *    (int) n_th
00046  *    float *input : input data values
00047  *    char *pix_locat_file (int *len_ploc_file)
00048  *      name of a FOTRAN unformatted file containing the following vectors
00049  *      of length n1 * n2 * num_annuli:
00050  *      (int) num_w
00051  *      (int) num_e
00052  *      (int) num_n
00053  *      (int) num_s
00054  *      (short) quad_w_x[num_w+1]
00055  *      (short) quad_w_y[num_w+1]
00056  *      (short) quad_e_x[num_e+1]
00057  *      (short) quad_e_y[num_e+1]
00058  *      (short) quad_n_x[num_n+1]
00059  *      (short) quad_n_y[num_n+1]
00060  *      (short) quad_s_x[num_s+1]
00061  *      (short) quad_s_y[num_s+1]
00062  *    float *out1
00063  *      array for the output of the do_fitting_() procedure
00064  *    float *out2
00065  *      array for the output of the gbtimes-2_() procedure
00066  *    Subsidiary FORTRAN code is in the following files:
00067  *    subroutine      location        called in
00068  *  C_CORRELATE sub_correlate_BLAS.f    MAINSUB1, MAINSUB2
00069  *  DO_FITTING  sub_do_fitting.f    MAINSUB1, MAINSUB2
00070  *  FILTER      sub_filter_HMI_ppline.f MAINSUB1, MAINSUB2
00071  *  GBTIMES02   MAINSUB1, MAINSUB2
00072  *  LMFIT       sub_lmfit.f     DO_FITTING
00073  *  SAXPY       library         MAINSUB1, MAINSUB2, DO_FITTING
00074  *  SCOPY       library         MAINSUB1, MAINSUB2, DO_FITTING
00075  *  SHIFT       library         MAINSUB1, MAINSUB2
00076  *  SSCAL       library         MAINSUB1, MAINSUB2, DO_FITTING
00077  *    Typical processing time on a single processor for a 256*256*640 cube
00078  *  ~
00079  *
00080  *  Bugs:
00081  *    Few parameter options; many of the parameters are read from pre-calculated
00082  *  data in the "pxloc: record. In particular, the input and output array
00083  *  sizes are hardcoded, but there is no checking of the actual input data.
00084  *      It is not clear that the module will work with any but a 512*512*640
00085  *  tracked cube
00086  *    Scaled output not properly supported; no setting of array parameter
00087  *    FORTRAN Code writes progress/status to stdout
00088  *    Missing input data are unconditionally set to 0
00089  *    The value written into CRLT_Mid is geocentric
00090  *    No documentation
00091  *
00092  *  Revision history is at the end of the file.
00093  *
00094  */
00095 
00096 #include <jsoc_main.h>
00097 #include "keystuff.c"
00098 #include "earth_ephem.c"
00099                               /*   module identifier  */
00100 char *module_name = "travel_times_loop_wGB";
00101 char *version_id = "0.9.2";
00102 
00103 ModuleArgs_t module_args[] = {
00104   {ARG_STRING, "in",
00105       "hmi_test.tdVtrack_synop[2010.07.08_04:00:00_TAI][][0.0][0.0]",
00106       "input record set"}, 
00107   {ARG_STRING, "pxloc", "hmi.tdpixlist[]", "pixel locater record set"}, 
00108   {ARG_STRING, "out", "hmi_test.tdVtimes_synop", "output series"}, 
00109   {ARG_STRING, "copy",  "+",
00110       "comma separated list of keys to propagate forward"},
00111   {ARG_FLAG, "v", "0", "verbose output"},      
00112   {}
00113 };
00114 
00115        /*  list of keywords to propagate (if possible) from input to output  */
00116 char *propagate[] = {"CarrRot", "CMLon", "LonHG", "LatHG", "LonCM", "MidTime",
00117     "CRVAL1", "CRVAL2", "CTYPE1", "CTYPE2", "CUNIT1", "CUNIT2",
00118     "T_START", "T_STOP", "LonSpan", "Coverage", "Duration", "MapProj", "MapPA",
00119     "Zonal Trk", "ZonalVel"};
00120 
00121 #define NUM_ANNULI_MAX  (20)
00122 
00123 extern void mainsub1_ (float *velo, float *wid, int *num_annuli, float *time,
00124     float *coef, int *n_start, int *n_th, float *input,
00125     char *pix_locat_file, int *len_ploc_file, float *out1, float *out2);
00126 extern void mainsub2_ (float *velo, float *wid, int *num_annuli, float *time,
00127     float *coef, int *n_start, int *n_th, float *input,
00128     char *pix_locat_file, int *len_ploc_file, float *out1, float *out2);
00129 
00130 int DoIt (void) {
00131   CmdParams_t *params = &cmdparams;
00132   DRMS_RecordSet_t *plds, *inds;
00133   DRMS_Record_t *ploc, *trck, *orec;
00134   DRMS_Segment_t *pseg, *dseg, *seg_gabor, *seg_gb;
00135   DRMS_Array_t *data_array, *res_gabor, *res_gb;
00136   FILE *timesfile, *fout;
00137   double mapscale_factor;
00138   float **ttimes;
00139   float *data, *input, *input2, *input3, *output, *output2, *out1, *out2;
00140   float *annmin, *annmax, *velo, *wid, *coef0, *coef1, *coef2, *coef3, *coef4;
00141   float coef[5];
00142   float newctr, newscale;
00143   long long n, ntot;
00144   int *num_annuli, *n_start;
00145   int trck_axis[2];
00146   int ann, i, j, k, *axes;
00147   int len_ploc_file, outok, ploc_version, status;
00148   int kstat, key_n, propct, propstd;
00149   int rec, rgnct, smpl, smpls;
00150   int need_ephem;
00151   char **pix_locat_file, **copykeylist;
00152   char *genstr, *errmsg;
00153   char plocfile[DRMS_MAXQUERYLEN];
00154   char source[DRMS_MAXQUERYLEN], recid[DRMS_MAXQUERYLEN];
00155   char module_ident[64], keyname[32];
00156 
00157   int keyct = sizeof (propagate) / sizeof (char *);
00158 
00159   char *plocds = strdup (params_get_str (params, "pxloc"));
00160   char *trckds = strdup (params_get_str (params, "in"));
00161   char *outser =  strdup (params_get_str (params, "out"));
00162   char *propagate_req = strdup (params_get_str (params, "copy"));
00163   int verbose = params_isflagset (params, "v");
00164 
00165   snprintf (module_ident, 64, "%s v %s", module_name, version_id);
00166   if (verbose) printf ("%s:\n", module_ident);
00167          /*  read info from correlation pixel-locator records to be used  */
00168   plds = drms_open_records (drms_env, plocds, &status);
00169   if (status) {
00170     fprintf (stderr, "Error: unable to open pixel target record set %s\n",
00171     plocds);
00172     return 1;
00173   }
00174   smpls = plds->n;
00175   if (verbose)
00176     printf ("  calculating travel time maps for %d annulus sets\n", smpls);
00177   if (smpls != 11)
00178     fprintf (stderr, "Warning: untested case of %d annulus sets\n", smpls);
00179   num_annuli = (int *)malloc (smpls * sizeof (int));
00180   annmin = (float *)malloc (smpls * sizeof (float));
00181   annmax = (float *)malloc (smpls * sizeof (float));
00182   velo = (float *)malloc (smpls * sizeof (float));
00183   wid = (float *)malloc (smpls * sizeof (float));
00184   n_start = (int *)malloc (smpls * sizeof (int));
00185   coef0 = (float *)malloc (smpls * sizeof (float));
00186   coef1 = (float *)malloc (smpls * sizeof (float));
00187   coef2 = (float *)malloc (smpls * sizeof (float));
00188   coef3 = (float *)malloc (smpls * sizeof (float));
00189   coef4 = (float *)malloc (smpls * sizeof (float));
00190   ttimes = (float **)malloc (smpls * sizeof (float *));
00191   pix_locat_file = (char **)malloc (smpls * sizeof (char **));
00192   for (smpl = 0; smpl < smpls; smpl++) {
00193     ploc = plds->records[smpl];
00194     drms_sprint_rec_query (source, ploc);
00195                    /*  check version of pixel locator record  */
00196     if (smpl) {
00197       if (ploc_version != drms_getkey_int (ploc, "Version", &status)) {
00198     if (ploc_version) {
00199       fprintf (stderr,
00200           "Warning: inconsistent versions of pixel locator records\n");
00201       fprintf (stderr,
00202           "         Version number will be set to 0 in output\n");
00203     }
00204     ploc_version = 0;
00205       }
00206     } else {
00207       ploc_version = drms_getkey_int (ploc, "Version", &status);
00208       if (status) {
00209     fprintf (stderr,
00210         "Warning: unknown version of pixel locator record(s)\n");
00211         ploc_version = 0;
00212       }
00213     }
00214                      /*  get additional key info from record  */
00215     num_annuli[smpl] = drms_getkey_int (ploc, "Annuli", &status);
00216     if (num_annuli[smpl] > NUM_ANNULI_MAX) {
00217       fprintf (stderr,
00218           "Warning: num_annuli exceeds fixed mem alloc in Fortran code\n");
00219       fprintf (stderr, "         reduced to %d\n", NUM_ANNULI_MAX);
00220       num_annuli[smpl] = NUM_ANNULI_MAX;
00221     }
00222     annmin[smpl] = drms_getkey_float (ploc, "RadInner", &status);
00223     annmax[smpl] = drms_getkey_float (ploc, "RadOuter", &status);
00224     velo[smpl] = drms_getkey_float (ploc, "PhaseVel", &status);
00225     wid[smpl] = drms_getkey_float (ploc, "FiltWid", &status);
00226     coef0[smpl] = drms_getkey_float (ploc, "AmplInit", &status);
00227     coef1[smpl] = drms_getkey_float (ploc, "DFreqIni", &status);
00228     coef2[smpl] = drms_getkey_float (ploc, "GrpTTIni", &status);
00229     coef3[smpl] = drms_getkey_float (ploc, "FreqInit", &status);
00230     coef4[smpl] = drms_getkey_float (ploc, "PhsTTIni", &status);
00231     n_start[smpl] = drms_getkey_float (ploc, "NStart", &status);
00232     ttimes[smpl] = (float *)calloc (NUM_ANNULI_MAX, sizeof (float));
00233     pseg = drms_segment_lookup (ploc, "times");
00234     if (!pseg) {
00235       fprintf (stderr,
00236           "Error: record %s\n  does not contain required segment \"times\"\n",
00237       source);
00238       drms_close_records (plds, DRMS_FREE_RECORD);
00239       return 1;
00240     }
00241     drms_segment_filename (pseg, plocfile);
00242     timesfile = fopen (plocfile, "r");
00243     for (n = 0; n < num_annuli[smpl]; n++)
00244       fscanf (timesfile, "%f", &(ttimes[smpl][n]));
00245           /*  should not be necessary, calloc should take care of it  */
00246     for (; n < NUM_ANNULI_MAX; n++) ttimes[smpl][n] = 0.0;
00247     fclose (timesfile);
00248                                   /*  read times  */
00249     pseg = drms_segment_lookup (ploc, "pxloc.fortran");
00250     if (!pseg) {
00251       fprintf (stderr,
00252           "Error: record %s\n  does not contain required segment \"pxloc.fortran\"\n",
00253       source);
00254       drms_close_records (plds, DRMS_FREE_RECORD);
00255       return 1;
00256     }
00257     drms_segment_filename (pseg, plocfile);
00258     pix_locat_file[smpl]= (char *)malloc (strlen (plocfile) + 1);
00259     strncpy (pix_locat_file[smpl], plocfile, strlen (plocfile));
00260     pix_locat_file[smpl][strlen (plocfile)] = '\0';
00261   }
00262                            /*  allocate the output array  */
00263   axes = (int *)malloc (4 * sizeof (int));
00264   axes[0] = axes[1] = 256;
00265   axes[2] = 4;
00266   axes[3] = smpls;
00267   output = (float *)malloc (smpls * 4 * 256 * 256 * sizeof (float));
00268   output2 = (float *)malloc (smpls * 4 * 256 * 256 * sizeof (float));
00269   out1 = (float *)malloc (4 * 256 * 256 *sizeof(float));
00270   out2 = (float *)malloc (4 * 256 * 256 *sizeof(float));
00271   res_gabor = drms_array_create (DRMS_TYPE_FLOAT, 4, axes, output, &status);
00272   res_gb = drms_array_create (DRMS_TYPE_FLOAT, 4, axes, output2, &status);
00273   drms_close_records (plds, DRMS_FREE_RECORD);
00274 
00275   input2 = (float *)malloc (256 * 256 * 640 * sizeof(float));
00276   input3 = (float *)malloc (256 * 256 * 640 * sizeof(float));
00277                              /*  open input data set  */
00278   inds = drms_open_records (drms_env, trckds, &status);
00279   if (status) {
00280     fprintf (stderr, "Error: unable to open input record set %s\n", trckds);
00281     drms_free_array (res_gabor);
00282     drms_free_array (res_gb);
00283     return 1;
00284   }
00285   rgnct = inds->n;
00286   if (verbose)
00287     printf ("  in each of %d regions\n", rgnct);
00288 
00289   outok = 0;
00290   for (rec = 0; rec < rgnct; rec++) {
00291     trck = inds->records[rec];
00292     dseg = drms_segment_lookup (trck, "V");
00293     if (!dseg) {
00294       fprintf (stderr, "Error: no segment \"V\" in record %d; skipped\n", rec);
00295       continue;
00296     }
00297     for (i = 0; i < 2; i++) {
00298       trck_axis[i] = dseg->axis[i];
00299       if (i) {
00300         if (trck_axis[i] - trck_axis[i-1]) {
00301       fprintf (stderr,
00302           "Error: code as implemented requires square tracked regions\n");
00303       drms_free_array (res_gabor);
00304       drms_free_array (res_gb);
00305       drms_close_records (inds, DRMS_FREE_RECORD);
00306       return 1;
00307     }
00308       }
00309     }
00310     mapscale_factor = 256.0 / (double)trck_axis[0];
00311                         /*  create the output record  */
00312     orec = drms_create_record (drms_env, outser, DRMS_PERMANENT, &status);
00313     if (!orec) {
00314       fprintf (stderr, "Error: unable to create record %d in series %s\n", rec,
00315       outser);
00316       fprintf (stderr, "       remainder of processing skipped\n");
00317       drms_close_records (inds, DRMS_FREE_RECORD);
00318       return 0;
00319     }
00320      /*  check the output data series struct (only needs to be done once  */
00321     seg_gabor = drms_segment_lookup (orec, "gabor");
00322     seg_gb = drms_segment_lookup (orec, "GB");
00323     if (!outok) {
00324                       /*  check for the appropriate segments  */
00325       if (!seg_gabor || !seg_gb) {
00326     fprintf (stderr, "Error: output data series %s\n", outser);
00327     fprintf (stderr, "       lacks segment \"GB\" and/or \"gabor\"\n");
00328     drms_close_records (inds, DRMS_FREE_RECORD);
00329     return 1;
00330 
00331       }
00332     /*  check structure of segments: in particular naxis[3] must match smpls  */
00333       outok = 1;
00334     }
00335     res_gabor->bscale = res_gb->bscale = 1.0;
00336     res_gabor->bzero = res_gb->bzero = 0.0;
00337                      /*  copy designated keywords from input  */
00338     propct = construct_stringlist (propagate_req, ',', &copykeylist);
00339     kstat = 0;
00340     propstd = 0;
00341                     /*  copy specifically requested keys  */
00342     for (key_n = 0; key_n < propct; key_n++) {
00343       if (strcmp (copykeylist[key_n], "+"))
00344         kstat += check_and_copy_key (orec, trck, copykeylist[key_n]);
00345       else propstd = 1;
00346     }
00347     if (propstd) {
00348                               /*  copy standard keys  */
00349       kstat += propagate_keys (orec, trck, propagate, keyct);
00350                        /*  and any additional prime keys  */
00351       kstat += copy_prime_keys (orec, trck);
00352     }
00353     if (kstat) {
00354       fprintf (stderr, "Error writing key value(s) to record in series %s:\n",
00355       outser);
00356       fprintf (stderr, "      series may not have appropriate structure;\n");
00357       fprintf (stderr, "      record %d skipped\n", rec);
00358       drms_close_record (orec, DRMS_FREE_RECORD);
00359       continue;
00360     }
00361     drms_sprint_rec_query (source, trck);
00362             /*  to write CRLT_Mid to output, require MidTime as well  */
00363     need_ephem = (drms_keyword_lookup (orec, "CRLT_Mid", 1) != NULL);
00364     if (need_ephem)
00365       need_ephem = (drms_keyword_lookup (orec, "MidTime", 1) != NULL);
00366     if (need_ephem) {
00367       double rsun, clat = 0.0 / 0.0, clon, vr, vn, vw;
00368       TIME tmid =  drms_getkey_time (orec, "MidTime", &status);
00369       earth_ephemeris (tmid, &rsun, &clat, &clon, &vr, &vn, &vw);
00370       kstat += check_and_set_key_float (orec, "CRLT_Mid", clat);
00371     }
00372                                 /*  set WCS keys  */
00373     kstat += check_and_set_key_int (orec, "WCSAXES", 2);
00374     newctr = drms_getkey_float (trck, "CRPIX1", &status);
00375     newctr = mapscale_factor * (newctr - 0.5) + 0.5;
00376     kstat += check_and_set_key_float (orec, "CRPIX1", newctr);
00377     newctr = drms_getkey_float (trck, "CRPIX2", &status);
00378     newctr = mapscale_factor * (newctr - 0.5) + 0.5;
00379     kstat += check_and_set_key_float (orec, "CRPIX2", newctr);
00380     newscale = drms_getkey_float (trck, "CDELT1", &status);
00381     newscale /= mapscale_factor;
00382     kstat += check_and_set_key_float (orec, "CDELT1", newscale);
00383     newscale = drms_getkey_float (trck, "CDELT2", &status);
00384     newscale /= mapscale_factor;
00385     kstat += check_and_set_key_float (orec, "CDELT2", newscale);
00386     genstr = drms_getkey_string (trck, "WCSNAME", &status);
00387     if (!status) {
00388           /*  strip off the trailing '/Time' from a tracked cube WCSNAME  */
00389       strtok (genstr, "/");
00390       kstat += check_and_set_key_str (orec, "WCSNAME", genstr);
00391     }
00392                          /*  set segment descriptor keys  */
00393     kstat += check_and_set_key_str (orec, "TTDiff_1",
00394     "mean (outgoing,ingoing)");
00395     kstat += check_and_set_key_str (orec, "TTDiff_2", "outgoing-ingoing");
00396     kstat += check_and_set_key_str (orec, "TTDiff_3", "west-east");
00397     kstat += check_and_set_key_str (orec, "TTDiff_4", "north-south");
00398                          /*  update the MapScale key  */
00399     newscale = drms_getkey_float (trck, "MapScale", &status);
00400     if (!status) {
00401       newscale /= mapscale_factor;
00402       kstat += check_and_set_key_float (orec, "MapScale", newscale);
00403     }
00404                          /*  set annulus descriptor keys  */
00405     for (smpl = 0; smpl < smpls; smpl++) {
00406       sprintf (keyname, "AnnMin%02d", smpl + 1);
00407       kstat += check_and_set_key_float (orec, keyname, annmin[smpl]);
00408       sprintf (keyname, "AnnMax%02d", smpl + 1);
00409       kstat += check_and_set_key_float (orec, keyname, annmax[smpl]);
00410       sprintf (keyname, "PhsVel%02d", smpl + 1);
00411       kstat += check_and_set_key_float (orec, keyname, velo[smpl]);
00412       sprintf (keyname, "FltWid%02d", smpl + 1);
00413       kstat += check_and_set_key_float (orec, keyname, wid[smpl]);
00414     }
00415     if (kstat) {
00416       fprintf (stderr, "Error writing key value(s) to record in series %s:\n",
00417       outser);
00418       fprintf (stderr, "      series may not have appropriate structure;\n");
00419       fprintf (stderr, "      record %d skipped\n", rec);
00420       drms_close_record (orec, DRMS_FREE_RECORD);
00421       continue;
00422     }
00423                              /*  read the input data  */
00424     data_array = drms_segment_read (dseg, DRMS_TYPE_FLOAT, &status);
00425     input = (float *)data_array->data;
00426     if (!input) {
00427       fprintf (stderr, "Error: unable to read input data record-segment %d\n",
00428           rec);
00429       drms_close_record (orec, DRMS_FREE_RECORD);
00430       continue;
00431     }
00432     for (i=0; i<640; i++) 
00433       for (j=0; j<256; j++)
00434     for (k=0; k<256; k++)
00435       input2[i*256*256 + j*256 + k] =
00436           0.25*(input[i*512*512 + (2*j*512) + (2*k)] + 
00437           input[i*512*512 + (2*j)*512 + (2*k+1)] +
00438           input[i*512*512 + (2*j+1)*512 + (2*k)] + 
00439           input[i*512*512 + (2*j+1)*512 + (2*k+1)]);
00440             /*  zero out any missing data; Fourier transforms object  */
00441     ntot = 1;
00442     for (i = 0; i < data_array->naxis; i++) ntot *= data_array->axis[i];
00443     for (n = 0; n < ntot; n++) if (isnan (input[n])) input[n] = 0.0;
00444 
00445     for (ann=0; ann<3; ann++) {
00446       int nth = ann + 1;
00447       if (verbose) printf ("begin annulus %d\n", ann);
00448              /*  string lengths required for C/FORTRAN interface  */
00449       len_ploc_file = strlen (pix_locat_file[ann]);
00450                    /*  Call Fortran Code subroutine MAIN_SUB  */
00451       coef[0] = coef0[ann];
00452       coef[1] = coef1[ann];
00453       coef[2] = coef2[ann];
00454       coef[3] = coef3[ann];
00455       coef[4] = coef4[ann];
00456       mainsub1_ (&velo[ann], &wid[ann], &num_annuli[ann], ttimes[ann], coef,
00457       &n_start[ann], &nth, input, pix_locat_file[ann],
00458       &len_ploc_file, out1, out2);
00459 
00460       for (j=0; j<4*256*256; j++) {
00461         output[ann*4*256*256+j] = out1[j]; 
00462         output2[ann*4*256*256+j] = out2[j];
00463       }
00464          /*  need to reread segment because input array may have changed  */
00465       drms_free_array (data_array);
00466       data_array = drms_segment_read (dseg, DRMS_TYPE_FLOAT, &status);
00467       input = (float *)data_array->data;
00468       if (!input) {
00469     fprintf (stderr, "Error: unable to read input data record-segment %d\n",
00470             rec);
00471     drms_close_record (orec, DRMS_FREE_RECORD);
00472     ann = 11;
00473         continue;
00474       }
00475       for (n = 0; n < ntot; n++) if (isnan (input[n])) input[n] = 0.0;
00476     }
00477 
00478     ntot = 256 * 256 * 640;
00479     for (n = 0; n < ntot; n++) if (isnan (input2[n])) input2[n] = 0.0;
00480     for (; ann<11; ann++) {
00481       int nth = ann + 1;
00482       if (verbose) printf ("begin annulus %d\n", ann);
00483       for (n = 0; n < ntot; n++) input3[n] = input2[n];
00484              /*  string lengths required for C/FORTRAN interface  */
00485       len_ploc_file = strlen (pix_locat_file[ann]);
00486       coef[0] = coef0[ann];
00487       coef[1] = coef1[ann];
00488       coef[2] = coef2[ann];
00489       coef[3] = coef3[ann];
00490       coef[4] = coef4[ann];
00491       mainsub2_ (&velo[ann], &wid[ann], &num_annuli[ann], ttimes[ann], coef,
00492       &n_start[ann], &nth, input3, pix_locat_file[ann],
00493       &len_ploc_file, out1, out2);
00494       for (j=0; j<4*256*256; j++) {
00495         output[ann*4*256*256 + j] = out1[j]; 
00496         output2[ann*4*256*256 + j] = out2[j];
00497       }
00498     }
00499     status = drms_segment_write (seg_gabor, res_gabor, 0);
00500     if (status) {
00501       drms_sprint_rec_query (recid, orec);
00502       fprintf (stderr, "Error writing data to record %s\n", recid);
00503       fprintf (stderr, "      drms_segment_write returned %d\n", status);
00504       fprintf (stderr, "      series may not have appropriate structure\n");
00505       fprintf (stderr, "      record %d skipped\n", rec);
00506       drms_close_record (orec, DRMS_FREE_RECORD);
00507       continue;
00508     }
00509     status = drms_segment_write (seg_gb, res_gb, 0);
00510     if (status) {
00511 
00512       drms_sprint_rec_query (recid, orec);
00513       fprintf (stderr, "Error writing data to record %s\n", recid);
00514       fprintf (stderr, "      series may not have appropriate structure\n");
00515       fprintf (stderr, "      record %d skipped\n", rec);
00516       drms_close_record (orec, DRMS_FREE_RECORD);
00517       continue;
00518     }
00519                             /*  set a few final keys  */
00520     kstat += check_and_set_key_str (orec, "Module", module_ident);
00521     kstat += check_and_set_key_str (orec, "BLD_VERS", jsoc_version);
00522     kstat += check_and_set_key_str (orec, "Source", source);
00523     kstat += check_and_set_key_str (orec, "Input", trckds);
00524     kstat += check_and_set_key_str (orec, "PixLocate", plocds);
00525     kstat += check_and_set_key_int (orec, "PxLocVer", ploc_version);
00526     kstat += check_and_set_key_time (orec, "Created", CURRENT_SYSTEM_TIME);
00527     if (kstat) {
00528       fprintf (stderr, "Error writing key value(s) to record in series %s:\n",
00529       outser);
00530       fprintf (stderr, "      series may not have appropriate structure;\n");
00531       fprintf (stderr, "      record %d skipped\n", rec);
00532       drms_close_record (orec, DRMS_FREE_RECORD);
00533       continue;
00534     }
00535 
00536     drms_close_record (orec, DRMS_INSERT_RECORD);
00537   }
00538   drms_close_records (inds, DRMS_FREE_RECORD);
00539   return 0;
00540 }
00541 /*
00542  *  Revision History (all mods by Rick Bogart unless otherwise noted)
00543  *
00544  *  10.09.24    created this file as prototype
00545  *  10.10.08    copied from working and tested version file
00546  *  /tmp00/junwei/HMI_script/travel_time/travel_times_loops_wGB.c
00547  *      added comments; minor mod to dup returned string pointers;
00548  *      incorporated get_data function in file from original file
00549  *  /tmp00/junwei/HMI_script/travel_time/get_data.c
00550  *  v 0.5 frozen 2010.10.11
00551  *  10.10.11    removed extraneous (commented) code
00552  *      replaced hard-coded parameter file names with hard-coded
00553  *  parameter values
00554  *  10.11.15    added verbose argument and some diagnostic output; added
00555  *  argument for and processing of DRMS data records for correlation pixel
00556  *  sets; input read from DRMS data records, with option of processing
00557  *  multiple records
00558  *  v 0.6 frozen 2011.01.06
00559  *  11.01.28    stripped old test code, debugging lines; changed output from
00560  *  named files to DRMS series; added some error checking on input and
00561  *  output; removed get_data, need for cfitsio; fixed possible indexing
00562  *  error; removed included fftw3 source from sub_filter_HMI_ppline.f;
00563  *  changed output to use drms utilities rather than private fits code;
00564  *  added several standard output keys; incorporated updated version of
00565  *  sub_GB_fitting_2002.f; removed debug prints from fortran routines
00566  *  11.02.25    fixed bug in setting of filename for fortran
00567  *  11.03.07    free data_array before rereading; added debug statements
00568  *  11.05.27    zero out any missing input data
00569  *  v 0.7 frozen 2011.06.08
00570  *  11.06.10    added propagation of version number of pixel locator records
00571  *  11.07.18    added setting of CRLT_Mid; added WCS keys plus numerous others
00572  *  to default propagation list
00573  *  v 0.8 frozen 2011.07,18
00574  *  11.12.03    removed setting of MapScale key, pending availability from
00575  *  pxloc data
00576  *  12.01.13    fixed setting of MapScale key; added setting of annuli info
00577  *  keys
00578  *  v 0.9 frozen 2012.01.14
00579  *  v 0.91 frozen 2012.08.04
00580  *  12.09.28    added option of adding to or overriding standard list of
00581  *  keywords from input to output
00582  *  13.03.12    corrected setting of  CRPIX1,2 and CDELT1,2 to values consistent
00583  *  with pixel locator file (supposedly); added MapProj and MapPA to the
00584  *  default list of propagated keywords; added modified propagation of
00585  *  WCSNAME value
00586  */
proj/timed/apps/travel_times.c
