JSOC_Documentation: proj/rings/apps/rdvinv.c Source File

00001 /*****************************************************************
00002  *  rdvinv.c                        ~rick/src/rdvinv
00003  *
00004  *  Responsible: Charles Baldner, Rick Bogart, Sarbani Basu
00005  *                          RBogart@spd.aas.org
00006  *
00007  *  Velocity inversion code for ring diagram power spectra fits
00008  *
00009  *  Parameters: (type   default         description)
00010  *  in  string  -       Input data set descriptor or filename
00011  *  out string  fort.10.hmixy   Output filename
00012  *  seg string  fit.out     Input data segment name
00013  *  cr  int NS      Carrington rotation common to dataset
00014  *  clon    float   NS      CM Longitude common to dataset
00015  *  lon float   NS      Carrington longitude common to dataset
00016  *  lat float   NS      Heliographic latitude common to dataset
00017  *  amu double  0.005       Error trade-off parameter
00018  *  ob  double  1.0     Lower frequency limit (mHz)
00019  *  oe  double  5.2     Upper frequency limit (mHz)
00020  *  rb  double  0.97        Lower radius limit
00021  *  re  double  1.00        Upper radius limit
00022  *  num int 40      Number of target inversion points
00023  *  kernel  string  -       Data record or pathname of file
00024  *              containing mode kernels for inversions
00025  *  ave string  Not Specified   Output file for averaging kernels
00026  *              (if not set, kernels are not written)
00027  *  coef    string  Not Specified   Output file for ?? coefficients
00028  *              (if not set, coefficients are not written)
00029  *
00030  *  Flags
00031  *  -v  run verbose
00032  *  
00033  *  Notes:
00034  *    Output is specified by argment 'out'.  If 'out' is a drms data series
00035  *  that can be opened, the inversions will be  saved to the appropriate
00036  *  drms records; otherwise, 'out' will be used as the base filename for
00037  *  the inversion outputs.
00038  *
00039  *  Bugs:
00040  *    If the input data is not in a drms series, an output series can't be
00041  *  created - rdvinv will exit with an error message.
00042  *    Argument 'out', when a file name, cannot contain any path information.
00043  *    When writing output to a file, rather than a drms record set, only one
00044  *  set of files can be specified; if the input contains more than one
00045  *  data record, each inversion will overwrite the last.
00046  *    There is no verification that the requested output segments have the
00047  *  correct protocol; probably doesn't matter anyway
00048  *    The constant value of amu is written needlessly for each target location
00049  *  for backward compatability
00050  *    The reading of the kernel file is still in the FORTRAN subroutine
00051  *
00052  *  Revision history is at end of file
00053  ******************************************************************/
00054 
00055 
00056 #include <stdio.h>
00057 #include <string.h>
00058 #include <math.h>
00059 #include <jsoc_main.h>
00060 #include "keystuff.c"
00061 #include "old_rdutil.c"
00062 #include "ola_xy.c"
00063 
00064 /* prototypes */
00065 
00066 extern void ola_(double *, int *, double *, double *, double *,
00067       double *, double *, int *, char *, double *, double *, double *,
00068       double *, double *, double *, double *, double *, double *,
00069       char *, char *, int *, int *, int *, double *, double *,
00070       int *, double *, double *, double *, int *, int, int, int);
00071 
00072 char *module_name = "rdvinv";
00073 char *version_id = "0.91";
00074 
00075 ModuleArgs_t module_args[] = {
00076   {ARG_STRING, "in", "", "Input data series or recordset"},
00077   {ARG_STRING, "seg", "fit.out", "Input data segment name"},
00078   {ARG_STRING, "uxseg", "Ux", "Output data segment name for Ux inversions"},
00079   {ARG_STRING, "uyseg", "Uy", "Output data segment name for Uy inversions"},
00080   {ARG_INT,    "cr", "Not Specified", "Carrington rotation for all regions"},
00081   {ARG_FLOAT,  "clon", "Not Specified",
00082       "Carrington longitude of CM for all regions"},
00083   {ARG_FLOAT,  "lon", "Not Specified", "Carrington longitude of all regions"},
00084   {ARG_FLOAT,  "lat", "Not Specified", "Carrington latitude of all regions"},
00085   {ARG_DOUBLE, "amu", "0.005", "Error trade-off parameter"},
00086   {ARG_DOUBLE, "ob", "1.0", "Lower frequency limit (mHz)"},
00087   {ARG_DOUBLE, "oe", "5.2", "Upper frequency limit (mHz)"},
00088   {ARG_DOUBLE, "rb", "0.97", "Lower radius limit"},
00089   {ARG_DOUBLE, "re", "1.00", "Upper radius limit"},
00090   {ARG_INT,    "num", "40", "Number of target inversion points"},
00091   {ARG_STRING, "out", "fort.10.hmixy", "Output filename"},
00092   {ARG_STRING, "kernel", "", ""},
00093   {ARG_STRING, "ave",  "Not Specified",
00094       "output file for averaging kernels (if not set, kernels are not written)"},
00095   {ARG_STRING, "coef", "Not Specified", ""},
00096   {ARG_FLAG,   "v",    "", "run in verbose mode"},
00097   {ARG_END}
00098 };
00099 
00100        /*  list of keywords to propagate (if possible) from input to output  */
00101 char *propagate[] = {"CarrTime", "CarrRot", "CMLon", "LonHG", "LatHG", "LonCM",
00102     "MidTime", "Duration", "MapProj", "MapScale", "Map_PA", "Width", "Height",
00103     "Size", "Cadence", "ZonalTrk", "ZonalVel", "MeridTrk", "MeridVel", "MAI"};
00104 
00105 int process_record (const char *filename, int drms_output, char *outfilex,
00106     char *outfiley, char *kernel, char *kername, char *ave, char *coef,
00107     int qave, int qcoef, int verbose, double ob, double oe, int num,
00108     double rb, double re, double amu, char *sourcerec, char *codever) {
00109   FILE *infile = fopen (filename, "r");
00110   FILE *filex = fopen (outfilex, "w");
00111   FILE *filey = fopen (outfiley, "w");
00112   double *rcgx, *rcgy, *quartx, *quarty, *solnx, *solny, *errx, *erry;
00113   double *l, *f, *ef, *ux, *eux, *uy, *euy, *rtrg;
00114   int *n, *mask;
00115   int npts, i, j, status;
00116   int lenkern, lenave, lencoef;
00117   status = read_fit_v (infile, &npts, &n, &l, &f, &ef, &ux, &eux, &uy, &euy);
00118   fclose (infile);
00119   if(status)    {
00120     fprintf(stderr, "File %s could not be read\n", filename);
00121     return 1;
00122   }
00123 
00124   rtrg = (double *)malloc (num * sizeof (double));
00125   rcgx = (double *)malloc (num * sizeof (double));
00126   rcgy = (double *)malloc (num * sizeof (double));
00127   quartx = (double *)malloc (3 * num * sizeof (double));
00128   quarty = (double *)malloc (3 * num * sizeof (double));
00129   solnx = (double *)malloc (num * sizeof (double));
00130   solny = (double *)malloc (num * sizeof (double));
00131   errx = (double *)malloc (num * sizeof (double));
00132   erry = (double *)malloc (num * sizeof (double));
00133   mask = (int *)malloc (npts * sizeof(int));
00134   interp_vel (n, l, f, ef, ux, eux, uy, euy, npts);
00135   autoweed_vel (n, l, ux, uy, mask, npts);
00136   j = 0;
00137   for (i=0; i<npts; i++) {
00138     if (mask[i]) {
00139       n[j] = n[i];
00140       l[j] = l[i];
00141       f[j] = f[i];
00142       ef[j] = ef[i];
00143       ux[j] = ux[i];
00144       eux[j] = eux[i];
00145       uy[j] = uy[i];
00146       euy[j] = euy[i];
00147       j++;
00148     }
00149   }
00150   npts = j;
00151 
00152   fprintf (filex, "# Flow inversion of %s\n", sourcerec);
00153   fprintf (filey, "# Flow inversion of %s\n", sourcerec);
00154   fprintf (filex, "# against kernel %s\n", kername);
00155   fprintf (filey, "# against kernel %s\n", kername);
00156   fprintf (filex, "# %s\n", codever);
00157   fprintf (filey, "# %s\n", codever);
00158   fprintf (filex, "# amu = %.5f ob = %.3f mHz, oe = %.3f mHz\n", amu, ob, oe);
00159   fprintf (filey, "# amu = %.5f ob = %.3f mHz, oe = %.3f mHz\n", amu, ob, oe);
00160   fprintf (filex, "# rb = %.3f, re= %.3f, num = %d\n",rb, re, num);
00161   fprintf (filey, "# rb = %.3f, re= %.3f, num = %d\n",rb, re, num);
00162   fprintf (filex,
00163       "# Rtrg    amu    CG of av. kernel, quartiles                soln          err\n");
00164   fprintf (filex, "#                           1       2       3     3-1\n");
00165   fprintf (filey,
00166       "# Rtrg    amu    CG of av. kernel, quartiles                soln          err\n");
00167   fprintf (filey, "#                           1       2       3     3-1\n");
00168 
00169   lenkern = strlen (kernel);
00170   lenave = (qave) ? strlen (ave) : 0;
00171   lencoef = (qcoef) ? strlen (coef) : 0;
00172 
00173   ola_ (l, n, f, ux, eux, uy, euy, &npts, kernel, rtrg, rcgx, rcgy,
00174       quartx, quarty, solnx, solny, errx, erry, ave, coef, &qave, &qcoef,
00175       &verbose, &ob, &oe, &num, &rb, &re, &amu, &status,
00176       lenkern, lenave, lencoef);
00177   if (status) return status;
00178 
00179   for (i = 0; i < num; i++)
00180     fprintf (filex, "%7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %12.5e %12.5e\n",
00181     rtrg[i], amu, rcgx[i], quartx[2 + 3*i], quartx[1 + 3*i], quartx[3*i],
00182     fabs (quartx[3*i] - quartx[2+ 3*i]), solnx[i], errx[i]);
00183   for (i = 0; i < num; i++)
00184     fprintf (filey, "%7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %7.5f %12.5e %12.5e\n",
00185     rtrg[i], amu, rcgy[i], quarty[2 + 3*i], quarty[1 + 3*i], quarty[3*i],
00186     fabs (quarty[3*i] - quarty[2+ 3*i]), solny[i], erry[i]);
00187 
00188   free (n);
00189   free (mask);
00190   free (l);
00191   free (f);
00192   free (ef);
00193   free (ux);
00194   free (eux);
00195   free (uy);
00196   free (euy);
00197 
00198   free (rtrg);
00199   free (rcgx);
00200   free (rcgy);
00201   free (quartx);
00202   free (quarty);
00203   free (solnx);
00204   free (solny);
00205   free (errx);
00206   free (erry);
00207   fclose (filex);
00208   fclose (filey);
00209 
00210   return 0;
00211 }
00212          /*  these functions should probably be added to keystuff.c  */
00213                         /*  they are also in rdsinv_v07  */
00214 int key_is_prime (DRMS_Record_t *rec, const char *key) {
00215   DRMS_Keyword_t *pkey, *keywd;
00216   int n, npct = rec->seriesinfo->pidx_num;
00217 
00218   if (!(keywd = drms_keyword_lookup (rec, key, 1))) return 0;
00219   if (npct) {
00220     for (n = 0; n < npct; n++) {
00221       pkey = rec->seriesinfo->pidx_keywords[n];
00222       if (pkey->info->recscope > 1) pkey = drms_keyword_slotfromindex (pkey);
00223       if (strcmp (key, pkey->info->name)) continue;
00224       return 1;
00225     }
00226   }
00227   return 0;
00228 }
00229 
00230 int unique_key_values (const char *dsqry, const char *key, DRMS_Record_t *rec) {
00231   int status, values;
00232   DRMS_Array_t *value;
00233   DRMS_Keyword_t *pkey, *keywd;
00234 
00235   if (!(keywd = drms_keyword_lookup (rec, key, 1))) return 0;
00236   value = drms_record_getvector (drms_env, dsqry, key,
00237       DRMS_TYPE_FLOAT, 1, &status);
00238   if (!value) return 0;
00239   values = value->axis[1];
00240   drms_free_array (value);
00241   return values;
00242 }
00243 
00244 char *printcrcl_loc (int cr, double cl, double lnhg, double lncm, double lat) {
00245   static char fnam[32];
00246   double dlon = lncm;
00247 
00248   if (cr > 0) {
00249     if (isnan (cl)) {
00250       if (isnan (lncm)) {
00251                        /*  construct a filename CR_lonHG+lat  */
00252     sprintf (fnam, "%04d_%04.1f%+05.1f", cr, lnhg, lat);
00253       } else {
00254                       /*  construct a filename CR_lonEWlatNS  */
00255     if (dlon >= 0) {
00256       if (lat >= 0) sprintf (fnam, "%04d_%04.1fW%04.1fN", cr, dlon, lat);
00257       else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fW%04.1fN",
00258           cr, dlon, -lat);
00259           else sprintf (fnam, "%04d_%04.1fW%04.1fS", cr, dlon, -lat);
00260     } else if (dlon > -0.05) {
00261       if (lat >= 0) sprintf (fnam, "%04d_%04.1fW%04.1fN", cr, -dlon, lat);
00262       else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fW%04.1fN",
00263           cr, -dlon, -lat);
00264           else sprintf (fnam, "%04d_%04.1fW%04.1fS", cr, -dlon, -lat);
00265     } else {
00266       if (lat >= 0) sprintf (fnam, "%04d_%04.1fE%04.1fN", cr, -dlon, lat);
00267       else if (lat > -0.05) sprintf (fnam, "%04d_%04.1fE%04.1fN",
00268           cr, -dlon, -lat);
00269           else sprintf (fnam, "%04d_%04.1fE%04.1fS", cr, -dlon, -lat);
00270     }
00271       }
00272       return fnam;
00273     }
00274                    /*  construct a filename CR:CL_lonEWlatNS  */
00275     while (cl < 0) {
00276       cl += 360;
00277       cr++;
00278     } while (cl > 360) {
00279       cl -= 360;
00280       cr--;
00281     }
00282     if (isnan (lncm)) {
00283       dlon = lnhg - cl;
00284       while (dlon > 180) dlon -= 360;
00285       while (dlon < -180) dlon += 360;
00286     }
00287     if (dlon >= 0) {
00288       if (lat >= 0)
00289         sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, dlon, lat);
00290       else if (lat > -0.05)
00291         sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, dlon, -lat);
00292       else sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fS", cr, cl, dlon, -lat);
00293     } else if (dlon > -0.05) {
00294       if (lat >= 0)
00295         sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, -dlon, lat);
00296       else if (lat > -0.05)
00297         sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fN", cr, cl, -dlon, -lat);
00298       else sprintf (fnam, "%04d:%05.1f_%04.1fW%04.1fS", cr, cl, -dlon, -lat);
00299     } else {
00300       if (lat >= 0)
00301         sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fN", cr, cl, -dlon, lat);
00302       else if (lat > -0.05)
00303         sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fN", cr, cl, -dlon, -lat);
00304       else sprintf (fnam, "%04d:%05.1f_%04.1fE%04.1fS", cr, cl, -dlon, -lat);
00305     }
00306   } else {
00307     if (isnan (dlon)) {
00308                       /*  construct a filename lonHG+lat  */
00309       sprintf (fnam, "%04.1f%+05.1f", lnhg, lat);
00310       return fnam;
00311     }
00312                      /*  construct a filename lonEWlatNS  */
00313     if (dlon >= 0) {
00314       if (lat >= 0) sprintf (fnam, "%04.1fW%04.1fN", dlon, lat);
00315       else if (lat > -0.05) sprintf (fnam, "%04.1fW%04.1fN", dlon, -lat);
00316       else sprintf (fnam, "%04.1fW%04.1fS", dlon, -lat);
00317     } else if (dlon > -0.05) {
00318       if (lat >= 0) sprintf (fnam, "%04.1fW%04.1fN", -dlon, lat);
00319       else if (lat > -0.05) sprintf (fnam, "%04.1fW%04.1fN", -dlon, -lat);
00320       else sprintf (fnam, "%04.1fW%04.1fS", -dlon, -lat);
00321     } else {
00322       if (lat >= 0) sprintf (fnam, "%04.1fE%04.1fN", -dlon, lat);
00323       else if (lat > -0.05) sprintf (fnam, "%04.1fE%04.1fN", -dlon, -lat);
00324       else sprintf (fnam, "%04.1fE%04.1fS", -dlon, -lat);
00325     }
00326   }
00327   return fnam;
00328 }
00329                                  /*  module body  */
00330 int DoIt(void)  {
00331   CmdParams_t *params = &cmdparams;
00332   int status = 0;
00333   int drms_input, drms_output;
00334   int nrec, rec_i;
00335   DRMS_RecordSet_t *recordSet = NULL, *kern_set = NULL;
00336   DRMS_Record_t *irec, *orec;
00337   DRMS_Segment_t *segment, *oseg;
00338   FILE *fpt;
00339   double latc, loncm, lonhg, loncCar;
00340   float fval;
00341   int ival;
00342   int n;
00343   char odir[DRMS_MAXPATHLEN], filename[DRMS_MAXPATHLEN+5];
00344   char buffer[1024], outfilex[DRMS_MAXPATHLEN], outfiley[DRMS_MAXPATHLEN];
00345   char outdir[DRMS_MAXPATHLEN], outdirx[DRMS_MAXPATHLEN],
00346       outdiry[DRMS_MAXPATHLEN], kernfile[DRMS_MAXPATHLEN];
00347   char *carstr;
00348   char rec_query[DRMS_MAXQUERYLEN], source[DRMS_MAXQUERYLEN];
00349   char module_ident[64];
00350 
00351   int twosegs = 0;
00352 
00353   char *in = strdup (params_get_str (params, "in"));
00354   int cr = params_get_int (params, "cr");
00355   float cl = params_get_float (params, "clon");
00356   float lat = params_get_float (params, "lat");
00357   float lon = params_get_float (params, "lon");
00358   char *seg = strdup (params_get_str  (params, "seg"));
00359   char *uxseg = strdup (params_get_str  (params, "uxseg"));
00360   char *uyseg = strdup (params_get_str  (params, "uyseg"));
00361   char *kernel = strdup (params_get_str (params, "kernel"));
00362   char *out = strdup (params_get_str (params, "out"));
00363   char *ave = strdup (params_get_str (params, "ave"));
00364   char *coef = strdup (params_get_str (params, "coef"));
00365   double ob = params_get_double (params, "ob");
00366   double oe = params_get_double (params, "oe");
00367   double rb = params_get_double (params, "rb");
00368   double re = params_get_double (params, "re");
00369   double amu = params_get_double (params, "amu");
00370   int num = params_get_int (params, "num");
00371   int verbose = params_isflagset (params, "v");
00372 
00373   int qave = strcmp (ave, "Not Specified");
00374   int qcoef = strcmp (coef, "Not Specified");
00375   int anycr = (cr <= 0);
00376   int anycl = isnan (cl);
00377   int anylon = isnan (lon);
00378   int anylat = isnan (lat);
00379 
00380   snprintf (module_ident, 64, "%s v %s", module_name, version_id);
00381   if (verbose) printf ("%s:\n", module_ident);
00382 
00383   drms_input = 1;
00384              /*  attempt to open the kernel as a DRMS record set  */
00385   kern_set = drms_open_records (drms_env, kernel, &status);
00386   if (status) {
00387                    /*  attempt to open the kernel as a file name  */
00388     fpt = fopen (kernel, "r");
00389     if (!fpt) {
00390       fprintf (stderr, "Error: kernel specification \"%s\"\n", kernel);
00391       fprintf (stderr,
00392           "       does not specify a readable data record nor file\n");
00393       return 1;
00394     }
00395     fclose (fpt);
00396     strcpy (kernfile, kernel);
00397   } else {
00398                /*  check that the kernel set specifies unique record  */
00399     if (kern_set->n != 1) {
00400       fprintf (stderr, "Error: no data records in set %s\n", kernel);
00401       return 1;
00402     }
00403     irec = kern_set->records[0];
00404     drms_record_directory (irec, kernfile, 1);
00405     segment = drms_segment_lookup (irec, "kernel");
00406     if (!segment) {
00407                    /*  and that it contains the segment "kernel"  */
00408       fprintf (stderr, "Error: kernel record \"%s\"\n", kernel);
00409       fprintf (stderr, "       does not contain a segment \"kernel\"\n");
00410       drms_close_records (kern_set, DRMS_FREE_RECORD);
00411       return 1;
00412     }
00413     sprintf (kernfile, "%s/S%05i/kernel", segment->record->su->sudir,
00414     segment->record->slotnum);
00415     drms_close_records (kern_set, DRMS_FREE_RECORD);
00416   }
00417 
00418   if (anycr && anycl && anylon && anylat) {
00419          /*  no target range specified, assume implicit in input set  */
00420              /*  but it might be implicit in output series spec!  */
00421     recordSet = drms_open_records (drms_env, in, &status);
00422     if (status) {
00423       fpt = fopen (in, "r");
00424       if (!fpt) {
00425     fprintf (stderr, "Error: in specification \"%s\"\n", in);
00426     fprintf (stderr,
00427         "       does not specify a readable data set nor file\n");
00428     return 1;
00429       }
00430       fclose (fpt);
00431       drms_input = 0;
00432       nrec = 1;
00433     } else nrec = recordSet->n;
00434     if (nrec < 1) {
00435       fprintf (stderr, "Error: no data records in set %s\n", in);
00436       return 1;
00437     }
00438   } else {
00439     strncpy (rec_query, in, DRMS_MAXQUERYLEN);
00440     if (!anycr) {
00441       snprintf (source, DRMS_MAXQUERYLEN, "[?CarrRot=%d?]", cr);
00442       strncat (rec_query, source, DRMS_MAXQUERYLEN);
00443     }
00444     if (!anycl) {
00445       float clmin = cl - 0.01;
00446       float clmax = cl + 0.01;
00447                         /*  need to deal with 0/360  */
00448       if (clmin < 0.0) {
00449         clmin += 360.0;
00450     snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g or CMLon<%g?]",
00451         clmin, clmax);
00452       } else if (clmax > 360.0) {
00453         clmax -= 360.0;
00454     snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g or CMLon<%g?]",
00455         clmin, clmax);
00456       } else snprintf (source, DRMS_MAXQUERYLEN, "[?CMLon>%g and CMLon<%g?]",
00457       clmin, clmax);
00458       strncat (rec_query, source, DRMS_MAXQUERYLEN);
00459     }
00460     if (!anylon) {
00461       float lonmin = lon - 0.01;
00462       float lonmax = lon + 0.01;
00463                         /*  need to deal with 0/360  */
00464       if (lonmin < 0.0) {
00465         lonmin += 360.0;
00466     snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g or LonHG<%g?]",
00467         lonmin, lonmax);
00468       } else if (lonmax > 360.0) {
00469         lonmax -= 360.0;
00470     snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g or LonHG<%g?]",
00471         lonmin, lonmax);
00472       } else snprintf (source, DRMS_MAXQUERYLEN, "[?LonHG>%g and LonHG<%g?]",
00473       lonmin, lonmax);
00474       strncat (rec_query, source, DRMS_MAXQUERYLEN);
00475     }
00476     if (!anylat) {
00477       float latmin = lat - 0.01;
00478       float latmax = lat + 0.01;
00479       snprintf (source, DRMS_MAXQUERYLEN, "[?LatHG>%g and LatHG<%g?]",
00480       latmin, latmax);
00481       strncat (rec_query, source, DRMS_MAXQUERYLEN);
00482     }
00483     if (!(recordSet = drms_open_records (drms_env, rec_query, &status))) {
00484       fprintf (stderr, "Error: unable to open input data set %s\n", rec_query);
00485       fprintf (stderr, "       status = %d\n", status);
00486       return 1;
00487     }
00488     nrec = recordSet->n;
00489   }
00490   if (nrec < 1) {
00491     fprintf (stderr, "Error: no records found in input data set\n");
00492     fprintf (stderr, "       %s\n", rec_query);
00493   }
00494   if (drms_input) {
00495         /*  check that range is restricted by dataset specification  */
00496     irec = recordSet->records[0];
00497     if (anycr) {
00498       anycr = unique_key_values (in, "CarrRot", irec) - 1;
00499       ival = drms_getkey_int (irec, "CarrRot", &status);
00500       if (!status) cr = ival;
00501     }
00502     if (anycl) {
00503       anycl = unique_key_values (in, "CMLon", irec) - 1;
00504       fval = drms_getkey_float (irec, "CMLon", &status);
00505       if (!status) cl = fval;
00506     }
00507     if (anylon) {
00508       anylon = unique_key_values (in, "LonHG", irec) - 1;
00509       fval = drms_getkey_float (irec, "LonHG", &status);
00510       if (!status) lon = fval;
00511     }
00512     if (anylat) {
00513       anylat = unique_key_values (in, "LatHG", irec) - 1;
00514       fval = drms_getkey_float (irec, "LatHG", &status);
00515       if (!status) lat = fval;
00516     }
00517     if (verbose) printf ("Processing %d records\n", nrec);
00518   }
00519                    /*  Attempt to create output drms record  */
00520   drms_output = 1;
00521   orec = drms_create_record (drms_env, out, DRMS_PERMANENT, &status);
00522   if (status) {
00523     drms_output = 0;
00524     fprintf (stderr,
00525         "Warning: drms_create_record() returned %d for data series:\n", status);
00526     fprintf (stderr,
00527     "       %s\n       will be interpreted as directory name\n", out);
00528     strncpy (odir, out, DRMS_MAXPATHLEN);
00529   } else if (drms_record_directory (orec, odir, 1)) {
00530     fprintf (stderr,
00531         "Error: drms_record_directory() failure: SUMS may not be available\n");
00532     fprintf (stderr,
00533         "       or series %s may contain no segments\n", out);
00534     return 1;
00535   }
00536   if (drms_output) {
00537     int kstat = 0;
00538     int keyct = sizeof (propagate) / sizeof (char *);
00539 
00540     kstat += check_and_set_key_str (orec, "Module", module_ident);
00541     kstat += check_and_set_key_time (orec, "Created", CURRENT_SYSTEM_TIME);
00542                            /*  module specific keywords  */
00543     kstat += check_and_set_key_float (orec, "amu", amu);
00544     kstat += check_and_set_key_float (orec, "freqmin", ob);
00545     kstat += check_and_set_key_float (orec, "freqmax", oe);
00546     kstat += check_and_set_key_float (orec, "radmin", rb);
00547     kstat += check_and_set_key_float (orec, "radmax", re);
00548     kstat += check_and_set_key_str (orec, "BLD_VERS", jsoc_version);
00549                  /*  set possible prime keys as appropriate  */
00550     if (anycr) {
00551       if (key_is_prime (orec, "CarrRot")) {
00552     fprintf (stderr,
00553         "Error: CarrRot is prime key for output series %s\n", out);
00554     fprintf (stderr,
00555         "       but input data set contains multiple values\n");
00556     drms_close_record (orec, DRMS_FREE_RECORD);
00557                      /*  this pretty much has to be true...  */
00558     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00559     return 1;
00560       }
00561     } else kstat += check_and_set_key_int (orec, "CarrRot", cr);
00562     if (anycl) {
00563       if (key_is_prime (orec, "CMLon")) {
00564     fprintf (stderr,
00565         "Error: CMLon is prime key for output series %s\n", out);
00566     fprintf (stderr,
00567         "       but input data set contains multiple values\n");
00568     drms_close_record (orec, DRMS_FREE_RECORD);
00569     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00570     return 1;
00571       }
00572     } else kstat += check_and_set_key_double (orec, "CMLon", cl);
00573     if (anycr || anycl) {
00574       if (key_is_prime (orec, "CarrTime")) {
00575     fprintf (stderr,
00576         "Error: CarrTime is prime key for output series %s\n", out);
00577     fprintf (stderr,
00578         "       but input data set contains multiple values\n");
00579     drms_close_record (orec, DRMS_FREE_RECORD);
00580     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00581     return 1;
00582       }
00583     } else {
00584       char crcl[32];
00585       sprintf (crcl, "%d:%03.0f", cr, cl);
00586       kstat += check_and_set_key_str (orec, "CarrTime", crcl);
00587     }
00588     if (anylon) {
00589       if (key_is_prime (orec, "LonHG")) {
00590     fprintf (stderr,
00591         "Error: LonHG is prime key for output series %s\n", out);
00592     fprintf (stderr,
00593         "       but input data set contains multiple values\n");
00594     drms_close_record (orec, DRMS_FREE_RECORD);
00595     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00596     return 1;
00597       }
00598     } else kstat += check_and_set_key_double (orec, "LonHG", lon);
00599     if (anylat) {
00600       if (key_is_prime (orec, "LatHG")) {
00601     fprintf (stderr,
00602         "Error: LatHG is prime key for output series %s\n", out);
00603     fprintf (stderr,
00604         "       but input data set contains multiple values\n");
00605     drms_close_record (orec, DRMS_FREE_RECORD);
00606     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00607     return 1;
00608       }
00609     } else kstat += check_and_set_key_double (orec, "LatHG", lat);
00610 
00611     if (kstat) {
00612       fprintf (stderr, "Error writing key value(s) to record in series %s:\n",
00613       out);
00614       fprintf (stderr, "      series may not have appropriate structure;\n");
00615       drms_close_record (orec, DRMS_FREE_RECORD);
00616       if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00617       return 1;
00618     }
00619                /*  search for separate output segments ux and uy  */
00620     if (drms_segment_lookup (orec, uxseg) && drms_segment_lookup (orec, uyseg))
00621       twosegs = 1;
00622     else {
00623                    /*  search for 1st generic segment and use it  */
00624       int segct = drms_record_numsegments (orec);
00625       if (segct < 1) {
00626     fprintf (stderr,
00627         "Error: no data segment in output series %s\n", out);
00628     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00629     drms_close_record (orec, DRMS_FREE_RECORD);
00630     return 1;
00631       }
00632       for (n = 0; n < segct; n++) {
00633     oseg = drms_segment_lookupnum (orec, n);
00634     if (oseg->info->protocol != DRMS_GENERIC) continue;
00635     break;
00636       }
00637       if (n == segct) {
00638     fprintf (stderr,
00639         "Error: no generic data segment in output series %s\n", out);
00640     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00641     drms_close_record (orec, DRMS_FREE_RECORD);
00642     return 1;
00643       }
00644       fprintf (stderr, "         writing to segment %s\n", oseg->info->name);
00645     }
00646   }
00647   if (twosegs) {
00648     sprintf (outdirx, "%s/%s", odir, uxseg);
00649     sprintf (outdiry, "%s/%s", odir, uyseg);
00650     mkdir (outdirx, 01755);
00651     mkdir (outdiry, 01755);
00652   } else {
00653     if (drms_output) {
00654       sprintf (outdir, "%s/%s", odir, oseg->info->name);
00655       mkdir (outdir, 01755);
00656     } else strcpy (outdir, out);
00657   }
00658                             /*  main processing loop  */
00659   for (rec_i=0; rec_i<nrec; rec_i++) {
00660     if (verbose) printf ("  Processing record %i\n", rec_i);
00661     if (drms_input) {
00662       irec = recordSet->records[rec_i];
00663               /*  should use call to drms_record_directory() instead  */
00664       if (irec->sunum != -1LL && irec->su == NULL) {
00665     irec->su = drms_getunit (irec->env, irec->seriesinfo->seriesname,
00666         irec->sunum, 1, &status);
00667       }
00668       segment = drms_segment_lookup(irec, seg);
00669       sprintf (filename, "%s/S%05i/%s", segment->record->su->sudir,
00670       segment->record->slotnum, seg);
00671       drms_sprint_rec_query (source, irec);
00672       cr = drms_getkey_int (irec, "CarrRot", &status);
00673       latc = drms_getkey_double (irec, "LatHG", &status);
00674       lonhg = drms_getkey_double (irec, "LonHG", &status);
00675       carstr = drms_getkey_string (irec, "CarrTime", &status);
00676       if (status) loncCar = drms_getkey_double (irec, "CMLon", &status);
00677       else sscanf (carstr, "%i:%lf", &cr, &loncCar);
00678       loncm = drms_getkey_double (irec, "LonCM", &status);
00679     } else {
00680       strcpy (filename, in);
00681       strcpy (source, in);
00682       latc = loncm = lonhg = 0.0;
00683     }
00684     if (twosegs) {
00685       sprintf (outfilex, "%s/%s/%s.Ux", odir, uxseg, 
00686       printcrcl_loc (cr, loncCar, lonhg, loncm, latc));
00687       sprintf (outfiley, "%s/%s/%s.Uy", odir, uyseg,
00688       printcrcl_loc (cr, loncCar, lonhg, loncm, latc));
00689     } else {
00690       sprintf (outfilex, "%s/%s.Ux", outdir, 
00691       printcrcl_loc (cr, loncCar, lonhg, loncm, latc));
00692       sprintf (outfiley, "%s/%s.Uy", outdir,
00693       printcrcl_loc (cr, loncCar, lonhg, loncm, latc));
00694     }
00695     status = process_record (filename, drms_output, outfilex, outfiley,
00696     kernfile, kernel, ave, coef, qave, qcoef, verbose, ob, oe, num, rb, re, amu,
00697     source, module_ident);
00698     if (status) {
00699       fprintf (stderr, "Error processing record %d (%s); aborted\n", rec_i,
00700       printcrcl_loc (cr, loncCar, lonhg, loncm, latc));
00701       if (drms_output) {
00702     if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00703     drms_close_record (orec, DRMS_FREE_RECORD);
00704     return 1;
00705       }
00706     }
00707   }
00708 
00709   if (drms_input) drms_close_records (recordSet, DRMS_FREE_RECORD);
00710   if (drms_output) drms_close_record (orec, DRMS_INSERT_RECORD);
00711 
00712   return 0;
00713 }
00714 
00715 /*
00716  *  Revision History
00717  *
00718  *  June 2009 - v 0.1 module created, no real drms functionality
00719  *  09.10.2009 - v 0.2 drms functionality added - input and output from drms
00720  *  records now  possible, but input/output to rooted filenames also
00721  *  possible
00722  *  0.3 2009.10.22  fixed format of name of output files
00723  *  0.4 2010.02.11  unified processing calls for named file and drms input
00724  *  branches
00725  *  0.5 2010.02.12  Combined branches into one, added some very obvious
00726  *  error handling (no given input argument, input cannot be opened)
00727  *          removed rdutil.h inclusion
00728  *  v 0.5 frozen 2010.03.08
00729  *    0.6   Added logic for specifying input grouping of data set records
00730  *      by Carr Rot, CM Longitude, Longitude, and/or Latitude
00731  *      Added keyword propagation list (but not propagation!)
00732  *      Output to multiple files named by lat and lon in segment
00733  *      directory
00734  *      Fixed processing of longitude queries near 0/360
00735  *  2010.03.10  Fixed setting of keys for run parameters; put in trap
00736  *      for empty input data sets and widened clon match range
00737  *  v 0.6 frozen 2010.04.23
00738  *    0.7       Output data files are now written into per segment
00739  *      subdirectories; added option for naming individual segments
00740  *      (if there are two)
00741  *  v 0.7 frozen 2010.08.19
00742  *    0.8   Created 2011.04.22; Fixed typos in ola_xy.c (not yet used)
00743  *      Fixed bugs in printing of filenames (particularly the
00744  *      reversal of east and west), and dropped requirement for CMLon;
00745  *      heliographic option for file naming
00746  *      Removed default for kernel
00747  *    0.9   Created 2011.12.04; Fixed typos in ola_xy.c (not yet used)
00748  *  v 0.8 frozen 2012.02.01 (but included in JSOC release 6.1 2011.12.05 in
00749  *      identical form except for comments)
00750  *    0.9   added (preferred) option for specifying kernel as DRMS data
00751  *      record rather than filename
00752  *  v 0.9 frozen 2012.04.24
00753  *    0.91  similar to 0.9, but includes fixes in FORTRAN code ola_xy_v12
00754  *      and provision for return status
00755  *  v 0.91 frozen 2012.04.24
00756  */