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

00001 /*
00002  *  datavg.c                        ~rick/src/drms
00003  *
00004  *  Construct mean and variance matrices of selected data segments from an
00005  *    input record set (same format)
00006  *  This is a DRMS version of SOI module doppavg generalized to a
00007  *    larger variety of data series and types
00008  *
00009  *  Responsible:  Rick Bogart               RBogart@spd.aas.org
00010  *
00011  *  Usage:
00012  *    datavg in= out= 
00013  *
00014  *  Parameters: (type   default         description)
00015  *  in  str -       Input dataset or data series;
00016  *          if only a data series name, additional selection
00017  *          parameters will be required
00018  *      out str -       Output data series
00019  *      count   str valid       Output series segment containing count
00020  *      mean    str mean        Output series segment containing mean
00021  *      power   str power       Output series segment containing variance
00022  *      log str Log     Output series segment containing run log
00023  *      tmid    str unspecified midpoint of averaging interval (CR:CL)
00024  *  length  float   unspecified length of averaging interval (in deg og
00025  *          Carrington rotation)
00026  *  reject  string  Not Specified   If specified, name of a file with a
00027  *              list of input images to be rejected regardless
00028  *              of quality mask values
00029  *  qmask   int 0x80000000  Quality mask for data acceptability;
00030  *              records rejected if (qmask & qkey:value) != 0
00031  *  cvok    str any     64-bit mask of acceptable values for
00032  *              calibration veraion of input; default -> all
00033  *              values acceptable
00034  *  cvno    str none        64-bit mask of acceptable values for
00035  *              calibration veraion of input; default -> no
00036  *              values unacceptable
00037  *  copy    str +       list of keys to be propagated as-is
00038  *  average str +       list of keys to be averaged
00039  *  pkey    str T_OBS       Name of keyword to be used as index
00040  *          over which input records are selected for averaging
00041  *  qual_key str    Quality     Key name of uint type keyword describing
00042  *          data quality
00043  *  cvkey   str CalVer64    Key name of 64-bit mask type keyword
00044  *          describing calibration version used
00045  *  roll_key str    CROTA2      Name of keyword describing roll angle
00046  *  pa  float   180.0       Centre of acceptable roll angles
00047  *  dpa float   1.0     Maximum deviation of acceptable roll
00048  *          angles from nominal centre
00049  *  mscale  float   (SEG default)   Output BSCALE factor for mean
00050  *  mzero   float   (SEG default)   Output BZERO offset for mean
00051  *  pscale  float   (SEG default)   Output BSCALE factor for power
00052  *  pzero   float   (SEG default)   Output BZERO offset for power
00053  *  setkey  str Not Specified   Name of special keyword to be set
00054  *  setval  double  Not Specified   Value for special keyword to be set, or
00055  *          name of key whose value is to be used
00056  *
00057  *  Flags
00058  *  -n  no ingestion of output; diagnostic only
00059  *  -o  correct individual images for orbital velocity
00060  *  -v  run verbose
00061  *
00062  *  Bugs:
00063  *    Although drms_segment_write() errors are flagged in stderr, they do
00064  *  not force an abort; consequently it is possible to have records
00065  *  created with incomplete SUMS data, or no SUMS data at all.
00066  *    Propagated keys are copied from the first input record, without checking
00067  *  for uniqueness
00068  *    If the input dataset is specified in the in string, the keywords for
00069  *  CarRot and CMLon are set to garbage, and there is no checking that
00070  *  the length is correct
00071  *    The MidTime value is geoecentric for midpoint Carrington longitude
00072  *    The correction for observer velocity only takes account of the radial
00073  *  component, and is only appropriate for Doppler data
00074  *    For log-base data, the LOG_BASE keyword is propagated to the output
00075  *  record, but it really only applies to the mean and variance; it should
00076  *  be a per-segment keyword
00077  *    The check on nominal roll angle and the averagang of CROTA* values
00078  *  assumes that the units are degrees
00079  *    The reported extremal values of mean and power are determined before
00080  *  possible output scaling; they may be outside the representable range
00081  *    Only a single value can be specified for cvok or cvno (which are treated
00082  *  as string arguments for local parsing, although in principle they
00083  *  should be arrays of ints)
00084  *    The protocol of the log segment, if present in the output series, is not
00085  *  verified to be generic
00086  *    The values of the primekey in the input set, as specified by pkey, are
00087  *  automatically averaged into a corresponding key (and D_key) in the
00088  *  output data set. If the default value of T_OBS is overridden, must
00089  *  provide appropriate keys in output series if desired; alternatively,
00090  *  T_OBS should be added to average list.
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 
00100 #define DO_NOTHING  (0)
00101 #define DO_ABSVAL   (1)
00102 #define DO_SQUARE   (2)
00103 
00104 char *module_name = "data average";
00105 char *version_id = "1.5";
00106 
00107 ModuleArgs_t module_args[] = {
00108   {ARG_STRING,  "in", "", "input data series or dataset"}, 
00109   {ARG_STRING,  "out", "", "output data series"}, 
00110   {ARG_STRING,  "tmid", "Not Specified",
00111       "midpoint of averaging interval (in Carrington longitude)"}, 
00112   {ARG_FLOAT,   "length", "Not Specified",
00113       "length of averaging interval (in degrees of Carrington rotation)"}, 
00114   {ARG_FLOAT,   "pa", "180.0", "centre of acceptable roll angles"}, 
00115   {ARG_FLOAT,   "dpa", "1.0", "maximum deviation of acceptable roll angles"}, 
00116   {ARG_INT, "qmask", "0x80000000", "quality bit mask for image rejection"},
00117   {ARG_STRING,  "reject", "Not Specified", "file containing rejection list"}, 
00118   {ARG_STRING,  "cvkey", "CalVer64", "keyname for Calibration Version key"}, 
00119   {ARG_STRING,  "cvok", "any", "Acceptable value of cvkey"},
00120   {ARG_STRING,  "cvno", "none", "Unacceptable value of cvkey"},
00121   {ARG_STRING,  "copy",  "+", "comma separated list of keys to propagate"},
00122   {ARG_STRING,  "average",  "+", "comma separated list of keys to average"},
00123   {ARG_STRING,  "pkey", "T_OBS",
00124       "keyname for index over which records are selected for averaging"}, 
00125   {ARG_STRING,  "qual_key", "Quality", "keyname for 32-bit image quality field"}, 
00126   {ARG_STRING,  "roll_key", "CROTA2", "keyname for input roll-angle field"}, 
00127   {ARG_STRING,  "count", "valid", "output data series segment containing count"}, 
00128   {ARG_STRING,  "mean", "mean", "output data series segment containing mean"}, 
00129   {ARG_STRING,  "power", "power", "output data series segment containing variance"}, 
00130   {ARG_STRING,  "log", "Log", "output data series segment containing run log"}, 
00131   {ARG_STRING,  "setkey", "Not Specified", "name of special extra key to be set"}, 
00132   {ARG_DOUBLE,  "setval", "Not Specified",
00133     "value of special extra key to be set; if invalid, name of key whose value is to be used"}, 
00134   {ARG_FLOAT,   "mscale", "Segment Default", "output BSCALE factor for mean"},
00135   {ARG_FLOAT,   "mzero", "Segment Default", "output BZERO offset for mean"},
00136   {ARG_FLOAT,   "pscale", "Segment Default", "output BSCALE factor for power"},
00137   {ARG_FLOAT,   "pzero", "Segment Default", "output BZERO offset for power"},
00138   {ARG_FLAG,    "n",    "", "no output records produced; diagnostics only"}, 
00139   {ARG_FLAG,    "o",    "", "remove effects of observer velocity"}, 
00140   {ARG_FLAG,    "v",    "", "verbose mode"}, 
00141   {ARG_END}
00142 };
00143 
00144     /*  list of keywords to propagate by default (if possible and unique)
00145                                 from input to output  */
00146 char *propagate[] = {"TELESCOP", "INSTRUME", "WCSNAME", "WCSAXES"};
00147          /*  list of keywords to average by default (if possible)
00148                                 from input to output  */
00149 char *average[] = {"T_OBS", "OBS_VR", "OBS_VW", "OBS_VN"};
00150               /*  the following belong in external utility files  */
00151 static long long params_get_mask (CmdParams_t *params, char *arg,
00152     long long defval) {
00153 /*
00154  *  This function parses the string associated with the command parameters
00155  *    argument "arg" as the hexadecimal representation of an unsigned 64-bit
00156  *    integer, which it returns. The string may consist of up to 16 hexadecimal
00157  *    characters. (They can optionally be preceded by '0x', but the string is
00158  *    treated as a hexadecimal representation regardless.)  If there are any
00159  *    extra or illegal characters in the string, the value "defval" is returned.
00160  */
00161   long long retval;
00162   const char *str = params_get_str (params, arg);
00163   char *ext;
00164 
00165   retval  = strtoull (str, &ext, 16);
00166   if (strlen (ext)) retval = defval;
00167   return retval;
00168 }
00169 
00170 static int fgetline (FILE *in, char *line, int max) {
00171   if (fgets (line, max, in) == NULL) return 0;
00172   else return (strlen (line));
00173 }
00174 
00175 static int read_reject_list (FILE *file, int **list) {
00176   int ds, sn, rec, last_rec;
00177   int allocd = 1024, ct = 0, gap = 0;
00178   char line[1024], t_str[64], estr[16];
00179 
00180   *list = (int *)malloc (allocd * sizeof (int));
00181   while (fgetline (file, line, 1024)) {
00182     if (strlen (line) == 1) continue;
00183     if (line[0] == '#') continue;
00184     if (sscanf (line, "%d %d %d %s", &ds, &sn, &rec, t_str) != 4) {
00185       if (sscanf (line, "%d %s", &rec, t_str) != 2) {
00186     sscanf (line, "%s", estr);
00187     if (strcmp (estr, "...")) continue;
00188     gap = 1;
00189     last_rec = rec;
00190     continue;
00191       }
00192     }
00193     if (gap) {
00194       while (rec > last_rec) {
00195     last_rec++;
00196     (*list)[ct++] = last_rec;
00197     if (ct >= allocd) {
00198       allocd += 1024;
00199       *list = (int *)realloc (*list, allocd * sizeof (int));
00200     }
00201       }
00202       gap = 0;
00203       continue;
00204     }
00205     (*list)[ct++] = rec;
00206     if (ct >= allocd) {
00207       allocd += 1024;
00208       *list = (int *)realloc (*list, allocd * sizeof (int));
00209     }
00210   }
00211   return ct;
00212 }
00213 
00214 static int key_params_from_dspec (const char *dspec) {
00215 /*
00216  *  Establish whether input data records are determined from dataset specifier
00217  *  assume that if a bracket is found in the dataset specifier it is a
00218  *  set of well-specified records containing a properly ordered input dataset;
00219  *  otherwise, a dataseries to be queried
00220  */
00221   int n, nt = strlen (dspec);
00222 
00223   for (n = 0; n < nt; n++) if (dspec[n] == '[') return 1;
00224   return 0;
00225 }
00226 
00227 int set_stats_keys (DRMS_Record_t *rec, DRMS_Array_t *vcts, DRMS_Array_t *mean,
00228     DRMS_Array_t *powr, long long ntot, int mscaled, int pscaled,
00229     double mrepmin, double mrepmax, double prepmin, double prepmax, int verbose) {
00230   double *vavg = (double *)mean->data;
00231   double *vvar = (powr) ? (double *)powr->data : NULL;
00232   int *vval = (int *)vcts->data;
00233   double mnvmin, mnvmax, pwrmin, pwrmax, norm, norm2, scale, sig, var;
00234   double sumv1, sump1, sum2, sum3, sum4;
00235   int sumv0, sump0;
00236   int valmin = ntot, valmax = 0;
00237   int n;
00238   int vminloc, vmaxloc, pminloc, pmaxloc;
00239   int morhict, morloct, morct, porhict, porloct, porct;
00240   int kstat = 0;
00241 
00242   mnvmin = pwrmin = DBL_MAX;
00243   mnvmax = pwrmax = -DBL_MAX;
00244   sumv1 = sump1 = 0.0;
00245   sumv0 = sump0 = 0;
00246   if (mscaled || pscaled) {
00247     morhict = morloct = 0;
00248     porhict = porloct = 0;
00249   }
00250 
00251   for (n = 0; n < ntot; n++) {
00252     if (vval[n] < valmin) valmin = vval[n];
00253     if (vval[n] > valmax) valmax = vval[n];
00254     if (vavg[n] < mnvmin) {
00255       mnvmin = vavg[n];
00256       vminloc = n;
00257     }
00258     if (vavg[n] > mnvmax) {
00259       mnvmax = vavg[n];
00260       vmaxloc = n;
00261     }
00262     if (isfinite (vavg[n])) {
00263       sumv0++;
00264       sumv1 += vavg[n];
00265     }
00266     if (mscaled) {
00267       if (vavg[n] < mrepmin)  morloct++;
00268       if (vavg[n] > mrepmax)  morhict++;
00269     }
00270     if (vvar) {
00271       if (vvar[n] < pwrmin) {
00272     pwrmin = vvar[n];
00273     pminloc = n;
00274       }
00275       if (vvar[n] > pwrmax) {
00276     pwrmax = vvar[n];
00277     pmaxloc = n;
00278       }
00279       if (isfinite (vvar[n])) {
00280     sump0++;
00281     sump1 += vvar[n];
00282       }
00283       if (pscaled) {
00284     if (vvar[n] < prepmin)  porloct++;
00285     if (vvar[n] > prepmax)  porhict++;
00286       }
00287     }
00288   }
00289   
00290   kstat += check_and_set_key_int (rec, "CountMIN", valmin);
00291   kstat += check_and_set_key_int (rec, "CountMAX", valmax);
00292   kstat += check_and_set_key_double (rec, "MeanMIN", mnvmin);
00293   kstat += check_and_set_key_double (rec, "MeanMAX", mnvmax);
00294   if (vvar) {
00295     kstat += check_and_set_key_double (rec, "PowrMIN", pwrmin);
00296     kstat += check_and_set_key_double (rec, "PowrMAX", pwrmax);
00297   }
00298 
00299   if (sumv0) {
00300     scale = 1.0 / sumv0;
00301     sumv1 *= scale;
00302     sum2 = sum3 = sum4 = 0.0;
00303     kstat += check_and_set_key_double (rec, "MeanAVG", sumv1);
00304     for (n = 0; n < ntot; n++) {
00305       if (isfinite (vavg[n])) {
00306         norm = vavg[n] - sumv1;
00307     norm2 = norm * norm;
00308     sum2 += norm2;
00309     sum3 += norm * norm2;
00310     sum4 += norm2 * norm2;
00311       }
00312     }
00313     kstat += check_and_set_key_double (rec, "MeanRMS", sqrt (sum2 / sumv0));
00314     if (sumv0 > 1) {
00315       var = sum2 / (sumv0 - 1);
00316       sig = sqrt (var);
00317       kstat += check_and_set_key_double (rec, "MeanSKEW",
00318       scale * sum3 / (var * sig));
00319       kstat += check_and_set_key_double (rec, "MeanKURT",
00320       scale * sum4 / (var * var) - 3.0);
00321     }
00322   }
00323 
00324   if (sump0) {
00325     scale = 1.0 / sump0;
00326     sump1 *= scale;
00327     sum2 = sum3 = sum4 = 0.0;
00328     kstat += check_and_set_key_double (rec, "PowrAVG", sump1);
00329     for (n = 0; n < ntot; n++) {
00330       if (isfinite (vvar[n])) {
00331         norm = vvar[n] - sump1;
00332     norm2 = norm * norm;
00333     sum2 += norm2;
00334     sum3 += norm * norm2;
00335     sum4 += norm2 * norm2;
00336       }
00337     }
00338     kstat += check_and_set_key_double (rec, "PowrRMS", sqrt (sum2 / sump0));
00339     if (sump0 > 1) {
00340       var = sum2 / (sump0 - 1);
00341       sig = sqrt (var);
00342       kstat += check_and_set_key_double (rec, "PowrSKEW",
00343       scale * sum3 / (var * sig));
00344       kstat += check_and_set_key_double (rec, "PowrKURT",
00345       scale * sum4 / (var * var) - 3.0);
00346     }
00347   }
00348 
00349   if (verbose) {
00350     int rank = mean->naxis;
00351     int rmnd;
00352 
00353     printf ("Extrema: minimum value = %.4e @ %d [", mnvmin, vminloc);
00354     rmnd = vminloc;
00355     for (n = 0; n < rank; n++) {
00356       if (n) printf (", ");
00357       printf ("%d", rmnd %  mean->axis[n]);
00358       rmnd /= mean->axis[n];
00359     }
00360     printf ("]\n");
00361 
00362     printf ("         maximum value = %.4e @ %d [", mnvmax, vmaxloc);
00363     rmnd = vmaxloc;
00364     for (n = 0; n < rank; n++) {
00365       if (n) printf (", ");
00366       printf ("%d", rmnd %  mean->axis[n]);
00367       rmnd /= mean->axis[n];
00368     }
00369     printf ("]\n");
00370 
00371     if (vvar) {
00372       printf ("         minimum power = %.4e @ %d [", pwrmin, pminloc);
00373       rmnd = pminloc;
00374       for (n = 0; n < rank; n++) {
00375     if (n) printf (", ");
00376     printf ("%d", rmnd %  mean->axis[n]);
00377     rmnd /= mean->axis[n];
00378       }
00379       printf ("]\n");
00380 
00381       printf ("         maximum power = %.4e @ %d [", pwrmax, pmaxloc);
00382       rmnd = pmaxloc;
00383       for (n = 0; n < rank; n++) {
00384     if (n) printf (", ");
00385     printf ("%d", rmnd %  mean->axis[n]);
00386     rmnd /= mean->axis[n];
00387       }
00388       printf ("]\n");
00389     }
00390   }
00391 
00392   if (mscaled || pscaled) {
00393     morct = morloct + morhict;
00394     porct = porloct + porhict;
00395     if (morct) {
00396       fprintf (stderr,
00397       "Warning: %d average values out of representable range\n", morct);
00398       fprintf (stderr,
00399       "         for BScale = %g and BZero = %g\n", mean->bscale,
00400       mean->bzero);
00401       fprintf (stderr,
00402       "         %d values < %g, %d values > %g\n", morloct, mrepmin,
00403       morhict, mrepmax);
00404     }
00405     if (porct) {
00406       fprintf (stderr,
00407       "Warning: %d variance values out of representable range\n", porct);
00408       fprintf (stderr,
00409       "         for PScale = %g and PZero = %g\n", powr->bscale,
00410       powr->bzero);
00411       fprintf (stderr,
00412       "         %d values < %g, %d values > %g\n", porloct, prepmin,
00413       porhict, prepmax);
00414     }
00415   }
00416 
00417   return kstat;
00418 }
00419 
00420 void report_pkey_value (FILE *out, DRMS_Record_t *rec, char *key, int recnum) {
00421   DRMS_Keyword_t *keyword;
00422   char buf[128];
00423 
00424   fprintf (out, "record ");
00425   keyword = drms_keyword_lookup (rec, key, 1);
00426   if (!keyword) {
00427     fprintf (out, "%d (unknown value) [#%lld]: ", recnum, rec->recnum);
00428     return;
00429   }
00430   drms_keyword_snprintfval (keyword, buf, sizeof (buf));
00431   fprintf (out, "%s (%d) [#%lld]: ", buf, recnum, rec->recnum);
00432   fflush (out);
00433 }
00434 
00435 int DoIt (void) {
00436   CmdParams_t *params = &cmdparams;
00437   DRMS_RecordSet_t *ids = NULL;
00438   DRMS_Record_t *irec, *orec = NULL;
00439   DRMS_Segment_t *iseg, *vseg = NULL, *mseg = NULL, *pseg = NULL;
00440   DRMS_Segment_t *logseg = NULL;
00441   DRMS_Array_t *data_array, *vcts, *mean, *powr;
00442   DRMS_Keyword_t *keywd;
00443   TIME tmid, tobs_rec, /* tobs, */ tfirst, tlast;
00444   FILE *runlog;
00445   double *v, *vavg, *vvar;
00446   double *crpix, *crval, *cdelt, *crota, *avgval;
00447   double *crpixv, *crvalv, *cdeltv, *crotav, *avgvalv;
00448   double crpix_rec, crval_rec, cdelt_rec, crota_rec, avg_rec;
00449   double /* tobsv, */ vobs;
00450   double log_base;
00451   double mrepmin, mrepmax, prepmin, prepmax;
00452   double crlnobs, crlnval, crlnvalv;
00453   float clstrt, clmid, clstop;
00454   float pa_rec, dpa;
00455   long long *cvgolist, *cvnolist, *cvlist;
00456   long long ntot, calver;
00457   unsigned int quality;
00458   int *inaxis, *vval, *reject_list, *cvfound;
00459   int rec, recct, segct, segnum, maxct, imgct;
00460   int crstrt, crmid, crstop, carrot;
00461   int checkseg, mscaled, pscaled;
00462   int kstat, status;
00463   int log_status;
00464   int writelog;
00465   int n, naxis, wcsaxes;
00466   int propct, meanct, add_defaults, cvct, cvgoct, cvnoct, cvused, cvmaxct;
00467   char **copykeylist, **meankeylist;
00468   char *source, *keystr;
00469   char recset_query[DRMS_MAXQUERYLEN], keyname[DRMS_MAXKEYNAMELEN];
00470   char calverinfo[DRMS_MAXQUERYLEN];
00471   char logfilename[DRMS_MAXPATHLEN];
00472   char module_ident[64], tbuf[64], vbuf[128], strbuf[128];
00473 
00474   double fp_nan = 0.0 / 0.0;
00475   int badqual = 0, blacklist = 0, badpa = 0, badcv = 0, rejects = 0;
00476   int needpowr = 1;
00477   int propkeyct = sizeof (propagate) / sizeof (char *);
00478   int meankeyct = sizeof (average) / sizeof (char *);
00479 
00480   char *inset = strdup (params_get_str (params, "in"));
00481   char *out_series = strdup (params_get_str (params, "out"));
00482   char *vsegname = strdup (params_get_str (params, "count"));
00483   char *msegname = strdup (params_get_str (params, "mean"));
00484   char *psegname = strdup (params_get_str (params, "power"));
00485   char *logsegname = strdup (params_get_str (params, "log"));
00486   char *tmid_str = strdup (params_get_str (params, "tmid"));
00487   float intrvl = params_get_float (params, "length");
00488   float pa_nom = params_get_float (params, "pa");
00489   float dpa_max = params_get_float (params, "dpa");
00490   unsigned int qmask = cmdparams_get_int64 (params, "qmask", &status);
00491   long long cvaccept = params_get_mask (params, "cvok", -1);
00492   long long cvreject = params_get_mask (params, "cvno", 0);
00493   char *calverkey = strdup (params_get_str (params, "cvkey"));
00494   char *rejectfile = strdup (params_get_str (params, "reject"));
00495   char *propagate_req = strdup (params_get_str (params, "copy"));
00496   char *average_req = strdup (params_get_str (params, "average"));
00497   char *primekey = strdup (params_get_str (params, "pkey"));
00498   char *qual_key = strdup (params_get_str (params, "qual_key"));
00499   char *roll_key = strdup (params_get_str (params, "roll_key"));
00500   char *spec_key = strdup (params_get_str (params, "setkey"));
00501   double spec_val = params_get_double (params, "setval");
00502   char *specuse_key = strdup (params_get_str (params, "setval"));
00503   double mscale = params_get_double (params, "mscale");
00504   double mzero = params_get_double (params, "mzero");
00505   double pscale = params_get_double (params, "pscale");
00506   double pzero = params_get_double (params, "pzero");
00507   int no_save = params_isflagset (params, "n");
00508   int remove_obsvel = params_isflagset (params, "o");
00509   int verbose = params_isflagset (params, "v");
00510 
00511   int dispose = (no_save) ? DRMS_FREE_RECORD : DRMS_INSERT_RECORD;
00512   int check_pa = (dpa_max < 180.0);
00513   int set_extra_key = (strcmp (spec_key, "Not Specified")) ? 1 : 0;
00514   int use_other_key = isnan (spec_val);
00515   int filt_on_calver = (cvreject || ~cvaccept);
00516 
00517   snprintf (module_ident, 64, "%s v %s", module_name, version_id);
00518   if (verbose) printf ("%s:\n", module_ident);
00519              /*  initialize CalVer lists (single element for now  */
00520   cvused = 0;
00521   cvmaxct = 64;
00522   cvfound = (int *)malloc (cvmaxct * sizeof (int));
00523   cvlist = (long long *)malloc (cvmaxct * sizeof (long long));
00524   if (filt_on_calver) {
00525     cvgoct = 1;
00526     cvgolist = (long long *)malloc (cvgoct * sizeof (long long));
00527     cvgolist[0] = cvaccept;
00528     cvnoct = 1;
00529     cvnolist = (long long *)malloc (cvnoct * sizeof (long long));
00530     cvnolist[0] = cvreject;
00531     if (~cvaccept) cvgoct = 0;
00532     if (!cvreject) cvnoct = 0;
00533   }
00534   if (!isfinite (pa_nom)) pa_nom = 180.0;
00535                           /*  additional initializations  */
00536   prepmin = mrepmin = 1.0 / 0.0;
00537   prepmax = mrepmax = - mrepmin;
00538                             /*  create output record  */
00539   if (no_save) {
00540     orec = drms_create_record (drms_env, out_series, DRMS_TRANSIENT, &status);
00541     if (status) {
00542       fprintf (stderr, "Warning: drms_create_record returned %d for data series:\n",
00543       status);
00544       fprintf (stderr, "  %s\nKeyword setting will not be checked\n", out_series);
00545     }
00546   } else {
00547     orec = drms_create_record (drms_env, out_series, DRMS_PERMANENT, &status);
00548     if (status) {
00549       fprintf (stderr, "Error: drms_create_record returned %d for data series:\n",
00550       status);
00551       fprintf (stderr, "  %s\n", out_series);
00552       return 1;
00553     }
00554   }
00555   if ((segct = orec->segments.num_total) < 1) {
00556     fprintf (stderr, "Error: no data segments in output data series:\n");
00557     fprintf (stderr, "  %s\n", out_series);
00558     return 1;
00559   }
00560   vseg = drms_segment_lookup (orec, vsegname);
00561   mseg = drms_segment_lookup (orec, msegname);
00562   pseg = drms_segment_lookup (orec, psegname);
00563   if (!vseg && !mseg && !pseg) {
00564     fprintf (stderr,
00565     "Error: output series %s does not contain any of the required segments:\n",
00566     out_series);
00567     fprintf (stderr, "       \"%s\", \"%s\", \"%s\"\n", vsegname, msegname,
00568     psegname);
00569     drms_close_record (orec, DRMS_FREE_RECORD);
00570     return 1;
00571   }
00572   if (vseg) {
00573     switch (vseg->info->type) {
00574       case (DRMS_TYPE_CHAR):
00575     maxct = SCHAR_MAX;
00576     break;
00577       case (DRMS_TYPE_SHORT):
00578     maxct = SHRT_MAX;
00579     break;
00580       default:
00581     maxct = INT_MAX;
00582     }
00583     maxct += vseg->bzero;
00584     maxct /= vseg->bscale;
00585   } else {
00586     fprintf (stderr,
00587     "Warning: output series %s does not contain the segment %s\n",
00588     out_series, vsegname);
00589     maxct = INT_MAX;
00590   }
00591   if (mseg) {
00592     mscaled = (mseg->info->type == DRMS_TYPE_CHAR) ||
00593               (mseg->info->type == DRMS_TYPE_SHORT) ||
00594               (mseg->info->type == DRMS_TYPE_INT) ||
00595               (mseg->info->type == DRMS_TYPE_LONGLONG);
00596   } else {
00597     fprintf (stderr,
00598     "Warning: output series %s does not contain the segment %s\n",
00599     out_series, msegname);
00600     mscaled = 0;
00601   }
00602   if (pseg) {
00603     pscaled = (pseg->info->type == DRMS_TYPE_CHAR) ||
00604               (pseg->info->type == DRMS_TYPE_SHORT) ||
00605               (pseg->info->type == DRMS_TYPE_INT) ||
00606               (pseg->info->type == DRMS_TYPE_LONGLONG);
00607   } else {
00608     fprintf (stderr,
00609     "Warning: output series %s does not contain the segment %s\n",
00610     out_series, psegname);
00611     needpowr = pscaled = 0;
00612   }
00613   writelog = 0;
00614   logseg = drms_segment_lookup (orec, logsegname);
00615   if (logseg && dispose == DRMS_INSERT_RECORD) {
00616     drms_segment_filename (logseg, logfilename);
00617     runlog = fopen (logfilename, "w");
00618     if (runlog) writelog = 1;
00619   }
00620 
00621                          /*  process input specification  */
00622   if (key_params_from_dspec (inset)) {
00623                     /*  input specified as specific data set  */
00624     if (!(ids = drms_open_records (drms_env, inset, &status))) {
00625       fprintf (stderr, "Error: (%s) unable to open input data set %s\n",
00626         module_ident, inset);
00627       fprintf (stderr, "       status = %d\n", status);
00628       return 1;
00629     }
00630     if ((recct = ids->n) < 2) {
00631       fprintf (stderr, "Error: (%s) <2 records in selected input set\n",
00632       module_ident);
00633       fprintf (stderr, "       %s\n", inset);
00634       drms_close_records (ids, DRMS_FREE_RECORD);
00635       return 1;
00636     }
00637     source = strdup (inset);
00638     for (n = 0; n < strlen (source); n++) if (source[n] == '[') {
00639       source[n] = '\0';
00640       break;
00641     }
00642     tmid = clmid = fp_nan;
00643     crmid = -1;
00644   } else {
00645                 /*  only the input data series is named,
00646                     get record specifications from arguments  */
00647     source = strdup (inset);
00648     if (sscanf (tmid_str, "%d:%f", &crmid, &clmid) == 2) {
00649                          /*  tmid specified as CR:CL  */
00650       if (isnan (clmid)) {
00651     fprintf (stderr,
00652         "Error: target Carrington longitude = %f\n", clmid);
00653     return 1;
00654       }
00655       while (clmid <= 0) {
00656         clmid += 360;
00657     crmid++;
00658       }
00659       while (clmid > 360) {
00660         clmid -= 360;
00661     crmid--;
00662       }
00663       clstrt = clmid + 0.5 * intrvl;
00664       clstop = clmid - 0.5 * intrvl;
00665       crstrt = crstop = crmid;
00666       if (intrvl < 360.0 && clstrt <= 360.0 && clstop >= 0.0) {
00667     sprintf (recset_query,
00668         "%s[?Car_Rot=%d and CRLN_OBS<=%.2f and CRLN_OBS>=%.2f?]", source,
00669         crmid, clstrt, clstop);
00670       } else {
00671     while (clstrt > 360.0) {
00672       clstrt -= 360.0;
00673       crstrt--;
00674     }
00675     while (clstop < 0.0) {
00676       clstop += 360.0;
00677       crstop++;
00678     }
00679     if ((crstop - crstrt) > 1)
00680       sprintf (recset_query,
00681           "%s[?(Car_Rot=%d and CRLN_OBS<=%.2f) or (Car_Rot=%d and CRLN_OBS>=%.2f) or (CarRot>%d and CarRot<%d?]",
00682           source, crstrt, clstrt, crstop, clstop, crstrt, crstop);
00683     else
00684       sprintf (recset_query,
00685           "%s[?(Car_Rot=%d and CRLN_OBS<=%.2f) or (Car_Rot=%d and CRLN_OBS>=%.2f)?]",
00686           source, crstrt, clstrt, crstop, clstop);
00687       }
00688       if (!(ids = drms_open_records (drms_env, recset_query, &status))) {
00689     fprintf (stderr, "Error: (%s) unable to open input data set %s\n",
00690           module_ident, recset_query);
00691     fprintf (stderr, "       status = %d\n", status);
00692     return 1;
00693       }
00694       if ((recct = ids->n) < 2) {
00695     fprintf (stderr, "Error: (%s) <2 records in selected input set\n",
00696         module_ident);
00697     fprintf (stderr, "       %s\n", recset_query);
00698     drms_close_records (ids, DRMS_FREE_RECORD);
00699     return 1;
00700       }
00701       free (inset);
00702       inset = strdup (recset_query);
00703       tmid = earth_meridian_crossing (clmid, crmid);
00704     } else {
00705                   /*  tmid specified as normal date-time string  */
00706       tmid = sscan_time (tmid_str);
00707       clmid = fp_nan;
00708       crmid = -1;
00709        sprintf (recset_query, "%s[?%s between %f and %f?]", source,
00710       primekey, tmid - 0.5 * intrvl, tmid + 0.5 * intrvl);
00711       if (!(ids = drms_open_records (drms_env, recset_query, &status))) {
00712     fprintf (stderr, "Error: (%s) unable to open input data set %s\n",
00713           module_ident, recset_query);
00714     fprintf (stderr, "       status = %d\n", status);
00715     return 1;
00716       }
00717       if ((recct = ids->n) < 2) {
00718     fprintf (stderr, "Error: (%s) <2 records in selected input set\n",
00719         module_ident);
00720     fprintf (stderr, "       %s\n", recset_query);
00721     drms_close_records (ids, DRMS_FREE_RECORD);
00722     return 1;
00723       }
00724       free (inset);
00725       inset = strdup (recset_query);
00726    }
00727   }
00728   if (verbose) printf ("processing %d input records\n", recct);
00729   propct = construct_stringlist (propagate_req, ',', &copykeylist);
00730                            /*  replace '+' with defaults  */
00731   if (propkeyct) {
00732     add_defaults = 0;
00733     for (n = 0; n < propct; n++) {
00734       if (!strcmp (copykeylist[n], "+")) {
00735     add_defaults = 1;
00736     strcpy (copykeylist[n], propagate[0]);
00737     n = propct;
00738       }
00739     }
00740     if (add_defaults) {
00741       int newct = propct + propkeyct - 1;
00742       copykeylist = (char **)realloc (copykeylist, newct * sizeof (char **));
00743       for (n = 1; n < propkeyct; n++) {
00744     copykeylist[n + propct - 1] = propagate[n];
00745       }
00746       propct = newct;
00747     }
00748   }
00749   if (verbose && propct) {
00750     printf ("propagating values for up to %d key(s):\n", propct);
00751     for (n = 0; n < propct; n++) {
00752       if (drms_keyword_lookup (orec, copykeylist[n], 0))
00753     printf ("  %s\n", copykeylist[n]);
00754       else
00755     printf ("  %s (not in output series)\n", copykeylist[n]);
00756     }
00757   }
00758 
00759   meanct = construct_stringlist (average_req, ',', &meankeylist);
00760                            /*  replace '+' with defaults  */
00761   if (meankeyct) {
00762     add_defaults = 0;
00763     for (n = 0; n < meanct; n++) {
00764       if (!strcmp (meankeylist[n], "+")) {
00765     add_defaults = 1;
00766     strcpy (meankeylist[n], average[0]);
00767     n = meanct;
00768       }
00769     }
00770     if (add_defaults) {
00771       int newct = meanct + meankeyct - 1;
00772       meankeylist = (char **)realloc (meankeylist, newct * sizeof (char **));
00773       for (n = 1; n < meankeyct; n++)
00774     meankeylist[n + meanct - 1] = average[n];
00775       meanct = newct;
00776     }
00777   }
00778   if (verbose) {
00779     printf ("averaging values for up to %d key(s):\n", meanct);
00780     for (n = 0; n < meanct; n++) {
00781       if (drms_keyword_lookup (orec, meankeylist[n], 0)) {
00782     printf ("  %s\n", meankeylist[n]);
00783     sprintf (keyname, "D_%s", meankeylist[n]);
00784     if (!drms_keyword_lookup (orec, keyname, 0))
00785       printf ("    Warning: %s not in output series\n", keyname);
00786       } else
00787     printf ("  %s (not in output series)\n", meankeylist[n]);
00788     }
00789   }
00790   avgval = (double *)calloc (meanct, sizeof (double));
00791   avgvalv = (double *)calloc (meanct, sizeof (double));
00792 
00793   irec = ids->records[0];
00794                          /*  establish input segment  */
00795   segct = drms_record_numsegments (irec);
00796   segnum = 0;
00797   if (segct > 1) {
00798     fprintf (stderr,
00799     "input records contain multiple segments, segment must be specified\n");
00800     drms_close_records (ids, DRMS_FREE_RECORD);
00801     return 0;
00802   }
00803   iseg = drms_segment_lookupnum (irec, segnum);
00804   naxis = iseg->info->naxis;
00805            /*  will dimensions of input segments have to be checked?  */
00806   checkseg = 0;
00807   if (iseg->info->scope == DRMS_VARDIM) {
00808     if (verbose) printf ("Warning: input segment sizes may vary\n");
00809     checkseg = 1;
00810     inaxis = (int *)malloc (naxis * sizeof (int));
00811   }
00812   ntot = 1;
00813   for (n = 0; n < naxis; n++) {
00814     ntot *= iseg->axis[n];
00815     if (checkseg) inaxis[n] = iseg->axis[n];
00816   }
00817                          /*  check for required keywords  */
00818   if (remove_obsvel) {
00819     keywd = drms_keyword_lookup (irec, "OBS_VR", 1);
00820     if (!keywd) {
00821       fprintf (stderr, "Warning: required keyword %s not found\n", "OBS_VR");
00822       fprintf (stderr, "         no correction applied for observer velocity\n");
00823       remove_obsvel = 0;
00824     }
00825   }
00826   if (filt_on_calver) {
00827     keywd = drms_keyword_lookup (irec, calverkey, 1);
00828     if (!keywd) {
00829       fprintf (stderr, "Warning: required keyword %s not found\n", calverkey);
00830       fprintf (stderr, "         no calibration version filtering applied\n");
00831       filt_on_calver = 0;
00832     }
00833   }
00834                        /*  set up WCS keys for averaging  */
00835   wcsaxes = drms_getkey_int (irec, "WCSAXES", &status);
00836   if (status) wcsaxes = naxis;
00837   crpix = (double *)calloc (wcsaxes, sizeof (double));
00838   crval = (double *)calloc (wcsaxes, sizeof (double));
00839   cdelt = (double *)calloc (wcsaxes, sizeof (double));
00840   crota = (double *)calloc (wcsaxes, sizeof (double));
00841   crpixv = (double *)calloc (wcsaxes, sizeof (double));
00842   crvalv = (double *)calloc (wcsaxes, sizeof (double));
00843   cdeltv = (double *)calloc (wcsaxes, sizeof (double));
00844   crotav = (double *)calloc (wcsaxes, sizeof (double));
00845                        /*  create new data maps for each of three statistics  */
00846   vval = (int *)calloc (ntot, sizeof (int));
00847   vavg = (double *)calloc (ntot, sizeof (double));
00848   if (needpowr) {
00849     vvar = (double *)calloc (ntot, sizeof (double));
00850     if (!vvar) {
00851       fprintf (stderr,
00852       "Error: unable to allocate %lld bytes required for output arrays\n",
00853       ntot * (sizeof (int) + 2 * sizeof (double)));
00854       drms_close_records (ids, DRMS_FREE_RECORD);
00855       return 1;
00856     }
00857   }
00858   if (!vval || !vavg) {
00859     fprintf (stderr,
00860     "Error: unable to allocate %lld bytes required for output arrays\n",
00861     ntot * (sizeof (int) + sizeof (double)));
00862     drms_close_records (ids, DRMS_FREE_RECORD);
00863     return 1;
00864   }
00865   vcts = drms_array_create (DRMS_TYPE_INT, naxis, iseg->axis, (void *)vval,
00866       &status);
00867   mean = drms_array_create (DRMS_TYPE_DOUBLE, naxis, iseg->axis, (void *)vavg,
00868       &status);
00869   if (needpowr) powr = drms_array_create (DRMS_TYPE_DOUBLE, naxis, iseg->axis,
00870       (void *)vvar, &status);
00871         /*  set the output scaling factors to be segment defaults
00872                     if not overridden as argument parameters  */
00873   if (mseg) {
00874     mean->bscale = (isnan (mscale) || mscale == 0.0) ? mseg->bscale : mscale;
00875     mean->bzero = (isnan (mzero)) ?  mseg->bzero : mzero;
00876   } else {
00877     mean->bscale = (isnan (mscale) || mscale == 0.0) ? 1.0 : mscale;
00878     mean->bzero = (isnan (mzero)) ?  0.0 : mzero;
00879   }
00880   if (needpowr) {
00881     if (pseg) {
00882       powr->bscale = (isnan (pscale) || pscale == 0.0) ? pseg->bscale : pscale;
00883       powr->bzero = (isnan (pzero)) ?  pseg->bzero : pzero;;
00884     } else {
00885       powr->bscale = (isnan (pscale) || pscale == 0.0) ? 1.0 : pscale;
00886       powr->bzero = (isnan (pzero)) ?  0.0 : pzero;;
00887     }
00888   }
00889            /*  determine range of representable values for scaled output  */
00890   if (mscaled) {
00891     unsigned long long maxval =
00892     (mseg->info->type == DRMS_TYPE_CHAR) ? SCHAR_MAX :
00893     (mseg->info->type == DRMS_TYPE_SHORT) ? SHRT_MAX :
00894     (mseg->info->type == DRMS_TYPE_INT) ? INT_MAX : LLONG_MAX;
00895     mrepmax = maxval;
00896     mrepmin = -mrepmax;
00897     mrepmax *= mean->bscale;
00898     mrepmin *= mean->bscale;
00899     mrepmax += mean->bzero;
00900     mrepmin += mean->bzero;
00901   }
00902   if (pscaled) {
00903     unsigned long long maxval =
00904     (pseg->info->type == DRMS_TYPE_CHAR) ? SCHAR_MAX :
00905     (pseg->info->type == DRMS_TYPE_SHORT) ? SHRT_MAX :
00906     (pseg->info->type == DRMS_TYPE_INT) ? INT_MAX : LLONG_MAX;
00907     prepmax = maxval;
00908     prepmin = -prepmax;
00909     prepmax *= powr->bscale;
00910     prepmin *= powr->bscale;
00911     prepmax += powr->bzero;
00912     prepmin += powr->bzero;
00913   }
00914 
00915   if (recct > maxct) {
00916     double bscale = 1.0;
00917     if (vseg) bscale = vseg->bscale;
00918     while (recct > maxct) {
00919       bscale *= 0.5;
00920       maxct *= 2;
00921     }
00922     vcts->bscale = bscale;
00923     fprintf (stderr,
00924         "Warning: more records in input data set than maximum valid count\n");
00925     if (vseg) fprintf (stderr,
00926         "BSCALE adjusted from %g to %g\n", vseg->bscale, bscale);
00927   }
00928      /*  write out some keys to make sure they're okay before proceeding  */
00929   kstat = 0;
00930   kstat += check_and_set_key_int   (orec, "CarrRot", crmid);
00931   kstat += check_and_set_key_float (orec, "CMLon", clmid);
00932   kstat += check_and_set_key_time  (orec, "MidTime", tmid);
00933   kstat += check_and_set_key_str   (orec, "Module", module_ident);
00934   kstat += check_and_set_key_str   (orec, "BLD_VERS", jsoc_version);
00935   kstat += check_and_set_key_str   (orec, "Source", source);
00936   kstat += check_and_set_key_str   (orec, "Input", inset);
00937   kstat += check_and_set_key_time  (orec, "Created", CURRENT_SYSTEM_TIME);
00938   kstat += check_and_set_key_float (orec, "Interval", intrvl);
00939   kstat += propagate_keys (orec, irec, copykeylist, propct);
00940                            /*  propagate WCS string keys  */
00941   for (n = 0; n < wcsaxes; n++) {
00942     sprintf (keyname, "CTYPE%d", n + 1);
00943     keystr = drms_getkey_string (irec, keyname, &status);
00944     if (!status) kstat += check_and_set_key_str (orec, keyname, keystr);
00945     sprintf (keyname, "CUNIT%d", n + 1);
00946     keystr = drms_getkey_string (irec, keyname, &status);
00947     if (!status) kstat += check_and_set_key_str (orec, keyname, keystr);
00948   }
00949   if (kstat) {
00950     fprintf (stderr, "Error writing key value(s) to %s\n", out_series);
00951     fprintf (stderr, "      output series may not have appropriate structure\n");
00952     if (!no_save) return 1;
00953   }
00954                /*  write argument list to Comment (or History)  */
00955 /*
00956   if (filt_on_calver) {
00957     calverinfo[0] = '\0';
00958     if (~cvaccept) {
00959       sprintf (calverinfo, "%s values accepted: %016llx", calverkey, cvaccept);
00960       if (cvreject) strcat (calverinfo, "\n");
00961     }
00962     if (cvreject) {
00963       sprintf (strbuf, "%s values rejected: %016llx", calverkey, cvreject);
00964       strcat (calverinfo, strbuf);
00965     }
00966   } else
00967     sprintf (calverinfo, "%s values accepted: ANY", calverkey);
00968 */
00969   if ((keywd = drms_keyword_lookup (orec, "COMMENT", 1))) {
00970     append_args_tokey (orec, "COMMENT");
00971 /*
00972     drms_appendhistory (orec, calverinfo, 1);
00973 */
00974   } else if ((keywd = drms_keyword_lookup (orec, "HISTORY", 1))) {
00975     append_args_tokey (orec, "HISTORY");
00976 /*
00977     drms_appendhistory (orec, calverinfo, 1);
00978 */
00979   }
00980           /*  support special hack of reading of rejection list file  */
00981   rejects = 0;
00982   if (strcmp (rejectfile, "Not Specified")) {
00983     FILE *rejectfp = fopen (rejectfile, "r");
00984     if (rejectfp) rejects = read_reject_list (rejectfp, &reject_list);
00985     else fprintf (stderr,
00986     "Warning: could not open rejection list %s; ignored\n", rejectfile);
00987   }
00988                       /*  initialize keys to be averaged  */
00989 /*
00990   tobs = tobsv = 0;
00991 */
00992   imgct = 0;
00993   crlnval = crlnvalv = 0.0;
00994                           /*  loop through input records  */
00995   for (rec = 0; rec < recct; rec++) {
00996     irec = ids->records[rec];
00997     if (verbose) report_pkey_value (stdout, irec, primekey, rec);
00998     if (writelog) report_pkey_value (runlog, irec, primekey, rec);
00999 
01000               /*  check for record quality, reject as applicable  */
01001     quality = drms_getkey_int (irec, qual_key, &status);
01002     if ((quality & qmask) && !status) {
01003       badqual++;
01004       if (verbose) printf ("skipped (quality)\n");
01005       if (writelog) fprintf (runlog, "skipped (quality)\n");
01006       continue;
01007     }
01008     if (rejects) {
01009                       /*  check against special rection list  */
01010       int idrec = drms_getkey_int (irec, "T_REC_index", &status);
01011       int match = 0;
01012       if (status) {
01013     fprintf (stderr, "Warning: \"T_REC_index\" keyword not found\n");
01014     fprintf (stderr, "         up to %d bad images could be processed\n",
01015       rejects);
01016     rejects = 0;
01017       }
01018       for (n = 0; n < rejects; n++) {
01019     if (idrec == reject_list[n]) {
01020       match = 1;
01021       break;
01022     }
01023       }
01024       if (match) {
01025         blacklist++;
01026         if (verbose) printf ("skipped (blacklist)\n");
01027         if (writelog) fprintf (runlog, "skipped (blacklist)\n");
01028         continue;
01029       }
01030     }
01031           /*  check for record calibration version; reject as applicable  */
01032     calver = drms_getkey_longlong (irec, calverkey, &status);
01033     if (filt_on_calver) {
01034       if (status) {
01035              /*  this probably never happens if the keyword exists  */
01036     if (verbose) printf ("skipped\n  (no %s key value)\n", calverkey);
01037     if (writelog) fprintf (runlog, "skipped\n  (no %s key value)\n",
01038         calverkey);
01039     continue;
01040       }
01041       for (cvct = 0; cvct < cvgoct; cvct++)
01042     if (calver == cvgolist[cvct]) break;
01043       if (cvct >= cvgoct) {
01044     badcv++;
01045     if (verbose)
01046       printf ("skipped\n  (calibration version key %s = %016llx)\n",
01047           calverkey, calver);
01048     if (writelog)
01049       fprintf (runlog, "skipped\n  (calibration version key %s = %016llx)\n",
01050           calverkey, calver);
01051     continue;
01052       }
01053       for (cvct = 0; cvct < cvnoct; cvct++) {
01054     if (calver == cvnolist[cvct]) {
01055       badcv++;
01056       if (verbose)
01057         printf ("skipped\n  (calibration version key %s = %016llx)\n",
01058             calverkey, calver);
01059       if (writelog)
01060         fprintf (runlog, "skipped\n  (calibration version key %s = %016llx)\n",
01061             calverkey, calver);
01062       continue;
01063     }
01064       }
01065     }
01066     for (cvct = 0; cvct < cvused; cvct++) {
01067       if (calver == cvlist[cvct]) {
01068         cvfound[cvct]++;
01069         break;
01070       }
01071     }
01072     if (cvct == cvused) {
01073       if (cvused < cvmaxct) {
01074         cvlist[cvct] = calver;
01075     cvfound[cvct] = 1;
01076     cvused++;
01077       }
01078     }
01079                   /*  check for non-nominal image rotation angle  */
01080     if (check_pa) {
01081       pa_rec = drms_getkey_float (irec, roll_key, &status);
01082       if (status && check_pa) {
01083     fprintf (stderr, "Warning: \"%s\" keyword not found\n", roll_key);
01084     fprintf (stderr, "         no limits on rotation angle\n");
01085     check_pa = 0;
01086     dpa_max = 360.0;
01087       }
01088       if (isfinite (pa_rec)) {
01089     dpa = fabs (pa_rec - pa_nom);
01090     while (dpa > 180.0) dpa -= 360.0;
01091     while (dpa < 0.0) dpa += 360.0;
01092     if (dpa > dpa_max) {
01093       badpa++;
01094       if (verbose) printf ("skipped (rotated)\n");
01095       if (writelog) fprintf (runlog, "skipped (rotated)\n");
01096       continue;
01097     }
01098       }
01099     }
01100                              /*  check segment structure  */
01101     iseg = drms_segment_lookupnum (irec, segnum);
01102     if (checkseg) {
01103       int okay = 1;
01104       for (n = 0; n < naxis; n++) {
01105     if (inaxis[n] != iseg->axis[n]) {
01106       okay = 0;
01107       break;
01108     }
01109       }
01110       if (!okay) {
01111         if (verbose) printf ("skipped (dimension mismatch)\n");
01112         if (writelog) fprintf (runlog, "skipped (dimension mismatch)\n");
01113         continue;
01114       }
01115     }
01116                                /*  read input data image  */
01117     data_array = drms_segment_read (iseg, DRMS_TYPE_DOUBLE, &status);
01118     if (status) {
01119       if (data_array) drms_free_array (data_array);
01120       if (verbose) printf ("skipped (segment read error)\n");
01121       if (writelog) fprintf (runlog, "skipped (segment read error)\n");
01122       else {
01123     fprintf (stderr, "Error reading file for record %d: ", rec);
01124     keywd = drms_keyword_lookup (irec, primekey, 1);
01125     if (!keywd) fprintf (stderr, "unknown value\n");
01126     else {
01127       drms_keyword_snprintfval (keywd, vbuf, sizeof (vbuf));
01128       fprintf (stderr, "%s\n", vbuf);
01129     }
01130       }
01131       continue;
01132     }
01133 
01134     tobs_rec = drms_getkey_time (irec, primekey, &status);
01135     if (status || time_is_invalid (tobs_rec)) {
01136       if (writelog) fprintf (runlog, "skipped (time invalid)\n");
01137       if (verbose) printf ("skipped (time invalid)\n");
01138       else fprintf (stderr, "error reading %s from record #%d\n", primekey, rec);
01139       if (data_array) drms_free_array (data_array);
01140       continue;
01141     }
01142     if (!imgct) tfirst = tobs_rec;
01143     tlast = tobs_rec;
01144 /*
01145     tobs += tobs_rec;
01146     tobsv += tobs_rec * tobs_rec;
01147 */
01148     imgct++;
01149 
01150     v = (double *)data_array->data;
01151     vval = (int *)vcts->data;
01152     vavg = (double *)mean->data;
01153     if (needpowr) vvar = (double *)powr->data;
01154     vobs = (remove_obsvel) ? drms_getkey_double (irec, "OBS_VR", &status) : 0.0;
01155                      /*  exponentiate logarithm as necessary  */
01156     log_base = drms_getkey_double (irec, "LOG_BASE", &status);
01157     if (!status && isfinite (log_base)) {
01158       log_base = log (log_base);
01159       for (n = 0; n < ntot; n++) v[n] = exp (log_base * v[n]);
01160       log_status = 1;
01161     } else log_status = 0;
01162                               /*  process valid data  */
01163     n = ntot;
01164     if (needpowr) {
01165       while (n--) {
01166     if (isfinite (*v)) {
01167       (*vval)++;
01168       *v -= vobs;
01169 /*
01170     if (function == DO_ABSVAL) *v = fabs (*v);
01171     if (function == DO_SQUARE) *v *= *v;
01172 */
01173           *vavg += *v;
01174           *vvar += *v * *v;
01175         }
01176     v++;
01177     vval++;
01178     vavg++;
01179     vvar++;
01180       }
01181     } else {
01182       while (n--) {
01183     if (isfinite (*v)) {
01184       (*vval)++;
01185       *v -= vobs;
01186 /*
01187     if (function == DO_ABSVAL) *v = fabs (*v);
01188     if (function == DO_SQUARE) *v *= *v;
01189 */
01190           *vavg += *v;
01191         }
01192     v++;
01193     vval++;
01194     vavg++;
01195       }
01196     }
01197     drms_free_array (data_array);
01198                           /*  get WCS keys for averaging  */
01199     for (n = 0; n < wcsaxes; n++) {
01200       sprintf (keyname, "CRPIX%d", n + 1);
01201       crpix_rec = drms_getkey_double (irec, keyname, &status);
01202       if (!status && isfinite (crpix_rec)) {
01203     crpix[n] += crpix_rec;
01204     crpixv[n] += crpix_rec * crpix_rec;
01205       }
01206       sprintf (keyname, "CRVAL%d", n + 1);
01207       crval_rec = drms_getkey_double (irec, keyname, &status);
01208       if (!status && isfinite (crval_rec)) {
01209     crval[n] += crval_rec;
01210     crvalv[n] += crval_rec * crval_rec;
01211       }
01212       sprintf (keyname, "CDELT%d", n + 1);
01213       cdelt_rec = drms_getkey_double (irec, keyname, &status);
01214       if (!status && isfinite (cdelt_rec)) {
01215     cdelt[n] += cdelt_rec;
01216     cdeltv[n] += cdelt_rec * cdelt_rec;
01217       }
01218       sprintf (keyname, "CROTA%d", n + 1);
01219       crota_rec = drms_getkey_double (irec, keyname, &status);
01220       if (!status && isfinite (crota_rec)) {
01221         crota_rec -= pa_nom;
01222     while (crota_rec < -180.0) crota_rec += 360.0;
01223     while (crota_rec > 180.0) crota_rec -= 360.0;
01224     crota[n] += crota_rec;
01225     crotav[n] += crota_rec * crota_rec;
01226       }
01227     }
01228                    /*  averaged CRLN_OBS and CarrRot/CAR_ROT  */
01229     crlnobs = drms_getkey_double (irec, "CRLN_OBS", &status);
01230     if (!status && isfinite (crlnobs)) {
01231       carrot = drms_getkey_int (irec, "CAR_ROT", &status);
01232       if (status || (carrot < 1)) {
01233     carrot = drms_getkey_int (irec, "CarrRot", &status);
01234     if (status || (carrot < 1)) carrot = 0;
01235       }
01236       crlnobs = 360.0 * carrot - crlnobs;
01237       crlnval += crlnobs;
01238       crlnvalv += crlnobs * crlnobs;
01239     }
01240                         /*  get requested keys for averaging  */
01241     for (n = 0; n < meanct; n++) {
01242       avg_rec = drms_getkey_double (irec, meankeylist[n], &status);
01243       if (!status && isfinite (avg_rec)) {
01244     avgval[n] += avg_rec;
01245     avgvalv[n] += avg_rec * avg_rec;
01246       }
01247     }
01248     if (verbose) printf ("ok\n");
01249     if (writelog) fprintf (runlog, "ok\n");
01250   }
01251   if (checkseg) free (inaxis);
01252   vval = (int *)vcts->data;
01253   vavg = (double *)mean->data;
01254   if (needpowr) {
01255     vvar = (double *)powr->data;
01256     for (n = 0; n < ntot; n++) {
01257       if (vval[n]) {
01258     vavg[n] /= vval[n];
01259     vvar[n] /= vval[n];
01260     vvar[n] -= vavg[n] * vavg[n];
01261       } else vavg[n] = vvar[n] = fp_nan;
01262     }
01263   } else {
01264     for (n = 0; n < ntot; n++) {
01265       if (vval[n]) vavg[n] /= vval[n];
01266       else vavg[n] = fp_nan;
01267     }
01268   }
01269   if (log_status) {
01270     double scale = 1.0 / log_base;
01271     if (needpowr) {
01272       for (n = 0; n < ntot; n++) {
01273     vavg[n] = scale * log (vavg[n]);
01274     vvar[n] = scale * log (vvar[n]);
01275       }
01276     } else {
01277       for (n = 0; n < ntot; n++) vavg[n] = scale * log (vavg[n]);
01278     }
01279   }
01280   if (vseg)
01281     if (drms_segment_write (vseg, vcts, 0))
01282       fprintf (stderr, "Warning: unable to write to count segment\n");
01283   if (mseg)
01284     if (drms_segment_write (mseg, mean, 0))
01285       fprintf (stderr, "Warning: unable to write to mean segment\n");
01286   if (pseg)
01287     if (drms_segment_write (pseg, powr, 0))
01288       fprintf (stderr, "Warning: unable to write to variance segment\n");
01289                               /*  set remaining keys  */
01290   kstat = 0;
01291   kstat += check_and_set_key_int  (orec, "DataRecs", imgct);
01292   kstat += check_and_set_key_int  (orec, "MissRecs", recct - imgct);
01293                              /*  statistics keys  */
01294   kstat += set_stats_keys (orec, vcts, mean, powr, ntot, mscaled, pscaled,
01295       mrepmin, mrepmax, prepmin, prepmax, verbose);
01296                                /*  averaged WCS keys  */
01297   if (imgct) {
01298     kstat += check_and_set_key_time  (orec, "T_FIRST", tfirst);
01299     kstat += check_and_set_key_time  (orec, "T_LAST", tlast);
01300 /*
01301     tobs /= imgct;
01302     kstat += check_and_set_key_time  (orec, primekey, tobs);
01303     tobsv /= imgct;
01304     tobsv -= tobs * tobs;
01305     sprintf (keyname, "D_%s", primekey);
01306     kstat += check_and_set_key_time  (orec, keyname, sqrt (tobsv));
01307 */
01308     for (n = 0; n < wcsaxes; n++) {
01309       crpix[n] /= imgct;
01310       sprintf (keyname, "CRPIX%d", n + 1);
01311       kstat += check_and_set_key_double (orec, keyname, crpix[n]);
01312       crpixv[n] /= imgct;
01313       crpixv[n] -= crpix[n] * crpix[n];
01314       sprintf (keyname, "D_CRPIX%d", n + 1);
01315       kstat += check_and_set_key_double (orec, keyname, sqrt (crpixv[n]));
01316       crval[n] /= imgct;
01317       sprintf (keyname, "CRVAL%d", n + 1);
01318       kstat += check_and_set_key_double  (orec, keyname, crval[n]);
01319       crvalv[n] /= imgct;
01320       crvalv[n] -= crval[n] * crval[n];
01321       sprintf (keyname, "D_CRVAL%d", n + 1);
01322       kstat += check_and_set_key_double (orec, keyname, sqrt (crvalv[n]));
01323       cdelt[n] /= imgct;
01324       sprintf (keyname, "CDELT%d", n + 1);
01325       kstat += check_and_set_key_double (orec, keyname, cdelt[n]);
01326       cdeltv[n] /= imgct;
01327       cdeltv[n] -= cdelt[n] * cdelt[n];
01328       sprintf (keyname, "D_CDELT%d", n + 1);
01329       kstat += check_and_set_key_double (orec, keyname, sqrt (cdeltv[n]));
01330       crota[n] /= imgct;
01331       crotav[n] /= imgct;
01332       crotav[n] -= crota[n] * crota[n];
01333       crota[n] += pa_nom;
01334       sprintf (keyname, "CROTA%d", n + 1);
01335       kstat += check_and_set_key_double (orec, keyname, crota[n]);
01336       sprintf (keyname, "D_CROTA%d", n + 1);
01337       kstat += check_and_set_key_double (orec, keyname, sqrt (crotav[n]));
01338     }
01339                  /*  averaged CRLN_OBS (and CarrRot/CAR_ROT)  */
01340     crlnval /= imgct;
01341     crlnobs = 360.0 - fmod (crlnval, 360.0);
01342     kstat += check_and_set_key_double (orec, "CRLN_OBS", crlnobs);
01343     crlnvalv /= imgct;
01344     crlnvalv -= crlnval * crlnval;
01345     kstat += check_and_set_key_double (orec, "D_CRLN_OBS", sqrt (crlnvalv));
01346                              /*  other averaged keys  */
01347     for (n = 0; n < meanct; n++) {
01348       avgval[n] /= imgct;
01349       kstat += check_and_set_key_double  (orec,  meankeylist[n], avgval[n]);
01350       avgvalv[n] /= imgct;
01351       avgvalv[n] -= avgval[n] * avgval[n];
01352       sprintf (keyname, "D_%s", meankeylist[n]);
01353       kstat += check_and_set_key_double (orec, keyname, sqrt (avgvalv[n]));
01354     }
01355     if (log_status) {
01356       sprintf (keyname, "LOG_BASE");
01357       kstat += check_and_set_key_double (orec, "LOG_BASE", exp (log_base));
01358     }
01359   }
01360   if (set_extra_key) {
01361     if (use_other_key) {
01362       spec_val = drms_getkey_double (orec, specuse_key, &status);
01363       if (status) {
01364         fprintf (stderr,
01365         "Warning: key %s does not exist in output series or is of wront type\n",
01366         specuse_key);
01367     fprintf (stderr, "         key %s will not be set\n", spec_key);
01368     set_extra_key = 0;
01369       }
01370     }
01371     kstat += check_and_set_key_double (orec, spec_key, spec_val);
01372   }
01373   if (kstat) {
01374     fprintf (stderr, "Error writing key value(s) to %s\n", out_series);
01375     fprintf (stderr, "      output series may not have appropriate structure\n");
01376     if (!no_save) return 1;
01377   }
01378                                     /*  clean up  */
01379   drms_close_records (ids, DRMS_FREE_RECORD);
01380   if (verbose) {
01381     printf ("record %s[:#%lld] ", out_series, orec->recnum);
01382     if (dispose == DRMS_FREE_RECORD) printf ("not ");
01383     printf ("written\n\n");
01384     if (badqual) printf
01385     ("    %d input records rejected for quality matching %08x\n",
01386     badqual, qmask);
01387     if (blacklist) printf
01388     ("    %d input records rejected from rejection list\n", blacklist);
01389     if (badpa) printf
01390     ("    %d input records rejected for roll difference exceeding %.2f\n",
01391     badpa, dpa_max);
01392     if (badcv) {
01393       printf ("    %d input records rejected for calib version matching %016llx\n",
01394       badcv, cvreject);
01395       printf ("                     or failing to match %016llx\n", cvaccept);
01396     }
01397     printf ("%s values used:", calverkey);
01398     for (cvct = 0; cvct < cvused; cvct++) {
01399       if (cvct) printf (",");
01400       printf (" %016llx (%d)", cvlist[cvct], cvfound[cvct]);
01401     }
01402     printf ("\n");
01403   }
01404   if (writelog) {
01405     if (badqual) fprintf (runlog,
01406     "    %d input records rejected for quality matching %08x\n",
01407     badqual, qmask);
01408     if (blacklist) fprintf (runlog,
01409     "    %d input records rejected from rejection list\n", blacklist);
01410     if (badpa) fprintf (runlog,
01411     "    %d input records rejected for roll difference exceeding %.2f\n",
01412     badpa, dpa_max);
01413     if (badcv) {
01414       fprintf (runlog,
01415       "    %d input records rejected for calib version matching %016llx\n",
01416       badcv, cvreject);
01417       fprintf (runlog,
01418       "                     or failing to match %016llx\n", cvaccept);
01419     }
01420     fprintf (runlog, "%s values used:", calverkey);
01421     for (cvct = 0; cvct < cvused; cvct++) {
01422       if (cvct) fprintf (runlog, ",");
01423       fprintf (runlog, " %016llx (%d)", cvlist[cvct], cvfound[cvct]);
01424     }
01425     fprintf (runlog, "\n");
01426     fclose (runlog);
01427   }
01428   drms_close_record (orec, dispose);
01429   return 0;
01430 }
01431 
01432 /*
01433  *  Revision history
01434  *  (all mods by Rick Bogart unless otherwise indicated)
01435  *
01436  *  v 0.0  09.11.01 created file, based on version 5.0 of SOI module
01437  *      doppavg (CM version of 2004.10.25)
01438  *  v 0.1  09.12.01 include a certain minimal amount of input processing
01439  *  v 0.2  10.01.12 include minimal processing and output to DRMS
01440  *  v 0.3  10.01.13 checks on input size
01441  *  v 0.4  10.02.09 include time selection by target CR:CL and rot width;
01442  *      setting of a few keys
01443  *  v 0.5   some cleanup, added setting of a few more keys, quality
01444  *      checking from keyword and blacklist, stubs for propagation of
01445  *      fixed keys and averaging of varying keys
01446  *  v 0.6   Stripped commented code from doppavg; see v05 for them
01447  *      Added support for key propagation, output keys MidTime,
01448  *        averaged, WCS keys, statistics keys
01449  *      Added optional correction for observer velocity (radial only)
01450  *      Moved construct_stringlist() function to keystuff
01451  *      Added trap for bad segment read
01452  *      Added option for no output
01453  *      Added diagnostics at end
01454  *  v 0.6 frozen 10.08.09
01455  *  v 0.7
01456  *  10.08.11    Removed default values for in/out arguments and
01457  *        unspecified defaults for input set specifers
01458  *      Fixed bug that was preventing HMI rejection lists from being
01459  *        read properly
01460  *  10.08.18    CRVALn, CRPIXn, CDELTn, & CROTn always averaged with
01461  *        standard deviations; CTYPEn and CUNITn always propagated;
01462  *        generalized to arbitrary number of WCS axes
01463  *      Implemented expansion of keyword requests for copied and
01464  *        averaged keywords to include defaults if requested
01465  *      Implemented averaging of requested keywords
01466  *  v 0.7 frozen 10.08.19
01467  *  v 0.8
01468  *  10.09.13    S Chakraborty added code to support log-base data
01469  *  10.11.16    added test for acceptable roll-angles, with defaults
01470  *      appropriate to HMI
01471  *  v 0.8 frozen 10.11.16
01472  *  v 0.9
01473  *  11.01.11    added reporting of locations of extremal values, checks
01474  *      for out-of-range values
01475  *  11.02.10    added traps for inappropriate output series structure
01476  *  v 0.9 frozen 11.02.28
01477  *  v 1.0
01478  *  11.04.23    changed argument tobs_key to pkey, generalized treatment
01479  *      of prime key to different types
01480  *  11.06.06    added keywords for output segment names to override the
01481  *      default values
01482  *  11.06.17    allow for output series without all three segments
01483  *  11.06.21    added setting of T_FIRST, T_LAST keys
01484  *  11.06.22    added averaging of CRLN_OBS
01485  *  v 1.0 frozen 11.11.14
01486  *  12.05.24    fixed averaging of CROTA* for cases where values may
01487  *      wrap; SAT_ROT values no longer averaged by default; provide
01488  *      option for overriding roll angle keyword; added setkey, setval
01489  *      options
01490  *  v 1.1 frozen 12.08.03
01491  *  12.09.12    added support for acceptance and/or rejection of
01492  *      records with certain values of CalVer64 (or other equivalent
01493  *      key)
01494  *  12.10.04    added recording of calling params info to comment or
01495  *      history keyword
01496  *  12.10.16    added writing of verbose output to log segment if
01497  *      present, and recording of input recnum's; added logging of
01498  *      CalVer64 values used; added support for specification of
01499  *      tmid in date_time format
01500  *  v 1.2 frozen 13.04.08
01501  *  14.03.01    added T_OBS to default list of keys to be averaged;
01502  *      added verbose warnings of missing output keys for propagation
01503  *      or averaging; removed automatic averaging of prime key
01504  *  v 1.3 frozen 14.07.07
01505  *  14.10.09    removed functions drms_appendstr_tokey() and
01506  *      append_args_tokey() (now in keystuff)
01507  *  v 1.4 frozen 15.01.13
01508  *  15.02.24    added initializations of several variables; made
01509  *      calculation of power values conditional on presence of the
01510  *      appropriate segment in the output record
01511  *  v 1.5 frozen 15.04.08
01512  */