JSOC_Documentation: proj/globalhs/apps/jtsfiddle.c Source File

00001 /* this JSOC module is a port of an SOI module written by Rasmus Larsen.
00002  * ported by Tim Larson.
00003  * tim doubts this works correctly with nskip or navg
00004  */
00005 
00006 /* 
00007  *  Detrending/gapfilling/differencing module
00008  *  ts_fiddle_rml upgrades ts_fiddle
00009  */
00010 
00011 #include <stdio.h>
00012 #include <stdlib.h>
00013 #include <strings.h>
00014 #include <time.h>
00015 #include <sys/time.h>
00016 #include <sys/resource.h>
00017 #include <math.h>
00018 #include <assert.h>
00019 #include <fftw3.h>
00020 
00021 #include "jsoc_main.h"
00022 #include "fitsio.h"
00023 
00024 #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
00025 #define kNOTSPECIFIED "not specified"
00026 
00027 char *module_name = "jtsfiddle";
00028 char *cvsinfo_jtsfiddle = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.15 2013/07/02 20:33:16 tplarson Exp $";
00029 
00030 /* Command line arguments: */
00031 ModuleArgs_t module_args[] =
00032 {
00033   {ARG_STRING,   "in", "", "input data records"},
00034   {ARG_STRING,   "tsout", kNOTSPECIFIED, "output dataseries for timeseries"},
00035   {ARG_STRING,   "powout", kNOTSPECIFIED, "output dataseries for power spectra"},
00036   {ARG_STRING,   "fftout", kNOTSPECIFIED, "output dataseries for fft's"},
00037   {ARG_STRING,   "fft1out", kNOTSPECIFIED, "output dataseries for fft's reordered"},
00038   {ARG_STRING,   "arpowout", kNOTSPECIFIED, "output dataseries for AR power"},
00039   {ARG_STRING,   "mavgout", kNOTSPECIFIED, "output dataseries for m-averaged power spectra"},
00040 //  {ARG_STRING,   "out", "", "output data series"},
00041   {ARG_STRING,   "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"},
00042   {ARG_STRING,   "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
00043   {ARG_STRING,   "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
00044   {ARG_STRING,   "srclink",  "SRCDATA", "name of link to source data"},
00045   {ARG_INT,      "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
00046   {ARG_INT,      "VERB", "1",  "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL},
00047   {ARG_STRING,   "TAG", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, usually used as an additional primekey"},
00048   {ARG_STRING,   "VERSION", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, useful for selecting obsolete records"},
00049 
00050   {ARG_STRING,   "LOGFILE", "stdout", ""},
00051   {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},
00052   {ARG_STRING,   "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"},
00053 
00054   {ARG_INT,      "NDT", "-1", "", ""},
00055 //  {ARG_FLOAT,    "TSAMPLE", "60", "", ""},
00056   {ARG_INT,      "IFILL", "1", "", ""},
00057   {ARG_INT,      "IFIRST", "0", "", ""},
00058 //  {ARG_INT,      "FLIPM", "0", "", ""},
00059   {ARG_INT,      "MAX_AR_ORDER", "360", "", ""},
00060   {ARG_INT,      "FU_FIT_EDGES", "1", "", ""},
00061   {ARG_INT,      "FU_ITER", "3", "", ""},
00062   {ARG_INT,      "FU_MIN_PERCENT", "90", "", ""},
00063   {ARG_FLOAT,    "CDETREND", "50.0", "", ""},
00064   {ARG_INT,      "DETREND_LENGTH", "1600", "", ""},
00065   {ARG_INT,      "DETREND_SKIP", "1440", "", ""},
00066   {ARG_INT,      "DETREND_FIRST", "0", "", ""},
00067 /*
00068   {ARG_INT,      "TSOUT", "0", "", ""},
00069   {ARG_INT,      "FFTOUT", "0", "", ""},
00070   {ARG_INT,      "FFT1OUT", "0", "", ""},
00071   {ARG_INT,      "POWOUT", "0", "", ""},
00072   {ARG_INT,      "ARPOWOUT", "0", "", ""},
00073   {ARG_INT,      "MAVGOUT", "0", "", ""},
00074 */
00075   {ARG_INT,      "NAVG", "0", "", ""},
00076   {ARG_INT,      "NSKIP", "0", "", ""},
00077   {ARG_INT,      "NOFF", "0", "", ""},
00078   {ARG_INT,      "IMIN", "0", "", ""},
00079   {ARG_INT,      "IMAX", "-1", "", ""},
00080   {ARG_END},
00081 };
00082 
00083 #include "saveparm.c"
00084 #include "timing.c"
00085 #include "set_history.c"
00086 
00087 extern void cdetrend_discontig( int n, _Complex float *data, int *isgood, 
00088                 int degree, int length, int skip, 
00089                 int m, int *sect_last, int detrend_first);
00090 
00091 int cfahlman_ulrych(int n, _Complex float *data, int *isgood, 
00092             int minpercentage, int maxorder, int iterations, 
00093             int padends, int *order, _Complex float *ar_coeff);
00094 
00095 //char *getdetrendversion(void);
00096 //char *getgapfillversion(void);
00097 
00098 /* global variables holding the values of the command line variables. */
00099 static char *logfile, *gapfile, *sectionfile;
00100 static int lmode, n_sampled_out, ifill, flipm;
00101 static int ar_order, max_ar_order, fu_iter, fu_fit_edges;
00102 static int fu_min_percent;
00103 static int detrend_length, detrend_skip, detrend_first;
00104 static float tsample, detrend_const;
00105 static int ifirst, tsout,fftout, fft1out, powout, arpowout, mavgout;
00106 static int navg, nskip, noff, imin, imax;
00107 
00108 
00109 /* Prototypes for local functions. */
00110 static double splitting(int l, int m);
00111 
00112 /* Prototypes for external functions */
00113 extern void cffti_(int *, float *);
00114 extern void cfftb_(int *, float *, float *);
00115 static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
00116           float tsample_in, float *tsample_out,
00117           int *num_sections, int *sections);
00118 static void extract_data(int n_in, int *gaps, float *data_in, float *data_out);
00119 static int read_section_file(char *filename, int *num_section, int **section);
00120 
00121 /************ Here begins the main module **************/
00122 int DoIt(void)
00123 {
00124   int i;
00125   int m, n_in, n_sampled_in;
00126   TIME epoch, step, start, stop;
00127   int gapsize, istart, istop, data_dim, detrend_order;
00128   int npow, mshift;
00129   float tsamplex, df1, c;
00130   int *agaps;
00131   int *gaps, *gaps2;
00132   float *msum, *pow;
00133   float *data, *out_data, *in_data;
00134   fftwf_plan plan;
00135   int num_sections, *sections;
00136   float *ar_coeff=NULL;
00137   float *datarr, *datptr, *outptr;
00138   float *datahold;
00139   char tstartstr[100];
00140   int lmin, lmax;
00141 
00142   int fstatus = 0;
00143   fitsfile *fitsptr;
00144   long fdims[1];
00145   long fpixel[]={1};
00146   int *anynul=0;
00147 
00148   int instart[2], inend[2];
00149   int tslength[2], tsstart[2], tsend[2], tstotal[2];
00150   int fft1length[2], fft1start[2], fft1end[2], fft1total[2];
00151   int powlength[2], powstart[2], powend[2], powtotal[2];
00152   int status=DRMS_SUCCESS;
00153   int newstat=0;
00154 
00155   DRMS_Segment_t *segin = NULL;
00156   DRMS_Segment_t *segout = NULL;
00157   DRMS_Array_t *inarr = NULL;
00158   DRMS_Array_t *outarr = NULL;
00159   DRMS_Record_t *inrec = NULL;
00160   DRMS_Record_t *outrec = NULL;
00161   DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
00162   DRMS_RecLifetime_t lifetime;
00163   long long histrecnum = -1;
00164   DRMS_Segment_t *gapseg = NULL;
00165   DRMS_Array_t *gaparr = NULL;
00166 
00167   HIterator_t outKey_list;
00168   DRMS_Keyword_t *outKey;
00169   TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
00170 
00171   int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
00172   int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
00173 
00174   double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout;
00175   double wt0, wt1, wt2, wt3, wt;
00176   double ut0, ut1, ut2, ut3, ut;
00177   double st0, st1, st2, st3, st;
00178   double ct0, ct1, ct2, ct3, ct;
00179 
00180   wt0=getwalltime();
00181   ct0=getcputime(&ut0, &st0);
00182 
00183   /* Read input parameters. */
00184   logfile     = (char *)cmdparams_save_str    (&cmdparams, "LOGFILE", &newstat);
00185   gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
00186   sectionfile = (char *)cmdparams_save_str    (&cmdparams, "SECTIONFILE", &newstat);
00187   n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);
00188 //  tsample     = cmdparams_save_float  (&cmdparams, "TSAMPLE", &newstat);
00189   ifill       = cmdparams_save_int    (&cmdparams, "IFILL", &newstat);
00190   max_ar_order= cmdparams_save_int    (&cmdparams, "MAX_AR_ORDER", &newstat);
00191   fu_fit_edges= cmdparams_save_int    (&cmdparams, "FU_FIT_EDGES", &newstat);
00192   fu_iter     = cmdparams_save_int    (&cmdparams, "FU_ITER", &newstat);
00193   fu_min_percent = cmdparams_save_int (&cmdparams, "FU_MIN_PERCENT", &newstat);
00194   ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
00195 //  flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);
00196   detrend_const = cmdparams_save_float(&cmdparams, "CDETREND", &newstat);
00197   detrend_length = cmdparams_save_int (&cmdparams, "DETREND_LENGTH", &newstat);
00198   detrend_skip = cmdparams_save_int   (&cmdparams, "DETREND_SKIP", &newstat);
00199   detrend_first = cmdparams_save_int  (&cmdparams, "DETREND_FIRST", &newstat);
00200 /*
00201   tsout       = cmdparams_save_int    (&cmdparams, "TSOUT", &newstat);
00202   fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);
00203   fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);
00204   powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);
00205   arpowout    = cmdparams_save_int    (&cmdparams, "ARPOWOUT", &newstat);
00206   mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);
00207 */
00208   navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);
00209   nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);
00210   noff        = cmdparams_save_int    (&cmdparams, "NOFF", &newstat);
00211   imin        = cmdparams_save_int    (&cmdparams, "IMIN", &newstat);
00212   imax        = cmdparams_save_int    (&cmdparams, "IMAX", &newstat);
00213 
00214   char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
00215 //  char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
00216   char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
00217   char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
00218   int seginflag = strcmp(kNOTSPECIFIED, segnamein);
00219   int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
00220   char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
00221   char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
00222   int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
00223   int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
00224   if (permflag)
00225     lifetime = DRMS_PERMANENT;
00226   else
00227     lifetime = DRMS_TRANSIENT;
00228 
00229   char *tag = (char *)cmdparams_save_str(&cmdparams, "TAG", &newstat);
00230   int tagflag = strcmp(kNOTSPECIFIED, tag);
00231   char *version = (char *)cmdparams_save_str(&cmdparams, "VERSION", &newstat);
00232   int verflag = strcmp(kNOTSPECIFIED, version);
00233 
00234   char *tsseries = (char *)cmdparams_save_str(&cmdparams, "tsout", &newstat);
00235   tsout = strcmp(kNOTSPECIFIED, tsseries);
00236   char *powseries = (char *)cmdparams_save_str(&cmdparams, "powout", &newstat);
00237   powout = strcmp(kNOTSPECIFIED, powseries);
00238   char *fftseries = (char *)cmdparams_save_str(&cmdparams, "fftout", &newstat);
00239   fftout = strcmp(kNOTSPECIFIED, fftseries);
00240   char *fft1series = (char *)cmdparams_save_str(&cmdparams, "fft1out", &newstat);
00241   fft1out = strcmp(kNOTSPECIFIED, fft1series);
00242   char *arpowseries = (char *)cmdparams_save_str(&cmdparams, "arpowout", &newstat);
00243   arpowout = strcmp(kNOTSPECIFIED, arpowseries);
00244   char *mavgseries = (char *)cmdparams_save_str(&cmdparams, "mavgout", &newstat);
00245   mavgout = strcmp(kNOTSPECIFIED, mavgseries);
00246 
00247   char *serieslist[6];
00248   enum seriesindex  {TSOUT, FFTOUT, FFT1OUT, POWOUT, ARPOWOUT, MAVGOUT};
00249   int  flagarr[6] = {tsout, fftout, fft1out, powout, arpowout, mavgout};
00250   int mfliparr[6] = {0, 0, 0, 0, 0, 0};
00251   int msignarr[6] = {0, 0, 0, 0, 0, 0};
00252   serieslist[TSOUT]=tsseries;
00253   serieslist[FFTOUT]=fftseries;
00254   serieslist[FFT1OUT]=fft1series;
00255   serieslist[POWOUT]=powseries;
00256   serieslist[ARPOWOUT]=arpowseries;
00257   serieslist[MAVGOUT]=mavgseries;
00258   long long histrecnumarr[6]={-1, -1, -1, -1, -1, -1};
00259   DRMS_RecordSet_t *outrecsetarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};
00260   DRMS_Record_t       *outrecarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};
00261   DRMS_Segment_t      *segoutarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};
00262 
00263   /**** Sanity check of input parameters. ****/
00264 
00265   if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0)) 
00266   {
00267     fprintf(stderr, "ERROR: must specify an output dataseries\n");
00268     return 1;
00269   }
00270 
00271   if (navg <= 0) 
00272     navg=1;
00273   if (nskip <= 0) 
00274     nskip=1;
00275   if ((navg != 1) && (nskip != 1))
00276   {
00277     fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");
00278     return 1;
00279   }
00280   if (noff < 0 || noff >= nskip) 
00281   {
00282     fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
00283     return 1;
00284   }
00285   if (arpowout && !ifill)
00286   {
00287     fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");
00288     return 1;
00289   }
00290   if (fu_min_percent <= 0 || fu_min_percent > 100)
00291   {
00292     fprintf(stderr, "ERROR: FU_MIN_PERCENT must be > 0 and <= 100\n");
00293     return 1;
00294   }
00295 
00296   if (newstat) 
00297   {
00298     fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
00299     cpsave_decode_error(newstat);
00300     return 1;
00301   }  
00302   else if (savestrlen != strlen(savestr)) 
00303   {
00304     fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
00305     return 1;
00306   }
00307 
00308 
00309  // these must be present in the output dataseries and variable, not links or constants
00310  // the loop is only necessary if the output series don't all link to the same history series, which they usually will, but just in case...
00311   char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
00312   int is, ishold, itest, mflip;
00313   char *holdseries="";
00314 /*
00315   char *cvsinfo;
00316   cvsinfo = (char *)malloc(1024);
00317   strcpy(cvsinfo,"$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.15 2013/07/02 20:33:16 tplarson Exp $");
00318   strcat(cvsinfo,"\n");
00319   strcat(cvsinfo,getdetrendversion());
00320   strcat(cvsinfo,"\n");
00321   strcat(cvsinfo,getgapfillversion());
00322 */
00323   for (is=0;is<6;is++)
00324   {
00325 
00326     if (!flagarr[is])
00327       continue;
00328 
00329     DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
00330                                                    serieslist[is],
00331                                                    DRMS_TRANSIENT, 
00332                                                    &status);
00333 
00334     if (status != DRMS_SUCCESS) 
00335     {
00336      fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
00337      return 1;
00338     }
00339 
00340     DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
00341 
00342 
00343     // set up ancillary dataseries for processing metadata
00344     if (histlink != NULL) 
00345     {
00346       if (strcmp(holdseries, histlink->info->target_series))
00347       {
00348         ishold=is;
00349         holdseries=strdup(histlink->info->target_series);
00350         histrecnumarr[is]=set_history(histlink);
00351         if (histrecnumarr[is] < 0)
00352         {
00353           fprintf(stderr, "ERROR: problem creating record in history dataseries for output dataseries %s\n", serieslist[is]);
00354           drms_close_record(tempoutrec, DRMS_FREE_RECORD);
00355           return 1;
00356         }
00357       }
00358       else
00359       {
00360         histrecnumarr[is]=histrecnumarr[ishold];
00361       }
00362     }
00363     else
00364     {
00365       fprintf(stderr,"WARNING: could not find history link in output dataseries %s\n", serieslist[is]);
00366     }
00367 
00368     for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
00369     {
00370       DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
00371       if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
00372       {
00373         fprintf(stderr, "ERROR: output keyword %s in series %s is either missing, constant, or a link\n", outchecklist[itest], serieslist[is]);
00374         drms_close_record(tempoutrec, DRMS_FREE_RECORD);
00375         return 1;
00376       }
00377     }
00378 
00379     mflip = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
00380     if (status == DRMS_SUCCESS)
00381       mfliparr[is]=mflip;    
00382 
00383     drms_close_record(tempoutrec, DRMS_FREE_RECORD);
00384   }
00385 
00386 // for efficiency require that tsseries and fftseries have the same value of MFLIPPED.  this restriction may be lifted as noted below.
00387   if ((tsout && fftout) && (mfliparr[TSOUT] != mfliparr[FFTOUT]))
00388   {
00389     fprintf(stderr, "ERROR: series %s and %s must have the same value of MFLIPPED\n", tsseries, fftseries);
00390     return 1;
00391   }
00392 
00393   if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
00394   {
00395     DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status);
00396     if (status != DRMS_SUCCESS || gaprecset == NULL) 
00397     {
00398       fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
00399       return 1;
00400     }
00401     if (gaprecset->n == 0)
00402     {
00403       fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
00404       return 1;
00405     }
00406     if (gaprecset->n > 1)
00407     {
00408       fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
00409       return 1;
00410     } 
00411     gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);
00412     if (gapseg != NULL)
00413       gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status);
00414     if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
00415     {
00416       fprintf(stderr, "problem reading gaps segment: status = %d\n", status);
00417       return 1;
00418     }
00419     agaps=(int *)(gaparr->data);
00420     gapsize=gaparr->axis[0];
00421   }
00422   else
00423   {
00424     fits_get_img_size(fitsptr, 1, fdims, &fstatus);
00425     gapsize=fdims[0];
00426     agaps=(int *)(malloc(gapsize*sizeof(int)));
00427     fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus); 
00428     fits_close_file(fitsptr, &fstatus);
00429     if (fstatus) 
00430     {
00431       fits_report_error(stderr, fstatus);
00432       return 1;
00433     }
00434   }
00435 
00436   if (verbflag)
00437     printf("gapfile read, gapsize = %d\n",gapsize);
00438 
00439   data_dim=gapsize;
00440   if (n_sampled_out>gapsize) 
00441     data_dim=n_sampled_out;
00442 
00443   /* allocate temporary storage */
00444   gaps=(int *)(malloc(gapsize*sizeof(int)));
00445   gaps2=(int *)(malloc(gapsize*sizeof(int)));
00446   data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
00447   ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float));
00448 // now done below
00449 //  if (fft1out || powout || mavgout || arpowout)
00450 //    pow=(float *)(malloc(n_sampled_out*sizeof(float)));
00451 
00452   /* Read location of jumps. */
00453   if (!strcmp(sectionfile,kNOTSPECIFIED))
00454   {
00455     /* No sections file specified. Assume that the whole 
00456        time series in one section. */
00457     num_sections=1;
00458     sections = malloc(2*sizeof(int));
00459     sections[0] = 0;
00460     sections[1] = data_dim-1;
00461   }
00462   else
00463   {
00464     int secstat;
00465     if (secstat=read_section_file(sectionfile, &num_sections, &sections))
00466     {
00467       DRMS_RecordSet_t *secrecset = drms_open_records(drms_env, sectionfile, &status);
00468       if (status != DRMS_SUCCESS || secrecset == NULL) 
00469       {
00470         fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status);
00471         return 1;
00472       }
00473       if (secrecset->n == 0)
00474       {
00475         fprintf(stderr,"ERROR: sectionfile recordset contains no records\n");
00476         return 1;
00477       }
00478       if (secrecset->n > 1)
00479       {
00480         fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n");
00481         return 1;
00482       }
00483       num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL);
00484       if (num_sections < 1)
00485       {
00486         fprintf(stderr,"ERROR: Invalid NSECS keywords\n");
00487         return 1;
00488       }
00489       sections = (int *)malloc(2*(num_sections)*sizeof(int));
00490       char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL);
00491       sscanf(strtok(sectkey," \n"),"%d",&(sections[0]));
00492       sscanf(strtok(NULL," \n"),"%d",&(sections[1]));
00493       int i;
00494       for (i=2 ;i<2*(num_sections); i+=2)
00495       {
00496         sscanf(strtok(NULL," \n"), "%d", &(sections)[i]);
00497         sscanf(strtok(NULL," \n"), "%d", &(sections)[i+1]);
00498 
00499         if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1]))
00500         {
00501           fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n");
00502           return 1;
00503         }
00504       }
00505       free(sectkey);
00506     }
00507   }
00508 
00509   if (verbflag)
00510     printf("num_sections = %d\n",num_sections);
00511 
00512   DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);
00513 
00514   if (status != DRMS_SUCCESS || inrecset == NULL) 
00515   {
00516     fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
00517     return 1;
00518   }
00519   if (inrecset->n == 0)
00520   {
00521     fprintf(stderr, "ERROR: input recordset contains no records\n");
00522     return 1;
00523   }
00524 
00525   inrec=inrecset->records[0];
00526   char *inchecklist[] = {"T_START", "T_STOP", "QUALITY", "LMIN", "LMAX", "T_STEP"};
00527   for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
00528   {
00529     DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
00530     if (inkeytest == NULL)
00531     {
00532       fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
00533       drms_close_records(inrecset, DRMS_FREE_RECORD);
00534       return 1;
00535     }
00536   }
00537 
00538   cadence=drms_getkey_float(inrec, "T_STEP", NULL); //assume constant across all input records
00539   tsample=cadence;
00540   tstartsave=drms_getkey_time(inrec, "T_START", NULL); //must be constant across all input records unless GAPFILE is implemented differently
00541 
00542   int mflipin=drms_getkey_int(inrec, "MFLIPPED", &status);
00543   if (status != DRMS_SUCCESS) mflipin=0;
00544 
00545 // already required mfliparr[TSOUT] == mfliparr[FFTOUT] above
00546   if (tsout && mfliparr[TSOUT] != mflipin)
00547     flipm=1;
00548   else if (fftout && mfliparr[FFTOUT] != mflipin)
00549     flipm=1;
00550   else
00551     flipm=0;
00552 
00553   for (is=2;is<5;is++)
00554   {
00555     if (!flagarr[is])
00556       continue;
00557     if (mfliparr[is] != (mflipin || flipm))
00558       msignarr[is]=-1;
00559     else
00560       msignarr[is]=1;
00561   }
00562 
00563   if (mflipin)
00564     msignarr[MAVGOUT]=1;
00565   else
00566     msignarr[MAVGOUT]=-1;
00567 
00568   // changed n_in to gapsize in following call
00569   if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
00570   {
00571     fprintf(stderr, "ERROR: problem in extract_gaps\n");
00572     return 1;
00573   }
00574 
00575   /* Adjust detrend parameters according to the number of
00576      available points. */
00577   if (n_sampled_in<detrend_length)
00578   {
00579     detrend_length = n_sampled_in;
00580     detrend_skip = detrend_length;
00581   }
00582 
00583   int idtf;
00584   if (detrend_first > 0)
00585   {
00586     for (idtf=0;idtf<detrend_first;idtf++)
00587       gaps[idtf]=0;
00588   }
00589 
00590 
00591   if (n_sampled_out == -1) 
00592     n_sampled_out = n_sampled_in;
00593   df1 = tsamplex*n_sampled_out;
00594 
00595   if (fftout || fft1out || powout || arpowout || mavgout)     
00596     plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, 
00597                 (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);
00598 
00599   /* Set sum to zero before starting */
00600   if (mavgout) 
00601     msum = (float *)malloc((n_sampled_out/2)*sizeof(float));
00602 
00603   if (fft1out || powout || mavgout || arpowout)
00604   {
00605     /* Set first and last output index if not set as a parameter. */
00606     if (imax < imin) 
00607     {
00608       imin=0;
00609       imax=n_sampled_out/2-1;
00610     }
00611     npow=imax-imin+1;
00612     pow=(float *)(malloc(n_sampled_out*sizeof(float)));
00613     datahold=(float *)malloc(2*npow*sizeof(float));
00614   }
00615 
00616 // needed to implement mfliparr[TSOUT] != mfliparr[FFTOUT] 
00617 //  float *datatemp=(float *)malloc(2*n_sampled_out*sizeof(float));
00618   float *tmpptr=data;
00619 
00620   if (tsout || fftout) 
00621   {
00622     tslength[0]=tstotal[0]=2*n_sampled_out;
00623     tslength[1]=1;
00624     tsstart[0]=0;
00625     tsend[0]=2*n_sampled_out-1;
00626   }
00627   if (fft1out)
00628   {
00629     fft1length[0]=fft1total[0]=2*npow;
00630     fft1length[1]=1;
00631     fft1start[0]=0;
00632     fft1end[0]=2*npow-1;
00633   }
00634   if (powout || arpowout || mavgout)
00635   {
00636     powlength[0]=powtotal[0]=npow;
00637     powlength[1]=1;
00638     powstart[0]=0;
00639     powend[0]=npow-1;
00640   }
00641   instart[0]=0;
00642   inend[0]=2*gapsize - 1;
00643 
00644   status=drms_stage_records(inrecset, 1, 0);
00645   if (status != DRMS_SUCCESS)
00646   {
00647     fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status);
00648     return 1;
00649   }
00650 
00651   for (is=0; is<6; is++)
00652   {
00653     if (!flagarr[is])
00654       continue;
00655 
00656     outrecsetarr[is] = drms_create_records(drms_env, inrecset->n, serieslist[is], DRMS_PERMANENT, &status);
00657 
00658     if (status != DRMS_SUCCESS || outrecsetarr[is] == NULL) 
00659     {
00660      fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
00661      return 1;
00662     }
00663   }
00664 
00665   int irec;
00666   for (irec=0;irec<inrecset->n;irec++)
00667   {
00668 
00669     inrec = inrecset->records[irec];
00670 
00671     lmin=drms_getkey_int(inrec, "LMIN", NULL);
00672     lmax=drms_getkey_int(inrec, "LMAX", NULL);
00673     if (lmin != lmax)
00674     {
00675       fprintf(stderr, "ERROR: lmin != lmax not yet supported.\n");
00676       return 0;
00677     }
00678     lmode=lmin;
00679 
00680     if (verbflag > 1) 
00681     {
00682       printf("starting l=%d\n", lmode);
00683     }
00684 
00685     tstart=drms_getkey_time(inrec, "T_START", NULL);
00686     tstop =drms_getkey_time(inrec, "T_STOP", NULL);
00687     cadence=drms_getkey_float(inrec, "T_STEP", NULL);
00688     if (tstart != tstartsave)
00689     {
00690 // to lift this restriction must be able to specify multiple gap files/records as input
00691       sprint_time(tstartstr, tstart, "TAI", 0);
00692       fprintf(stderr, "ERROR: all input records must have same T_START, record %d has %s\n", irec, tstartstr);
00693       return 0;
00694     }
00695 
00696     n_in=(tstop-tstart)/cadence;
00697     if (n_in != gapsize)
00698     {
00699       fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, n_in=%d\n", gapsize, n_in);
00700       return 0;
00701     }
00702 
00703     tstartout=tstart+ifirst*tsample;
00704     tstopout=tstartout+df1;
00705     sprint_time(tstartstr, tstartout, "TAI", 0);
00706 
00707     if (seginflag)
00708       segin = drms_segment_lookup(inrec, segnamein);
00709     else
00710       segin = drms_segment_lookupnum(inrec, 0);
00711     if (segin == NULL)
00712     {
00713       fprintf(stderr, "ERROR: problem looking up input segment: recnum = %lld\n", inrec->recnum);
00714       return 0;
00715     }
00716 
00717     for (is=0; is<6; is++)
00718     {
00719       if (!flagarr[is])
00720         continue;
00721 /*
00722       outrecarr[is] = drms_create_record(drms_env, serieslist[is], DRMS_PERMANENT, &status);
00723       if (status != DRMS_SUCCESS) 
00724       {
00725        fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
00726        return 0;
00727       }
00728 */
00729       outrec=outrecsetarr[is]->records[irec];
00730       if (histrecnumarr[is] > 0)
00731         drms_setlink_static(outrec, histlinkname,  histrecnumarr[is]);
00732       drms_setlink_static(outrec, srclinkname,  inrec->recnum);
00733 
00734       if (segoutflag)
00735         segoutarr[is] = drms_segment_lookup(outrec, segnameout);
00736       else
00737         segoutarr[is] = drms_segment_lookupnum(outrec, 0);
00738       if (segoutarr[is] == NULL)
00739       {
00740         fprintf(stderr, "ERROR: problem looking up outputput segment in series %s\n", serieslist[is]);
00741         return 0;
00742       }
00743 
00744 // the following is not needed if at least one of the first N-1 dimensions are 0 in the jsd, 
00745 // or if the first N-1 dimensions in the jsd are always the dimensions desired.
00746 // it *is* needed any time one must override the jsd defaults
00747 /*
00748       switch(is)
00749       {
00750         case TSOUT:
00751         case FFTOUT:
00752           segoutarr[is]->axis[0]=2*n_sampled_out;
00753           segoutarr[is]->axis[1]=lmode+1;
00754           break;
00755         case FFT1OUT:
00756           segoutarr[is]->axis[0]=2*npow;
00757           segoutarr[is]->axis[1]=2*lmode+1;
00758           break;
00759         case POWOUT:
00760         case ARPOWOUT:
00761           segoutarr[is]->axis[0]=npow;
00762           segoutarr[is]->axis[1]=2*lmode+1;
00763           break;
00764         case MAVGOUT:
00765           segoutarr[is]->axis[0]=npow;
00766           segoutarr[is]->axis[1]=1;
00767           break;
00768         default:
00769           ;
00770       }
00771 */
00772       switch(is)
00773       {
00774         case TSOUT:
00775         case FFTOUT:
00776           tstotal[1]=lmode+1;
00777           break;
00778         case FFT1OUT:
00779           fft1total[1]=2*lmode+1;
00780           break;
00781         case POWOUT:
00782         case ARPOWOUT:
00783           powtotal[1]=2*lmode+1;
00784           break;
00785         case MAVGOUT:
00786         default:
00787           ;
00788       }
00789 
00790     }
00791 
00792     if (mavgout)
00793       for (i=0;i<n_sampled_out/2;i++) 
00794         msum[i] = 0.0;
00795 
00796   /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
00797     for (m=0; m<=lmode; m++) 
00798     {
00799 
00800       if (verbflag > 2) 
00801       {
00802         printf("starting m=%d\n", m);
00803       }
00804 
00805       instart[1]=m;
00806       inend[1]=m;
00807       inarr = drms_segment_readslice(segin, usetype, instart, inend, &status);
00808       if (status != DRMS_SUCCESS || inarr == NULL )
00809       {
00810         fprintf(stderr, "ERROR: problem reading input segment: status = %d, recnum = %lld\n", status, inrec->recnum);
00811         return 0;
00812       }
00813       in_data=(float *)(inarr->data);
00814 
00815     /* Extracts data copy */
00816       extract_data(n_sampled_in, gaps, in_data, data);
00817 
00818     /* Detrend entire time series if required */
00819       if (detrend_const != 0) 
00820       {
00821         detrend_order = floor(detrend_length/detrend_const)+2;
00822         cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps, 
00823              detrend_order, detrend_length, detrend_skip, 
00824              num_sections, sections, detrend_first);
00825       }
00826 
00827     /* Fill gaps in entire time series if required */
00828       memcpy(gaps2, gaps, n_sampled_in*sizeof(int));
00829       if (ifill != 0 && max_ar_order > 0) 
00830       {
00831         cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2,
00832               fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, 
00833               &ar_order, (_Complex float *)ar_coeff);
00834       }
00835 
00836       drms_free_array(inarr);
00837 
00838       memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
00839       if ((ifirst+n_sampled_out) >= n_sampled_in) 
00840         memset(&data[2*(n_sampled_in-ifirst)], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));      
00841 
00842       if (tsout)
00843       {
00844 // code to flip m on this output separately.  in the present code this has already been accomplished by setting flipm.
00845 /*
00846         if (mfliparr[TSOUT] != mflipin)
00847         {
00848           for (i = 0; i < n_sampled_out; i++)
00849           {
00850             datatemp[2*i]=data[2*i];
00851             datatemp[2*i+1]=-data[2*i+1];
00852           }
00853           tmpptr=datatemp;
00854         }
00855         else
00856           tmpptr=data;
00857 */
00858         tsstart[1]=m;
00859         tsend[1]=m;
00860         outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
00861         outarr->bzero=segoutarr[TSOUT]->bzero;
00862         outarr->bscale=segoutarr[TSOUT]->bscale;
00863         status=drms_segment_writeslice_ext(segoutarr[TSOUT], outarr, tsstart, tsend, tstotal, 0);
00864         free(outarr);
00865         if (status != DRMS_SUCCESS)
00866         {
00867           fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[TSOUT], status, m, tstartstr, lmode, histrecnumarr[TSOUT]);
00868           return 0;
00869         }
00870       }
00871 
00872       if (fftout || fft1out || powout || mavgout)
00873         fftwf_execute(plan);
00874 
00875       if (fftout)
00876       {
00877 // code to flip m on this output separately.  in the present code this has already been accomplished by setting flipm.
00878 /*
00879         if (mfliparr[FFTOUT] != mflipin)
00880         {
00881           datatemp[0]=data[0];
00882           datatemp[1]=-data[1];
00883           for (i = 1; i < n_sampled_out; i++)
00884           {
00885             datatemp[2*i]=data[2*(n_sampled_out-i)];
00886             datatemp[2*i+1]=-data[2*(n_sampled_out-i)+1];
00887           }
00888           tmpptr=datatemp;
00889         }
00890         else
00891           tmpptr=data;
00892 */
00893         tsstart[1]=m;
00894         tsend[1]=m;
00895         outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
00896         outarr->bzero=segoutarr[FFTOUT]->bzero;
00897         outarr->bscale=segoutarr[FFTOUT]->bscale;
00898         status=drms_segment_writeslice_ext(segoutarr[FFTOUT], outarr, tsstart, tsend, tstotal, 0);
00899         free(outarr);
00900         if (status != DRMS_SUCCESS)
00901         {
00902           fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[FFTOUT], status, m, tstartstr, lmode, histrecnumarr[FFTOUT]);
00903           return 0;
00904         }
00905       }
00906 
00907       if (fft1out)
00908       {
00909 
00910         fft1start[1]=lmode+m*msignarr[FFT1OUT];
00911         fft1end[1]=lmode+m*msignarr[FFT1OUT];
00912         outarr = drms_array_create(usetype, 2, fft1length, data+2*imin, &status);
00913         outarr->bzero=segoutarr[FFT1OUT]->bzero;
00914         outarr->bscale=segoutarr[FFT1OUT]->bscale;
00915         status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
00916         free(outarr);
00917         if (status != DRMS_SUCCESS)
00918         {
00919           fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[FFT1OUT], status, m, tstartstr, lmode, histrecnumarr[FFT1OUT]);
00920           return 0;
00921         }
00922 
00923         if (m > 0) 
00924         {
00925           /* Do negative m */
00926           for (i=0; i<npow; i++) 
00927           {
00928             /* Conjugate for negative m */
00929             datahold[2*i] = data[2*(n_sampled_out-imin-i)];
00930             datahold[2*i+1] = -data[2*(n_sampled_out-imin-i)+1];
00931           }
00932           fft1start[1]=lmode-m*msignarr[FFT1OUT];
00933           fft1end[1]=lmode-m*msignarr[FFT1OUT];
00934           outarr = drms_array_create(usetype, 2, fft1length, datahold, &status);
00935           outarr->bzero=segoutarr[FFT1OUT]->bzero;
00936           outarr->bscale=segoutarr[FFT1OUT]->bscale;
00937           status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
00938           free(outarr);
00939           if (status != DRMS_SUCCESS)
00940           {
00941             fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[FFT1OUT], status, m, tstartstr, lmode, histrecnumarr[FFT1OUT]);
00942             return 0;
00943           }
00944         }
00945       }
00946 
00947       if (powout || mavgout)
00948         for (i = 0; i < n_sampled_out; i++) 
00949           pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];
00950 
00951       if (powout)
00952       {
00953         powstart[1]=lmode+m*msignarr[POWOUT];
00954         powend[1]=lmode+m*msignarr[POWOUT];
00955         outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status);
00956         outarr->bzero=segoutarr[POWOUT]->bzero;
00957         outarr->bscale=segoutarr[POWOUT]->bscale;
00958         status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
00959         free(outarr);
00960         if (status != DRMS_SUCCESS)
00961         {
00962           fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[POWOUT], status, m, tstartstr, lmode, histrecnumarr[POWOUT]);
00963           return 0;
00964         }
00965 
00966         if (m > 0) 
00967         {
00968           /* Do negative m */
00969           if (imin)
00970             datahold[0] = pow[n_sampled_out-imin];
00971           else
00972             datahold[0] = pow[0];
00973           for (i=1; i<npow;i++) 
00974             datahold[i] = pow[n_sampled_out-imin-i];
00975 
00976           powstart[1]=lmode-m*msignarr[POWOUT];
00977           powend[1]=lmode-m*msignarr[POWOUT];
00978           outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
00979           outarr->bzero=segoutarr[POWOUT]->bzero;
00980           outarr->bscale=segoutarr[POWOUT]->bscale;
00981           status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
00982           free(outarr);
00983           if (status != DRMS_SUCCESS)
00984           {
00985             fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[POWOUT], status, m, tstartstr, lmode, histrecnumarr[POWOUT]);
00986             return 0;
00987           }
00988         }
00989       }
00990 
00991       if (mavgout)
00992       {
00993         mshift = floor(df1*splitting(lmode,m*msignarr[MAVGOUT])+0.5f);
00994         if (mshift >= 0) 
00995           for (i=0; i<n_sampled_out/2-mshift; i++)
00996             msum[i] += pow[mshift+i];
00997         else
00998           for (i=0; i<n_sampled_out/2+mshift; i++)
00999             msum[i-mshift] += pow[i];
01000         if (m > 0) 
01001         {
01002           /* Do negative m */
01003           mshift = floor(df1*splitting(lmode,-m*msignarr[MAVGOUT])+0.5f);
01004           if (mshift >=0)
01005             for (i=1; i<n_sampled_out/2-mshift;i++)
01006               msum[i] += pow[n_sampled_out-mshift-i];
01007           else
01008             for (i=1; i<n_sampled_out/2+mshift; i++)
01009               msum[i-mshift] += pow[n_sampled_out-i];
01010         }
01011       }
01012 
01013       if (arpowout)
01014       {
01015         memmove(data, ar_coeff, 2*(ar_order+1)*sizeof(float));
01016         memset(&data[2*(ar_order+1)],0,2*(n_sampled_out-ar_order-1)*sizeof(float));
01017         fftwf_execute(plan);
01018         for (i = 0; i < n_sampled_out; i++) 
01019           pow[i] = 1.0/(data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]);
01020 
01021         powstart[1]=lmode+m*msignarr[ARPOWOUT];
01022         powend[1]=lmode+m*msignarr[ARPOWOUT];
01023         outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status);
01024         outarr->bzero=segoutarr[ARPOWOUT]->bzero;
01025         outarr->bscale=segoutarr[ARPOWOUT]->bscale;
01026         status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
01027         free(outarr);
01028         if (status != DRMS_SUCCESS)
01029         {
01030           fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[ARPOWOUT], status, m, tstartstr, lmode, histrecnumarr[ARPOWOUT]);
01031           return 0;
01032         }
01033 
01034         if (m > 0) 
01035         {
01036           /* Do negative m */
01037           if (imin)
01038             datahold[0] = pow[n_sampled_out-imin];
01039           else
01040             datahold[0] = pow[0];
01041           for (i=1; i<npow;i++) 
01042             datahold[i] = pow[n_sampled_out-imin-i];
01043 
01044           powstart[1]=lmode-m*msignarr[ARPOWOUT];
01045           powend[1]=lmode-m*msignarr[ARPOWOUT];
01046           outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
01047           outarr->bzero=segoutarr[ARPOWOUT]->bzero;
01048           outarr->bscale=segoutarr[ARPOWOUT]->bscale;
01049           status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
01050           free(outarr);
01051           if (status != DRMS_SUCCESS)
01052           {
01053             fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, m = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[ARPOWOUT], status, m, tstartstr, lmode, histrecnumarr[ARPOWOUT]);
01054              return 0;
01055           }
01056         }
01057       }
01058 
01059     } //end of loop over m
01060 
01061     if (mavgout)
01062     {
01063       c=1.0f/(2*lmode+1);
01064       for (i=0; i<npow; i++) 
01065         datahold[i] = msum[imin+i] * c;
01066       outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
01067       outarr->bzero=segoutarr[MAVGOUT]->bzero;
01068       outarr->bscale=segoutarr[MAVGOUT]->bscale;
01069       status=drms_segment_write(segoutarr[MAVGOUT], outarr, 0);
01070       free(outarr);
01071       if (status != DRMS_SUCCESS)
01072       {
01073         fprintf(stderr, "ERROR: problem writing output segment in series %s: status = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", serieslist[MAVGOUT], status, tstartstr, lmode, histrecnum);
01074         return 0;
01075       }
01076     }
01077 
01078     for (is=0;is<6;is++)
01079     {
01080       if (!flagarr[is])
01081         continue;
01082 
01083       outrec=outrecsetarr[is]->records[irec];
01084       drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
01085       drms_setkey_int(outrec, "QUALITY", 0);
01086       drms_setkey_time(outrec, "T_START", tstartout);
01087       drms_setkey_time(outrec, "T_STOP", tstopout);
01088       drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2);
01089       drms_setkey_float(outrec, "T_STEP", tsamplex);
01090       drms_setkey_int(outrec, "NDT", n_sampled_out);
01091       if (tagflag)
01092         drms_setkey_string(outrec, "TAG", tag);
01093       if (verflag)
01094          drms_setkey_string(outrec, "VERSION", version);
01095 
01096       tnow = (double)time(NULL);
01097       tnow += UNIX_epoch;
01098       drms_setkey_time(outrec, "DATE", tnow);
01099 
01100 //      drms_close_record(outrec, DRMS_INSERT_RECORD);
01101 
01102     }
01103 
01104   }
01105 
01106   drms_close_records(inrecset, DRMS_FREE_RECORD);
01107   for (is=0;is<6;is++)
01108   {
01109     if (!flagarr[is])
01110       continue;
01111     drms_close_records(outrecsetarr[is], DRMS_INSERT_RECORD);
01112   }
01113 
01114   if(ar_coeff != NULL)
01115     free(ar_coeff);
01116 
01117   if (fftout || fft1out || powout || arpowout || mavgout)   
01118     fftwf_destroy_plan(plan); 
01119 
01120   wt=getwalltime();
01121   ct=getcputime(&ut, &st);
01122   if (verbflag) 
01123   {
01124     printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
01125             wt-wt0, ct-ct0);
01126   }
01127 
01128   printf("module %s successful completion\n", cmdparams.argv[0]);
01129   return 0;
01130 }
01131 
01132 
01133 static double splitting(int l, int m)
01134 {
01135   double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0};
01136   double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8};
01137   double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5};
01138   /* f0 is the frequency used for generating the above a-coeff */
01139   double f0 = 1500.0;
01140   /* fcol is the frequency for which to optimize the collaps */
01141   double fcol = 800.0;
01142 
01143   double ll,x,x3,x5,a1,a2,a3,a5,frac,lx;
01144   int ix;
01145 
01146   if (l == 0) 
01147     ll = 1; 
01148   else 
01149     ll = sqrt(l*(l+1.));
01150   x = m/ll;
01151   x3 = x*x*x;
01152   x5 = x3*x*x;
01153   lx = 5*log10(l*f0/fcol)-4;
01154   ix = floor(lx);
01155   frac = lx-ix;
01156   if (lx <= 0) {
01157     ix = 0;
01158     frac = 0.0;
01159   }
01160   if (lx >=8) {
01161     ix = 7;
01162     frac = 1.0;
01163   }
01164   a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1];
01165   a2 = 0.0;
01166   a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1];
01167   a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1];
01168 
01169   return ll*(a1*x+a2*(1.5*x*x-0.5)+a3*(2.5*x3-1.5*x)+a5*(63./8.*x5-70./8.*x3+15./8.*x))/1e9;
01170 }
01171 
01172 
01173 
01174 /* Extract the data samples actually used by skipping or averaging
01175    data. Replace missing data that are not marked as gaps with zero. */
01176 void extract_data(int n_in, int *gaps, float *data_in, float *data_out)
01177 {
01178   int i,j, nmiss = 0;
01179 // the following check will be skipped in production executables since they will be built with NDEBUG defined.
01180 // it doesn't matter because the check is already done in DoIt().
01181   assert(nskip==1 || navg==1);
01182   if ((navg == 1) && (nskip == 1)) 
01183   {
01184     for (i=0; i<n_in; i++) 
01185     {
01186       if (gaps[i]==0 )
01187       {
01188     data_out[2*i] = 0.0f;
01189     data_out[2*i+1] = 0.0f;
01190       }
01191       else
01192       {
01193     if (isnan(data_in[2*i]) || isnan(data_in[2*i+1]))
01194     {
01195       data_out[2*i] = 0.0f;
01196       data_out[2*i+1] = 0.0f;
01197       gaps[i] = 0;
01198       nmiss++;
01199     }
01200     else
01201     {
01202       data_out[2*i] = data_in[2*i];
01203       data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1];
01204     }
01205       }
01206     }
01207   }
01208   else if (nskip != 1) 
01209   {
01210     for (i = 0; i<n_in; i++) 
01211     {
01212       if (gaps[i]==0 )
01213       {
01214     data_out[2*i] = 0.0f;
01215     data_out[2*i+1] = 0.0f;
01216       }
01217       else
01218       {
01219     int ix = nskip*i+noff;
01220     if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
01221     {
01222       data_out[2*i] = 0.0f;
01223       data_out[2*i+1] = 0.0f;
01224       gaps[i] = 0;
01225       nmiss++;
01226     }
01227     else
01228     {
01229       data_out[2*i] = data_in[2*ix];
01230       data_out[2*i+1] =  flipm ? -data_in[2*ix+1] : data_in[2*ix+1];
01231     }
01232       }
01233     }
01234   }
01235   else if (navg != 1) 
01236   {
01237     for (i = 0; i<n_in; i++) 
01238     {
01239       if (gaps[i]==0 )
01240       {
01241     data_out[2*i] = 0.0f;
01242     data_out[2*i+1] = 0.0f;
01243       }
01244       else
01245       {
01246     float avgr = 0.0f;
01247     float avgi = 0.0f;
01248     for (j=0; j<navg; j++) 
01249     {
01250       int ix = navg*i+j;
01251       if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
01252       {
01253         data_out[2*i] = 0.0f;
01254         data_out[2*i+1] = 0.0f;
01255         gaps[i] = 0;
01256         nmiss++;
01257         break;
01258       }
01259       else 
01260       {
01261         avgr += data_in[2*ix];
01262         avgi += data_in[2*ix+1];
01263       }
01264     }
01265     if (j == navg) 
01266     {
01267       data_out[2*i] = avgr/navg;
01268       data_out[2*i+1] = avgi/navg;
01269     }
01270       }
01271     }  
01272   }
01273 }
01274 
01275 /* Extract the windows function for the samples actually used after
01276    subsampling or averaging. Modify the list of continous sections
01277    accordingly, and make sure we do not average across section boundaries. */
01278 static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
01279           float tsample_in, float *tsample_out, 
01280           int *num_sections, int *sections)
01281 {
01282 // the following check will be skipped in production executables since they will be built with NDEBUG defined.
01283 // it doesn't matter because the check is already done in DoIt().
01284   assert(nskip==1 || navg==1);
01285   int i,j,sect, start,stop;
01286 
01287 
01288   if (*num_sections<1)
01289   {
01290     fprintf(stderr,"Apparantly no sections of data available.");
01291     return 1;
01292   }
01293   /* Mask out data that in between sections if this hasn't already been
01294      done in gapfile. */
01295   for (i=0; i<sections[0] && i<n_in; i++)
01296     gaps_in[i] = 0;
01297   for (sect=1; sect<(*num_sections); sect++)
01298   {
01299     for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++)
01300       gaps_in[i] = 0;
01301   }
01302   for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++)
01303     gaps_in[i] = 0;
01304     
01305 
01306   if ((navg == 1) && (nskip == 1)) 
01307   {
01308     *n_out = n_in;
01309     *tsample_out = tsample_in;
01310     for (i=0; i<*n_out; i++) 
01311       gaps_out[i] = gaps_in[i];
01312   }
01313   else if (nskip != 1) 
01314   {
01315     *n_out = n_in/nskip;
01316     *tsample_out = nskip*tsample_in;
01317     for (i=0; i<*n_out; i++) 
01318     {
01319       gaps_out[i] = gaps_in[nskip*i+noff];
01320     }
01321     for (i=0; i<2*(*num_sections); i+=2)
01322     {
01323       start = (int)ceilf(((float)(sections[i]-noff))/nskip);
01324       stop = (int)floorf(((float)(sections[i+1]-noff))/nskip);
01325       if (start <= stop)
01326       {
01327     sections[i] = start;
01328     sections[i+1] = stop;
01329       }
01330       else 
01331       {
01332     /* This section was skipped entirely. */
01333     for (j=i; j< 2*(*num_sections-1); j+=2) 
01334     {
01335       sections[j] = sections[j+2];
01336       sections[j+1] = sections[j+3];
01337     }
01338     *num_sections -= 1;
01339       }
01340     }
01341   }
01342   else  if (navg != 1) 
01343   {
01344     sect = 0;
01345     *n_out = n_in/navg;
01346     *tsample_out = tsample*navg;
01347     for (i=0; i<*n_out; i++) 
01348     {      
01349       int igx = 1;
01350       while (sect < *num_sections && 
01351          !(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1]))
01352     sect++;
01353 
01354       /* Don't allow averaging of data from different sections. */
01355       if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1])
01356       {
01357     for (j=0; j<navg; j++) 
01358       igx *= gaps_in[navg*i+j];
01359     gaps_out[i] = igx;
01360       }
01361       else
01362     gaps_out[i] = 0;
01363     }
01364     for (i=0; i<2*(*num_sections); i+=2)
01365     {
01366       start = (int)ceilf(((float)sections[i])/navg);
01367       stop = (int)floorf(((float)sections[i+1])/navg);
01368       if (start <= stop)
01369       {
01370     sections[i] = start;
01371     sections[i+1] = stop;
01372       }
01373       else 
01374       {
01375     /* This section was skipped entirely. */
01376     for (j=i; j< 2*(*num_sections-1); j+=2) 
01377     {
01378       sections[j] = sections[j+2];
01379       sections[j+1] = sections[j+3];
01380     }
01381     *num_sections -= 1;
01382       }
01383     }
01384   }  
01385   return 0;
01386 }
01387 
01388 
01389 /*
01390    Read a list of sections [start_sample:end_sample] that should
01391    contain continuous data. When detrending, jumps are allow to
01392    occur between sections. The sections file should be a text file 
01393    of the form:
01394 
01395    num
01396    start_1  end_1
01397    start_2  end_2
01398    ...
01399    start_num  end_num
01400 */
01401 static int read_section_file(char *filename, int *num_sections, int **sections)
01402 {
01403   FILE *fh;
01404   int i;
01405 
01406   fh = fopen(filename,"r");
01407   if (fh!=NULL)
01408   {
01409     fscanf(fh,"%d",num_sections);
01410     *sections = (int *)malloc(2*(*num_sections)*sizeof(int));
01411     
01412     for (i=0 ;i<2*(*num_sections); i+=2)
01413     {
01414       fscanf(fh,"%d",&(*sections)[i]);
01415       fscanf(fh,"%d",&(*sections)[i+1]);
01416     
01417       if ((*sections)[i]>(*sections)[i+1] ||
01418       (i>0 && (*sections)[i]<=(*sections)[i-1]))
01419       {
01420     fprintf(stderr,"Invalid sections file, sections overlap.\n");
01421     return 2;
01422       }
01423     }
01424     fclose(fh);
01425     return 0;
01426   }
01427   else 
01428     return 1;
01429 }  
01430