(file) Return to jtsslice.c CVS log (file) (dir) Up to [Development] / JSOC / proj / globalhs / apps

  1 tplarson 1.3 /*
  2              this module creates power spectra out of slices of timeseries
  3              based on jtsfiddle, but no detrending or gapfilling
  4              */
  5 tplarson 1.1 
  6              /* 
  7               *  Detrending/gapfilling/differencing module
  8               *  ts_fiddle_rml upgrades ts_fiddle
  9               */
 10              
 11              #include <stdio.h>
 12              #include <stdlib.h>
 13              #include <strings.h>
 14              #include <time.h>
 15              #include <sys/time.h>
 16              #include <sys/resource.h>
 17              #include <math.h>
 18              #include <assert.h>
 19              #include <fftw3.h>
 20              
 21              #include "jsoc_main.h"
 22 kehcheng 1.7 #include "fitsio.h"
 23 tplarson 1.1 
 24              #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
 25              #define kNOTSPECIFIED "not specified"
 26              
 27 tplarson 1.5 char *module_name = "jtsslice";
 28 tplarson 1.1 
 29              /* Command line arguments: */
 30              ModuleArgs_t module_args[] =
 31              {
 32                {ARG_STRING,   "in", "", "input data records"},
 33                {ARG_STRING,   "out", "", "output data series"},
 34                {ARG_STRING,   "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"},
 35                {ARG_STRING,   "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
 36                {ARG_STRING,   "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
 37                {ARG_STRING,   "srclink",  "SRCDATA", "name of link to source data"},
 38                {ARG_INT,      "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
 39                {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},
 40                {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},
 41              
 42                {ARG_INT,      "NDT", "-1", "", ""},
 43 tplarson 1.3   {ARG_STRING,   "TTOTAL",   NULL,          "total length of time in output"},
 44              
 45 tplarson 1.1   {ARG_INT,      "IFIRST", "0", "", ""},
 46 tplarson 1.3 //  {ARG_INT,      "FLIPM", "0", "", ""},
 47              
 48 tplarson 1.1   {ARG_INT,      "TSOUT", "0", "", ""},
 49                {ARG_INT,      "FFTOUT", "0", "", ""},
 50                {ARG_INT,      "FFT1OUT", "0", "", ""},
 51                {ARG_INT,      "POWOUT", "0", "", ""},
 52                {ARG_INT,      "MAVGOUT", "0", "", ""},
 53                {ARG_INT,      "NAVG", "0", "", ""},
 54                {ARG_INT,      "NSKIP", "0", "", ""},
 55                {ARG_INT,      "NOFF", "0", "", ""},
 56                {ARG_INT,      "IMIN", "0", "", ""},
 57                {ARG_INT,      "IMAX", "-1", "", ""},
 58                {ARG_END},
 59              };
 60              
 61              #include "saveparm.c"
 62              #include "timing.c"
 63 tplarson 1.3 #include "set_history.c"
 64 tplarson 1.1 
 65 tplarson 1.2 #define DIE(code) { fprintf(stderr,"jtsslice died with error code %d\n",(code)); return 0;}
 66 tplarson 1.1 
 67              /* global variables holding the values of the command line variables. */
 68 tplarson 1.3 static char *gapfile;
 69              static int lmode, n_sampled_out, flipm;
 70              static float tsample;
 71              static int ifirst, tsout,fftout, fft1out, powout, mavgout;
 72 tplarson 1.1 static int navg, nskip, noff, imin, imax;
 73              
 74              
 75              /* Prototypes for local functions. */
 76              static double splitting(int l, int m);
 77              
 78              /* Prototypes for external functions */
 79              extern void cffti_(int *, float *);
 80              extern void cfftb_(int *, float *, float *);
 81              static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
 82              		  float tsample_in, float *tsample_out,
 83              		  int *num_sections, int *sections);
 84              static void extract_data(int n_in, int *gaps, float *data_in, float *data_out);
 85              
 86              /************ Here begins the main module **************/
 87              int DoIt(void)
 88              {
 89                int i;
 90                int start_index[2], counts[2], intervals[2];
 91                int m, n_in, n_sampled_in;
 92                TIME epoch, step, start, stop;
 93 tplarson 1.1   int gapsize, istart, istop, data_dim, detrend_order;
 94                int npow, mshift;
 95                float tsamplex, df1, c;
 96                int *agaps;
 97 tplarson 1.3   int *gaps;
 98 tplarson 1.1   float *msum, *pow;
 99                float *data, *out_data, *in_data;
100                fftwf_plan plan;
101                int num_sections, *sections;
102 tplarson 1.2   char tstartstr[100];
103 tplarson 1.1 
104                int fstatus = 0;
105                fitsfile *fitsptr;
106                long fdims[1];
107                long fpixel[]={1};
108                int *anynul=0;
109              
110                int length[2], startind[2], endind[2];
111                float *datarr;
112                int status=DRMS_SUCCESS;
113                int newstat=0;
114 tplarson 1.3   char *ttotal;
115                double nseconds;
116 tplarson 1.1 
117                DRMS_Segment_t *segin = NULL;
118                DRMS_Segment_t *segout = NULL;
119                DRMS_Array_t *inarr = NULL;
120                DRMS_Array_t *outarr = NULL;
121                DRMS_Record_t *inrec = NULL;
122                DRMS_Record_t *outrec = NULL;
123                DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
124                DRMS_RecLifetime_t lifetime;
125                long long histrecnum = -1;
126                DRMS_Segment_t *gapseg = NULL;
127                DRMS_Array_t *gaparr = NULL;
128              
129                HIterator_t outKey_list;
130                DRMS_Keyword_t *outKey;
131                TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
132              
133                int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
134                int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
135              
136 tplarson 1.6   double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout;
137 tplarson 1.1   double wt0, wt1, wt2, wt3, wt;
138                double ut0, ut1, ut2, ut3, ut;
139                double st0, st1, st2, st3, st;
140                double ct0, ct1, ct2, ct3, ct;
141              
142                wt0=getwalltime();
143                ct0=getcputime(&ut0, &st0);
144              
145                /* Read input parameters. */
146 tplarson 1.2   gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
147 tplarson 1.1   n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);
148                ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
149 tplarson 1.3 //  flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);
150 tplarson 1.1   tsout       = cmdparams_save_int    (&cmdparams, "TSOUT", &newstat);
151                fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);
152                fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);
153                powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);
154                mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);
155                navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);
156                nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);
157                noff        = cmdparams_save_int    (&cmdparams, "NOFF", &newstat);
158                imin        = cmdparams_save_int    (&cmdparams, "IMIN", &newstat);
159                imax        = cmdparams_save_int    (&cmdparams, "IMAX", &newstat);
160              
161 tplarson 1.2   char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
162                char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
163                char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
164                char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
165 tplarson 1.1   int seginflag = strcmp(kNOTSPECIFIED, segnamein);
166                int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
167 tplarson 1.2   char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
168                char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
169 tplarson 1.1   int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
170                int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
171                if (permflag)
172                  lifetime = DRMS_PERMANENT;
173                else
174                  lifetime = DRMS_TRANSIENT;
175              
176 tplarson 1.3   ttotal=(char *)cmdparams_save_str(&cmdparams, "TTOTAL", &newstat);
177                status=drms_names_parseduration(&ttotal, &nseconds, 1);
178 tplarson 1.9   if (status != DRMS_SUCCESS)
179 tplarson 1.2   {
180 tplarson 1.3     fprintf(stderr, "ERROR: problem parsing TTOTAL, = %s\n", ttotal);
181                  return 1; 
182 tplarson 1.2   }
183              
184                if (newstat) 
185                {
186                  fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
187                  cpsave_decode_error(newstat);
188                  return 1;
189                }  
190                else if (savestrlen != strlen(savestr)) 
191                {
192                  fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
193                  return 1;
194                }
195              
196 tplarson 1.1 
197                /**** Sanity check of input parameters. ****/
198              
199                /* Default is to just output the timeseries. */
200 tplarson 1.3   if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0)) 
201 tplarson 1.1     tsout=1;
202                /* Only one type of output allowed at a time */
203 tplarson 1.3   if ((tsout+fftout+fft1out+powout+mavgout) !=1) 
204 tplarson 1.1   {
205                  fprintf(stderr, "ERROR: only one type of output allowed at a time\n");
206 tplarson 1.2     return 1;
207 tplarson 1.1   }
208                if (navg <= 0) 
209                  navg=1;
210                if (nskip <= 0) 
211                  nskip=1;
212                if ((navg != 1) && (nskip != 1))
213                {
214                  fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");
215 tplarson 1.2     return 1;
216 tplarson 1.1   }
217                if (noff < 0 || noff >= nskip) 
218                {
219                  fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
220 tplarson 1.2     return 1;
221 tplarson 1.1   }
222 tplarson 1.3 
223              
224                DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
225              					       outseries,
226              					       DRMS_TRANSIENT, 
227              					       &status);
228              
229                if (status != DRMS_SUCCESS) 
230 tplarson 1.1   {
231 tplarson 1.3     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
232 tplarson 1.2     return 1;
233 tplarson 1.1   }
234              
235 tplarson 1.3   DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
236                DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);
237              
238 tplarson 1.11   char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsslice.c,v 1.10 2012/09/20 17:33:45 tplarson Exp $");
239 tplarson 1.3    if (histlink != NULL) 
240                 {
241 tplarson 1.5      histrecnum=set_history(histlink, cvsinfo);
242 tplarson 1.3      if (histrecnum < 0)
243                   {
244                     drms_close_record(tempoutrec, DRMS_FREE_RECORD);
245                     return 1;
246                   }
247                 }
248                 else
249                 {
250                   fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
251                 }
252 tplarson 1.1  
253                // these must be present in the output dataseries and variable, not links or constants   
254 tplarson 1.5    char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
255 tplarson 1.1  
256                 int itest;
257                 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
258                 {
259                   DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
260 tplarson 1.9      if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
261 tplarson 1.1      {
262                     fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
263                     drms_close_record(tempoutrec, DRMS_FREE_RECORD);
264 tplarson 1.2        return 1;
265 tplarson 1.1      }
266                 }
267               
268 tplarson 1.6    cadence=drms_getkey_float(tempoutrec, "T_STEP", &status);
269                 double chunksecs = n_sampled_out*cadence;
270 tplarson 1.3    if (fmod(nseconds,chunksecs) != 0.0)
271                 {
272                   fprintf(stderr, "ERROR: input parameters seem incompatible (chunksecs must divide nseconds): nseconds = %f, chunksecs = %f\n", nseconds, chunksecs);
273                   drms_close_record(tempoutrec, DRMS_FREE_RECORD);
274                   return 1;
275                 }
276                 int ntimechunks = nseconds/chunksecs;
277 tplarson 1.6    tsample=cadence;
278 tplarson 1.3  
279                 int mflipout = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
280                 if (status != DRMS_SUCCESS) mflipout=0;
281               
282 tplarson 1.1    drms_close_record(tempoutrec, DRMS_FREE_RECORD);
283               
284               
285                 if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
286                 {
287                   DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status);
288                   if (status != DRMS_SUCCESS || gaprecset == NULL) 
289                   {
290                     fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
291 tplarson 1.2        return 1;
292 tplarson 1.1      }
293                   if (gaprecset->n == 0)
294                   {
295                     fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
296 tplarson 1.2        return 1;
297 tplarson 1.1      }
298                   if (gaprecset->n > 1)
299                   {
300                     fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
301 tplarson 1.2        return 1;
302 tplarson 1.1      } 
303                   gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);
304                   if (gapseg != NULL)
305                     gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status);
306                   if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
307                   {
308                     fprintf(stderr, "ERROR: problem reading gaps segment: status = %d\n", status);
309 tplarson 1.2        return 1;
310 tplarson 1.1      }
311                   agaps=(int *)(gaparr->data);
312                   gapsize=gaparr->axis[0];
313                 }
314                 else
315                 {
316                   fits_get_img_size(fitsptr, 1, fdims, &fstatus);
317                   gapsize=fdims[0];
318                   agaps=(int *)(malloc(gapsize*sizeof(int)));
319                   fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus); 
320                   fits_close_file(fitsptr, &fstatus);
321                   if (fstatus) 
322                   {
323                     fprintf(stderr, "ERROR: ");
324                     fits_report_error(stderr, fstatus);
325 tplarson 1.2        return 1;
326 tplarson 1.1      }
327                 }
328               
329 tplarson 1.9    if (verbflag)
330                   printf("gapfile read, gapsize = %d\n",gapsize);
331 tplarson 1.1  
332                 data_dim=gapsize;
333                 if (n_sampled_out>gapsize) 
334                   data_dim=n_sampled_out;
335               
336 tplarson 1.3    num_sections=1;
337                 sections = malloc(2*sizeof(int));
338                 sections[0] = 0;
339                 sections[1] = data_dim-1;
340 tplarson 1.1  
341                 /* allocate temporary storage */
342                 gaps=(int *)(malloc(gapsize*sizeof(int)));
343                 data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
344               
345 tplarson 1.2    DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);
346 tplarson 1.1  
347 tplarson 1.2    if (status != DRMS_SUCCESS || inrecset == NULL) 
348 tplarson 1.1    {
349                   fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
350 tplarson 1.2      return 1;
351 tplarson 1.1    }
352               
353 tplarson 1.2    if (inrecset->n == 0)
354 tplarson 1.1    {
355 tplarson 1.4      fprintf(stderr, "ERROR: input recordset contains no records\n");
356 tplarson 1.2      return 1;
357 tplarson 1.1    }
358               
359 tplarson 1.2    inrec=inrecset->records[0];
360 tplarson 1.1    char *inchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "T_STEP"};
361               
362                 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
363                 {
364                   DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
365 tplarson 1.9      if (inkeytest == NULL)
366 tplarson 1.1      {
367                     fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
368 tplarson 1.2        drms_close_records(inrecset, DRMS_FREE_RECORD);
369                     return 1;
370 tplarson 1.1      }
371                 }
372               
373 tplarson 1.6    tstartsave=drms_getkey_time(inrec, "T_START", NULL); //must be constant across all input records unless GAPFILE is implemented differently
374               
375 tplarson 1.3    int mflipin = drms_getkey_int(inrec, "MFLIPPED", &status);
376                 if (status != DRMS_SUCCESS) mflipin=0;
377               
378                 if (mflipout != mflipin)
379                   flipm=1;
380                 else
381                   flipm=0;
382               
383 tplarson 1.1    // changed n_in to gapsize in following call
384                 if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
385                 {
386                   fprintf(stderr, "ERROR: problem in extract_gaps\n");
387                   return 0;
388                 }
389               
390                 if (n_sampled_out == -1) 
391                   n_sampled_out = n_sampled_in;
392                 df1 = tsamplex*n_sampled_out;
393               
394 tplarson 1.3    if (fftout || fft1out || powout || mavgout)     
395 tplarson 1.1      plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, 
396               			    (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);
397               
398                 /* Set sum to zero before starting */
399                 if (mavgout) 
400                   msum = (float *)calloc(n_sampled_out/2,sizeof(float));
401               
402 tplarson 1.5    if (fft1out || powout || mavgout)
403 tplarson 1.1    {
404                   /* Set first and last output index if not set as a parameter. */
405                   if (imax < imin) 
406                   {
407                     imin=0;
408                     imax=n_sampled_out/2-1;
409                   }
410                   npow=imax-imin+1;
411                   pow=(float *)(malloc(n_sampled_out*sizeof(float)));
412                 }
413               
414               
415 tplarson 1.5    status=drms_stage_records(inrecset, 1, 0);
416                 if (status != DRMS_SUCCESS)
417                 {
418                   fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status);
419                   return 1;
420                 }
421 tplarson 1.1  
422 tplarson 1.5    int nrecsout = ntimechunks * inrecset->n;
423                 DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecsout, outseries, DRMS_PERMANENT, &status);
424                 if (status != DRMS_SUCCESS || outrecset == NULL) 
425 tplarson 1.1    {
426 tplarson 1.5     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
427                  return 1;
428 tplarson 1.1    }
429               
430 tplarson 1.9    int gapoffset, gapoffset0 = ifirst;
431 tplarson 1.8    double tstart0 = tstartsave+ifirst*tsample;
432 tplarson 1.6    int irec;
433 tplarson 1.5    for (irec=0;irec<inrecset->n;irec++)
434 tplarson 1.3    {
435               
436 tplarson 1.5      inrec = inrecset->records[irec];
437 tplarson 1.1  
438 tplarson 1.5      int lmin=drms_getkey_int(inrec, "LMIN", NULL);
439                   int lmax=drms_getkey_int(inrec, "LMAX", NULL);
440                   if (lmin != lmax)
441 tplarson 1.3      {
442 tplarson 1.5        fprintf(stderr, "ERROR: lmin != lmax not yet supported.\n");
443 tplarson 1.3        return 0;
444                   }
445 tplarson 1.5      lmode=lmin;
446 tplarson 1.3  
447 tplarson 1.9      if (verbflag)
448                     printf("starting l=%d\n",lmode);
449               
450 tplarson 1.6      tstart=drms_getkey_time(inrec, "T_START", NULL);
451                   tstop =drms_getkey_time(inrec, "T_STOP", NULL);
452                   cadence=drms_getkey_float(inrec, "T_STEP", NULL);
453                   if (tstart != tstartsave)
454                   {
455               // to lift this restriction must be able to specify multiple gap files/records as input
456                     sprint_time(tstartstr, tstart, "TAI", 0);
457                     fprintf(stderr, "ERROR: all input records must have same T_START, record %d has %s\n", irec, tstartstr);
458                     return 0;
459                   }
460 tplarson 1.1  
461 tplarson 1.6      n_in=(tstop-tstart)/cadence;
462 tplarson 1.5      if (n_in != gapsize)
463 tplarson 1.3      {
464 tplarson 1.5        fprintf(stderr, "ERROR: gaps seem incompatible with time-series, gapsize=%d, n_in=%d\n", gapsize, n_in);
465                     return 0;
466 tplarson 1.3      }
467 tplarson 1.1  
468 tplarson 1.5      /* Create output data set. */
469                   if (tsout || fftout) 
470                   {
471                     length[0]=2*n_sampled_out;
472                     length[1]=lmode+1;
473                   }
474                   else 
475                   {
476                     if (fft1out) 
477                       length[0]=2*npow;
478                     else 
479                       length[0]=npow;
480                     if (powout || fft1out) 
481                       length[1]=2*lmode+1;
482                     else 
483                       length[1]=1;
484                   }
485 tplarson 1.1  
486 tplarson 1.5      startind[1]=0;
487                   endind[1]=lmode;
488               
489                   if (seginflag)
490                     segin = drms_segment_lookup(inrec, segnamein);
491 tplarson 1.3      else
492 tplarson 1.5        segin = drms_segment_lookupnum(inrec, 0);
493                   if (segin == NULL)
494 tplarson 1.3      {
495 tplarson 1.5        fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld\n", inrec->recnum);
496                     return 0;
497 tplarson 1.3      }
498 tplarson 1.5  
499 tplarson 1.6      int iset;
500 tplarson 1.5      for (iset=0; iset < ntimechunks; iset++)
501                   {
502                     tstartout=tstart0 + iset * chunksecs;
503                     tstopout=tstartout+chunksecs;
504                     sprint_time(tstartstr, tstartout, "TAI", 0);
505 tplarson 1.9        gapoffset = gapoffset0 + iset*n_sampled_out;
506               
507                     if (verbflag>1)
508                     {
509                       printf("  starting tstart = %s, ",tstartstr);
510 tplarson 1.11         if (irec == 0)
511                       {
512                         int ig, gtotal=0;
513                         for (ig=0;ig<n_sampled_out;ig++)
514                           gtotal+=gaps[gapoffset+ig];
515                         float dutycycle = (float)gtotal/n_sampled_out;
516                         printf("dutycycle = %f\n", dutycycle);
517                       }
518 tplarson 1.9        }
519 tplarson 1.5  
520 tplarson 1.9  // set ifirst to gapoffset for reading slice
521 tplarson 1.6        ifirst=gapoffset;
522 tplarson 1.9        startind[0]=2*(ifirst);
523                     endind[0]=2*(ifirst + n_sampled_out) - 1;
524 tplarson 1.6  // set ifirst=0 because data is read starting at ifirst by drms_segment_readslice
525                     ifirst=0;
526 tplarson 1.5  
527                     inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);
528                     if (status != DRMS_SUCCESS || inarr == NULL)
529                     {
530                       fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld, status = %d\n", inrec->recnum, status);
531                       return 0;
532                     }
533                     datarr=(float *)(inarr->data);
534               
535               /*
536                     outrec = drms_create_record(drms_env, outseries, lifetime, &status);
537                     if (status != DRMS_SUCCESS || outrec==NULL)
538                     {
539                        fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
540                        return 0;
541                     }
542               */
543 tplarson 1.6  //      outrec = outrecset->records[ntimechunks*irec+iset];
544                     outrec = outrecset->records[iset*inrecset->n + irec];
545 tplarson 1.5  
546 tplarson 1.9        if (histlink != NULL)
547 tplarson 1.5          drms_setlink_static(outrec, histlinkname,  histrecnum);
548 tplarson 1.9        if (srclink != NULL)
549 tplarson 1.5          drms_setlink_static(outrec, srclinkname,  inrec->recnum);
550               
551                     if (segoutflag)
552                       segout = drms_segment_lookup(outrec, segnameout);
553                     else
554                       segout = drms_segment_lookupnum(outrec, 0);
555                     outarr = drms_array_create(usetype, 2, length, NULL, &status);
556                     if (status != DRMS_SUCCESS || outarr == NULL || segout == NULL)
557                     {
558                       fprintf(stderr,"ERROR: problem creating output array or segment: length = [%d, %d], status = %d", length[0], length[1], status);
559                       return 0; 
560                     }
561                     out_data = outarr->data;
562 tplarson 1.1  
563                 /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
564 tplarson 1.5        for (m=0; m<=lmode; m++) 
565                     {
566 tplarson 1.1  
567 tplarson 1.9          if (verbflag>2)
568                         printf("    starting m=%d\n",m);
569               
570 tplarson 1.5          in_data=datarr + m*2*n_sampled_out;
571 tplarson 1.1  
572 tplarson 1.5          /* Extracts data copy */
573                       extract_data(n_sampled_out, gaps+gapoffset, in_data, data);
574 tplarson 1.1  
575 tplarson 1.5          /* pad with zeros if the last point output (n_sampled_out+ifirst) 
576                          is after the last data points read in (nread). */
577                       memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
578                       if ((ifirst+n_sampled_out) >= n_sampled_in) 
579                         memset(&data[2*n_sampled_in], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));      
580 tplarson 1.1  
581 tplarson 1.5          if (fftout || fft1out || powout || mavgout) 
582                         fftwf_execute(plan);
583 tplarson 1.1  
584 tplarson 1.5          if (tsout || fftout) 
585                         memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));
586 tplarson 1.1  
587 tplarson 1.5          else if (fft1out) 
588 tplarson 1.3          {
589 tplarson 1.5          /* Do positive m */
590                         memcpy(&out_data[2*(lmode+m)*npow], &data[2*imin], 2*npow*sizeof(float));
591                         if (m > 0) 
592 tplarson 1.3            {
593 tplarson 1.5            /* Do negative m */
594               	    for (i=0; i<npow; i++) 
595               	    {
596                          /* Conjugate for negative m */
597               	      out_data[2*((lmode-m)*npow+i)] = data[2*(n_sampled_out-imin-i)];
598               	      out_data[2*((lmode-m)*npow+i)+1] = -data[2*(n_sampled_out-imin-i)+1];
599               	    }
600 tplarson 1.3            }
601                       }
602                       else 
603                       {
604 tplarson 1.5          /* Compute power spectrum from complex FFT. */      
605                         for (i = 0; i < n_sampled_out; i++) 
606                           pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];
607                         if (powout) 
608                         {
609                         /* Do positive m */
610                           memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));
611               	    if (m > 0) 
612                           {
613                           /* Do negative m */
614               	      if (imin)
615               	        out_data[(lmode-m)*npow] = pow[n_sampled_out-imin];
616               	      else
617               	        out_data[(lmode-m)*npow] = pow[0];
618               	      for (i=1; i<npow;i++) 
619               	        out_data[(lmode-m)*npow+i] = pow[n_sampled_out-imin-i];
620                           }
621                         }
622                         else 
623 tplarson 1.3            {
624 tplarson 1.5            /* m-averaged output */
625                         /* Sum all freqs, not only those in output. Should be fixed. */
626                         /* Maybe should allow for wrapping around Nyquist. */
627                         /* Do positive m */
628                           mshift = floor(df1*splitting(lmode,m)+0.5f);
629                           if (mshift >= 0) 
630                             for (i=0; i<n_sampled_out/2-mshift; i++)
631                               msum[i] += pow[mshift+i];
632 tplarson 1.3              else
633 tplarson 1.5                for (i=0; i<n_sampled_out/2+mshift; i++)
634                               msum[i-mshift] += pow[i];
635                           if (m > 0) 
636                           {
637                           /* Do negative m */
638                             mshift = floor(df1*splitting(lmode,-m)+0.5f);
639                             if (mshift >=0)
640                               for (i=1; i<n_sampled_out/2-mshift;i++)
641                                 msum[i] += pow[n_sampled_out-mshift-i];
642                             else
643                               for (i=1; i<n_sampled_out/2+mshift; i++)
644                                 msum[i-mshift] += pow[n_sampled_out-i];
645                           }
646 tplarson 1.3            }
647                       }
648               
649 tplarson 1.5        } /* End of m loop */
650               
651                     drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
652                     drms_setkey_int(outrec, "QUALITY", 0);
653                     drms_setkey_time(outrec, "T_START", tstartout);
654                     drms_setkey_time(outrec, "T_STOP", tstopout);
655                     drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2);
656                     drms_setkey_float(outrec, "T_STEP", tsamplex);
657                     drms_setkey_int(outrec, "NDT", n_sampled_out);
658               
659                     tnow = (double)time(NULL);
660                     tnow += UNIX_epoch;
661                     drms_setkey_time(outrec, "DATE", tnow);
662 tplarson 1.3  
663 tplarson 1.5        if (mavgout) 
664                     {
665                       c=1.0f/(2*lmode+1);
666                       for (i=0; i<npow; i++) 
667                         out_data[i] = msum[imin+i] * c;
668                     }
669 tplarson 1.3  
670 tplarson 1.9        outarr->bzero=segout->bzero;
671                     outarr->bscale=segout->bscale;
672 tplarson 1.5        status=drms_segment_write(segout, outarr, 0);
673                     if (status != DRMS_SUCCESS)
674                     {
675                       fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", status, tstartstr, lmode, histrecnum);
676                       return 0;
677                     }
678 tplarson 1.8        drms_free_array(outarr);
679                     drms_free_array(inarr);
680 tplarson 1.5  //    drms_close_record(outrec, DRMS_INSERT_RECORD);
681 tplarson 1.1  
682 tplarson 1.3      }
683 tplarson 1.1  
684 tplarson 1.3    }
685 tplarson 1.1  
686 tplarson 1.5    drms_close_records(inrecset,  DRMS_FREE_RECORD);
687                 drms_close_records(outrecset, DRMS_INSERT_RECORD);
688 tplarson 1.2  
689 tplarson 1.3    if (fftout || fft1out || powout || mavgout)   
690 tplarson 1.1      fftwf_destroy_plan(plan); 
691               
692                 wt=getwalltime();
693                 ct=getcputime(&ut, &st);
694                 if (verbflag) 
695                 {
696                   fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
697                           wt-wt0, ct-ct0);
698                 }
699               
700 tplarson 1.3    printf("module %s successful completion\n", cmdparams.argv[0]);
701 tplarson 1.1    return 0;
702               }
703               
704               
705               static double splitting(int l, int m)
706               {
707                 double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0};
708                 double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8};
709                 double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5};
710                 /* f0 is the frequency used for generating the above a-coeff */
711                 double f0 = 1500.0;
712                 /* fcol is the frequency for which to optimize the collaps */
713                 double fcol = 800.0;
714               
715                 double ll,x,x3,x5,a1,a2,a3,a5,frac,lx;
716                 int ix;
717               
718                 if (l == 0) 
719                   ll = 1; 
720                 else 
721                   ll = sqrt(l*(l+1.));
722 tplarson 1.1    x = m/ll;
723                 x3 = x*x*x;
724                 x5 = x3*x*x;
725                 lx = 5*log10(l*f0/fcol)-4;
726                 ix = floor(lx);
727                 frac = lx-ix;
728                 if (lx <= 0) {
729                   ix = 0;
730                   frac = 0.0;
731                 }
732                 if (lx >=8) {
733                   ix = 7;
734                   frac = 1.0;
735                 }
736                 a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1];
737                 a2 = 0.0;
738                 a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1];
739                 a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1];
740               
741                 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;
742               }
743 tplarson 1.1  
744               
745               
746               /* Extract the data samples actually used by skipping or averaging
747                  data. Replace missing data that are not marked as gaps with zero. */
748               void extract_data(int n_in, int *gaps, float *data_in, float *data_out)
749               {
750                 int i,j, nmiss = 0;
751 tplarson 1.2  // the following check will be skipped in production executables since they will be built with NDEBUG defined.
752               // it doesn't matter because the check is already done in DoIt().
753 tplarson 1.1    assert(nskip==1 || navg==1);
754                 if ((navg == 1) && (nskip == 1)) 
755                 {
756                   for (i=0; i<n_in; i++) 
757                   {
758                     if (gaps[i]==0 )
759                     {
760               	data_out[2*i] = 0.0f;
761               	data_out[2*i+1] = 0.0f;
762                     }
763                     else
764                     {
765               	if (isnan(data_in[2*i]) || isnan(data_in[2*i+1]))
766               	{
767               	  data_out[2*i] = 0.0f;
768               	  data_out[2*i+1] = 0.0f;
769               	  gaps[i] = 0;
770               	  nmiss++;
771               	}
772               	else
773               	{
774 tplarson 1.1  	  data_out[2*i] = data_in[2*i];
775               	  data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1];
776               	}
777                     }
778                   }
779                 }
780                 else if (nskip != 1) 
781                 {
782                   for (i = 0; i<n_in; i++) 
783                   {
784                     if (gaps[i]==0 )
785                     {
786               	data_out[2*i] = 0.0f;
787               	data_out[2*i+1] = 0.0f;
788                     }
789                     else
790                     {
791               	int ix = nskip*i+noff;
792               	if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
793               	{
794               	  data_out[2*i] = 0.0f;
795 tplarson 1.1  	  data_out[2*i+1] = 0.0f;
796               	  gaps[i] = 0;
797               	  nmiss++;
798               	}
799               	else
800               	{
801               	  data_out[2*i] = data_in[2*ix];
802               	  data_out[2*i+1] =  flipm ? -data_in[2*ix+1] : data_in[2*ix+1];
803               	}
804                     }
805                   }
806                 }
807                 else if (navg != 1) 
808                 {
809                   for (i = 0; i<n_in; i++) 
810                   {
811                     if (gaps[i]==0 )
812                     {
813               	data_out[2*i] = 0.0f;
814               	data_out[2*i+1] = 0.0f;
815                     }
816 tplarson 1.1        else
817                     {
818               	float avgr = 0.0f;
819               	float avgi = 0.0f;
820               	for (j=0; j<navg; j++) 
821               	{
822               	  int ix = navg*i+j;
823               	  if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
824               	  {
825               	    data_out[2*i] = 0.0f;
826               	    data_out[2*i+1] = 0.0f;
827               	    gaps[i] = 0;
828               	    nmiss++;
829               	    break;
830               	  }
831               	  else 
832               	  {
833               	    avgr += data_in[2*ix];
834               	    avgi += data_in[2*ix+1];
835               	  }
836               	}
837 tplarson 1.1  	if (j == navg) 
838               	{
839               	  data_out[2*i] = avgr/navg;
840               	  data_out[2*i+1] = avgi/navg;
841               	}
842                     }
843                   }  
844                 }
845               }
846               
847               /* Extract the windows function for the samples actually used after
848                  subsampling or averaging. Modify the list of continous sections
849                  accordingly, and make sure we do not average across section boundaries. */
850               static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
851               		  float tsample_in, float *tsample_out, 
852               		  int *num_sections, int *sections)
853               {
854 tplarson 1.2  // the following check will be skipped in production executables since they will be built with NDEBUG defined.
855               // it doesn't matter because the check is already done in DoIt().
856 tplarson 1.1    assert(nskip==1 || navg==1);
857                 int i,j,sect, start,stop;
858               
859               
860                 if (*num_sections<1)
861                 {
862                   fprintf(stderr,"Apparantly no sections of data available.");
863                   return 1;
864                 }
865                 /* Mask out data that in between sections if this hasn't already been
866                    done in gapfile. */
867                 for (i=0; i<sections[0] && i<n_in; i++)
868                   gaps_in[i] = 0;
869                 for (sect=1; sect<(*num_sections); sect++)
870                 {
871                   for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++)
872                     gaps_in[i] = 0;
873                 }
874                 for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++)
875                   gaps_in[i] = 0;
876                   
877 tplarson 1.1  
878                 if ((navg == 1) && (nskip == 1)) 
879                 {
880                   *n_out = n_in;
881                   *tsample_out = tsample_in;
882                   for (i=0; i<*n_out; i++) 
883                     gaps_out[i] = gaps_in[i];
884                 }
885                 else if (nskip != 1) 
886                 {
887                   *n_out = n_in/nskip;
888                   *tsample_out = nskip*tsample_in;
889                   for (i=0; i<*n_out; i++) 
890                   {
891                     gaps_out[i] = gaps_in[nskip*i+noff];
892                   }
893                   for (i=0; i<2*(*num_sections); i+=2)
894                   {
895                     start = (int)ceilf(((float)(sections[i]-noff))/nskip);
896                     stop = (int)floorf(((float)(sections[i+1]-noff))/nskip);
897                     if (start <= stop)
898 tplarson 1.1        {
899               	sections[i] = start;
900               	sections[i+1] = stop;
901                     }
902                     else 
903                     {
904               	/* This section was skipped entirely. */
905               	for (j=i; j< 2*(*num_sections-1); j+=2) 
906               	{
907               	  sections[j] = sections[j+2];
908               	  sections[j+1] = sections[j+3];
909               	}
910               	*num_sections -= 1;
911                     }
912                   }
913                 }
914                 else  if (navg != 1) 
915                 {
916                   sect = 0;
917                   *n_out = n_in/navg;
918                   *tsample_out = tsample*navg;
919 tplarson 1.1      for (i=0; i<*n_out; i++) 
920                   {      
921                     int igx = 1;
922                     while (sect < *num_sections && 
923               	     !(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1]))
924               	sect++;
925               
926                     /* Don't allow averaging of data from different sections. */
927                     if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1])
928                     {
929               	for (j=0; j<navg; j++) 
930               	  igx *= gaps_in[navg*i+j];
931               	gaps_out[i] = igx;
932                     }
933                     else
934               	gaps_out[i] = 0;
935                   }
936                   for (i=0; i<2*(*num_sections); i+=2)
937                   {
938                     start = (int)ceilf(((float)sections[i])/navg);
939                     stop = (int)floorf(((float)sections[i+1])/navg);
940 tplarson 1.1        if (start <= stop)
941                     {
942               	sections[i] = start;
943               	sections[i+1] = stop;
944                     }
945                     else 
946                     {
947               	/* This section was skipped entirely. */
948               	for (j=i; j< 2*(*num_sections-1); j+=2) 
949               	{
950               	  sections[j] = sections[j+2];
951               	  sections[j+1] = sections[j+3];
952               	}
953               	*num_sections -= 1;
954                     }
955                   }
956                 }  
957                 return 0;
958               }
959               

Karen Tian
Powered by
ViewCVS 0.9.4