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

Diff for /JSOC/proj/globalhs/apps/jtsslice.c between version 1.2 and 1.3

version 1.2, 2011/08/16 02:29:10 version 1.3, 2011/11/15 05:08:46
Line 1 
Line 1 
 /* this JSOC module is a port of an SOI module written by Rasmus Larsen. (Danish spelling)  /*
  * ported by Tim Larson. (Swedish spelling)  this module creates power spectra out of slices of timeseries
  * this version probably no longer works with nonzero NAVG or NSKIP  based on jtsfiddle, but no detrending or gapfilling
  * for those options use jtsfiddle instead of jtsslice  
  */  */
  
 /* /*
Line 17 
Line 16 
 #include <sys/resource.h> #include <sys/resource.h>
 #include <math.h> #include <math.h>
 #include <assert.h> #include <assert.h>
 //#include <complex.h>  
 #include <fftw3.h> #include <fftw3.h>
 #include "fahlman_ulrych_C99.h" #include "fahlman_ulrych_C99.h"
 #include "detrend_C99.h" #include "detrend_C99.h"
 //#include <xmem.h>  
  
 #include "jsoc_main.h" #include "jsoc_main.h"
 #include "/home/jsoc/include/fitsio.h" #include "/home/jsoc/include/fitsio.h"
Line 29 
Line 26 
 #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0])) #define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
 #define kNOTSPECIFIED "not specified" #define kNOTSPECIFIED "not specified"
  
 char *module_name = "jtsslice";  char *module_name = "jtsslice2";
  
 /* Command line arguments: */ /* Command line arguments: */
 ModuleArgs_t module_args[] = ModuleArgs_t module_args[] =
Line 42  ModuleArgs_t module_args[] =
Line 39  ModuleArgs_t module_args[] =
   {ARG_STRING,   "srclink",  "SRCDATA", "name of link to source data"},   {ARG_STRING,   "srclink",  "SRCDATA", "name of link to source data"},
   {ARG_INT,      "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},   {ARG_INT,      "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
   {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},   {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},
   {ARG_STRING,   "LOGFILE", "stdout", ""},  
   {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},   {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},
   {ARG_STRING,   "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"},  
  
   {ARG_INT,      "NDT", "-1", "", ""},   {ARG_INT,      "NDT", "-1", "", ""},
 //  {ARG_FLOAT,    "TSAMPLE", "60", "", ""},    {ARG_STRING,   "TTOTAL",   NULL,          "total length of time in output"},
   {ARG_INT,      "IFILL", "0", "", ""}, //changed from jtsfiddle  
   {ARG_INT,      "IFIRST", "0", "", ""},   {ARG_INT,      "IFIRST", "0", "", ""},
   {ARG_INT,      "FLIPM", "0", "", ""},  //  {ARG_INT,      "FLIPM", "0", "", ""},
   {ARG_INT,      "MAX_AR_ORDER", "360", "", ""},  
   {ARG_INT,      "FU_FIT_EDGES", "1", "", ""},  
   {ARG_INT,      "FU_ITER", "3", "", ""},  
   {ARG_INT,      "FU_MIN_PERCENT", "90", "", ""},  
   {ARG_FLOAT,    "CDETREND", "0.0", "", ""}, //changed from jtsfiddle  
   {ARG_INT,      "DETREND_LENGTH", "1600", "", ""},  
   {ARG_INT,      "DETREND_SKIP", "1440", "", ""},  
   {ARG_INT,      "DETREND_FIRST", "0", "", ""},  
   {ARG_INT,      "TSOUT", "0", "", ""},   {ARG_INT,      "TSOUT", "0", "", ""},
   {ARG_INT,      "FFTOUT", "0", "", ""},   {ARG_INT,      "FFTOUT", "0", "", ""},
   {ARG_INT,      "FFT1OUT", "0", "", ""},   {ARG_INT,      "FFT1OUT", "0", "", ""},
   {ARG_INT,      "POWOUT", "0", "", ""},   {ARG_INT,      "POWOUT", "0", "", ""},
   {ARG_INT,      "ARPOWOUT", "0", "", ""},  
   {ARG_INT,      "MAVGOUT", "0", "", ""},   {ARG_INT,      "MAVGOUT", "0", "", ""},
   {ARG_INT,      "NAVG", "0", "", ""},   {ARG_INT,      "NAVG", "0", "", ""},
   {ARG_INT,      "NSKIP", "0", "", ""},   {ARG_INT,      "NSKIP", "0", "", ""},
Line 75  ModuleArgs_t module_args[] =
Line 62  ModuleArgs_t module_args[] =
  
 #include "saveparm.c" #include "saveparm.c"
 #include "timing.c" #include "timing.c"
   #include "set_history.c"
  
 #define DIE(code) { fprintf(stderr,"jtsslice died with error code %d\n",(code)); return 0;} #define DIE(code) { fprintf(stderr,"jtsslice died with error code %d\n",(code)); return 0;}
  
 /* global variables holding the values of the command line variables. */ /* global variables holding the values of the command line variables. */
 static char *logfile, *gapfile, *sectionfile;  static char *gapfile;
 static int lmode, n_sampled_out, ifill, flipm;  static int lmode, n_sampled_out, flipm;
 static int ar_order, max_ar_order, fu_iter, fu_fit_edges;  static float tsample;
 static int fu_min_percent;  static int ifirst, tsout,fftout, fft1out, powout, mavgout;
 static int detrend_length, detrend_skip, detrend_first;  
 static float tsample, detrend_const;  
 static int ifirst, tsout,fftout, fft1out, powout, arpowout, mavgout;  
 static int navg, nskip, noff, imin, imax; static int navg, nskip, noff, imin, imax;
  
  
Line 99  static int extract_gaps(int n_in, int *g
Line 84  static int extract_gaps(int n_in, int *g
                   float tsample_in, float *tsample_out,                   float tsample_in, float *tsample_out,
                   int *num_sections, int *sections);                   int *num_sections, int *sections);
 static void extract_data(int n_in, int *gaps, float *data_in, float *data_out); static void extract_data(int n_in, int *gaps, float *data_in, float *data_out);
 static int read_section_file(char *filename, int *num_section, int **section);  
  
 /************ Here begins the main module **************/ /************ Here begins the main module **************/
 int DoIt(void) int DoIt(void)
Line 112  int DoIt(void)
Line 96  int DoIt(void)
   int npow, mshift;   int npow, mshift;
   float tsamplex, df1, c;   float tsamplex, df1, c;
   int *agaps;   int *agaps;
   int *gaps, *gaps2;    int *gaps;
   float *msum, *pow;   float *msum, *pow;
   float *data, *out_data, *in_data;   float *data, *out_data, *in_data;
   fftwf_plan plan;   fftwf_plan plan;
   int num_sections, *sections;   int num_sections, *sections;
   float *ar_coeff=NULL;  
   char tstartstr[100];   char tstartstr[100];
  
   int fstatus = 0;   int fstatus = 0;
Line 130  int DoIt(void)
Line 113  int DoIt(void)
   float *datarr;   float *datarr;
   int status=DRMS_SUCCESS;   int status=DRMS_SUCCESS;
   int newstat=0;   int newstat=0;
     char *ttotal;
     double nseconds;
  
   DRMS_Segment_t *segin = NULL;   DRMS_Segment_t *segin = NULL;
   DRMS_Segment_t *segout = NULL;   DRMS_Segment_t *segout = NULL;
Line 160  int DoIt(void)
Line 145  int DoIt(void)
   ct0=getcputime(&ut0, &st0);   ct0=getcputime(&ut0, &st0);
  
   /* Read input parameters. */   /* Read input parameters. */
   logfile     = (char *)cmdparams_save_str    (&cmdparams, "LOGFILE", &newstat);  
   gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);   gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
   sectionfile = (char *)cmdparams_save_str    (&cmdparams, "SECTIONFILE", &newstat);  
   n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);   n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);
 //  tsample     = cmdparams_save_float  (&cmdparams, "TSAMPLE", &newstat);  
   ifill       = cmdparams_save_int    (&cmdparams, "IFILL", &newstat);  
   max_ar_order= cmdparams_save_int    (&cmdparams, "MAX_AR_ORDER", &newstat);  
   fu_fit_edges= cmdparams_save_int    (&cmdparams, "FU_FIT_EDGES", &newstat);  
   fu_iter     = cmdparams_save_int    (&cmdparams, "FU_ITER", &newstat);  
   fu_min_percent = cmdparams_save_int (&cmdparams, "FU_MIN_PERCENT", &newstat);  
   ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);   ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
   flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);  //  flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);
   detrend_const = cmdparams_save_float(&cmdparams, "CDETREND", &newstat);  
   detrend_length = cmdparams_save_int (&cmdparams, "DETREND_LENGTH", &newstat);  
   detrend_skip = cmdparams_save_int   (&cmdparams, "DETREND_SKIP", &newstat);  
   detrend_first = cmdparams_save_int  (&cmdparams, "DETREND_FIRST", &newstat);  
   tsout       = cmdparams_save_int    (&cmdparams, "TSOUT", &newstat);   tsout       = cmdparams_save_int    (&cmdparams, "TSOUT", &newstat);
   fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);   fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);
   fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);   fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);
   powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);   powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);
   arpowout    = cmdparams_save_int    (&cmdparams, "ARPOWOUT", &newstat);  
   mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);   mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);
   navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);   navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);
   nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);   nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);
Line 203  int DoIt(void)
Line 175  int DoIt(void)
   else   else
     lifetime = DRMS_TRANSIENT;     lifetime = DRMS_TRANSIENT;
  
 //#include "histinclude.c"    ttotal=(char *)cmdparams_save_str(&cmdparams, "TTOTAL", &newstat);
     status=drms_names_parseduration(&ttotal, &nseconds, 1);
   DRMS_Record_t *tempoutrec = drms_create_record(drms_env,    if (status)
                                                outseries,  
                                                DRMS_TRANSIENT,  
                                                &status);  
   
   if (status != DRMS_SUCCESS)  
   {   {
     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);      fprintf(stderr, "ERROR: problem parsing TTOTAL, = %s\n", ttotal);
     return 1;     return 1;
   }   }
  
   DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);  
   DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);  
   
   if (histlink != NULL)  
   {  
     DRMS_Record_t *histrec = drms_create_record(drms_env,  
                                                   histlink->info->target_series,  
                                                   DRMS_PERMANENT,  
                                                   &status);  
     if (status != DRMS_SUCCESS)  
     {  
       fprintf(stderr,"ERROR: could not open a record in history dataseries %s, status = %d\n", histlink->info->target_series, status);  
       return 1;  
     }  
   
     histrecnum = histrec->recnum;  
     tnow = (double)time(NULL);  
     tnow += UNIX_epoch;  
     status = drms_setkey_time(histrec, "DATE", tnow);  
     if (status != DRMS_SUCCESS)  
     {  
       fprintf(stderr,"ERROR: problem writing keyword DATE in history dataseries, status = %d\n", status);  
       drms_close_record(histrec, DRMS_FREE_RECORD);  
       return 1;  
     }  
     status = drms_setkey_string(histrec, "MODNAME", cmdparams.argv[0]);  
     if (status != DRMS_SUCCESS)  
     {  
       fprintf(stderr,"ERROR: problem writing keyword MODNAME in history dataseries, status = %d\n", status);  
       drms_close_record(histrec, DRMS_FREE_RECORD);  
       return 1;  
     }  
     status = drms_setkey_string(histrec, "ARGSUSED", savestr);  
     if (status != DRMS_SUCCESS)  
     {  
       fprintf(stderr,"ERROR: problem writing keyword DATE in history dataseries, status = %d\n", status);  
       drms_close_record(histrec, DRMS_FREE_RECORD);  
       return 1;  
     }  
     drms_close_record(histrec, DRMS_INSERT_RECORD);  
   
   }  
   else  
   {  
     fprintf(stderr,"WARNING: could not find history link in output dataseries\n");  
   }  
   
   if (newstat)   if (newstat)
   {   {
     fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);     fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
Line 279  int DoIt(void)
Line 199  int DoIt(void)
   /**** Sanity check of input parameters. ****/   /**** Sanity check of input parameters. ****/
  
   /* Default is to just output the timeseries. */   /* Default is to just output the timeseries. */
   if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0))    if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0))
     tsout=1;     tsout=1;
   /* Only one type of output allowed at a time */   /* Only one type of output allowed at a time */
   if ((tsout+fftout+fft1out+powout+mavgout+arpowout) !=1)    if ((tsout+fftout+fft1out+powout+mavgout) !=1)
   {   {
     fprintf(stderr, "ERROR: only one type of output allowed at a time\n");     fprintf(stderr, "ERROR: only one type of output allowed at a time\n");
     return 1;     return 1;
Line 301  int DoIt(void)
Line 221  int DoIt(void)
     fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");     fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
     return 1;     return 1;
   }   }
   if(arpowout && !ifill)  
   
     DRMS_Record_t *tempoutrec = drms_create_record(drms_env,
                                                  outseries,
                                                  DRMS_TRANSIENT,
                                                  &status);
   
     if (status != DRMS_SUCCESS)
   {   {
     fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");      fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
     return 1;     return 1;
   }   }
  
     DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
     DRMS_Link_t *srclink = hcon_lookup_lower(&tempoutrec->links, srclinkname);
   
     if (histlink != NULL)
     {
       histrecnum=set_history(histlink);
       if (histrecnum < 0)
       {
         drms_close_record(tempoutrec, DRMS_FREE_RECORD);
         return 1;
       }
     }
     else
     {
       fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
     }
  
  // these must be present in the output dataseries and variable, not links or constants  // these must be present in the output dataseries and variable, not links or constants
   char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX"};   char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX"};
Line 323  int DoIt(void)
Line 266  int DoIt(void)
     }     }
   }   }
  
     tstep=drms_getkey_float(tempoutrec, "T_STEP", &status);
     double chunksecs = n_sampled_out*tstep;
     if (fmod(nseconds,chunksecs) != 0.0)
     {
       fprintf(stderr, "ERROR: input parameters seem incompatible (chunksecs must divide nseconds): nseconds = %f, chunksecs = %f\n", nseconds, chunksecs);
   drms_close_record(tempoutrec, DRMS_FREE_RECORD);   drms_close_record(tempoutrec, DRMS_FREE_RECORD);
       return 1;
     }
     int ntimechunks = nseconds/chunksecs;
  
 //  xmem_off();    int mflipout = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
     if (status != DRMS_SUCCESS) mflipout=0;
   
     drms_close_record(tempoutrec, DRMS_FREE_RECORD);
  
  
   if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))   if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
Line 378  int DoIt(void)
Line 332  int DoIt(void)
   if (n_sampled_out>gapsize)   if (n_sampled_out>gapsize)
     data_dim=n_sampled_out;     data_dim=n_sampled_out;
  
   /* Read location of jumps. */  
   if (!strcmp(sectionfile,kNOTSPECIFIED))  
   {  
     /* No sections file specified. Assume that the whole  
        time series in one section. */  
     num_sections=1;     num_sections=1;
     sections = malloc(2*sizeof(int));     sections = malloc(2*sizeof(int));
     sections[0] = 0;     sections[0] = 0;
     sections[1] = data_dim-1;     sections[1] = data_dim-1;
   }  
   else  
   {  
     int secstat;  
     if (secstat=read_section_file(sectionfile, &num_sections, &sections))  
     {  
       DRMS_RecordSet_t *secrecset = drms_open_records(drms_env, sectionfile, &status);  
       if (status != DRMS_SUCCESS || secrecset == NULL)  
       {  
         fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status);  
         return 1;  
       }  
       if (secrecset->n == 0)  
       {  
         fprintf(stderr,"ERROR: sectionfile recordset contains no records\n");  
         return 1;  
       }  
       if (secrecset->n > 1)  
       {  
         fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n");  
         return 1;  
       }  
       num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL);  
       if (num_sections < 1)  
       {  
         fprintf(stderr,"ERROR: Invalid NSECS keywords\n");  
         return 1;  
       }  
       sections = (int *)malloc(2*(num_sections)*sizeof(int));  
       char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL);  
       sscanf(strtok(sectkey," \n"),"%d",&(sections[0]));  
       sscanf(strtok(NULL," \n"),"%d",&(sections[1]));  
       int i;  
       for (i=2 ;i<2*(num_sections); i+=2)  
       {  
         sscanf(strtok(NULL," \n"), "%d", &(sections)[i]);  
         sscanf(strtok(NULL," \n"), "%d", &(sections)[i+1]);  
   
         if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1]))  
         {  
           fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n");  
           return 1;  
         }  
       }  
       free(sectkey);  
     }  
   }  
   
   printf("num_sections = %d\n",num_sections);  
  
   /* allocate temporary storage */   /* allocate temporary storage */
   gaps=(int *)(malloc(gapsize*sizeof(int)));   gaps=(int *)(malloc(gapsize*sizeof(int)));
   gaps2=(int *)(malloc(gapsize*sizeof(int)));  
   data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));   data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
   ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float));  
   
   
  
   DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);   DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);
  
Line 455  int DoIt(void)
Line 351  int DoIt(void)
  
   if (inrecset->n == 0)   if (inrecset->n == 0)
   {   {
     printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]);      printf("ERROR: input recordset contains no records\n");
     return 1;     return 1;
   }   }
  
   if (inrecset->n > 1)   if (inrecset->n > 1)
   {   {
     printf("ERROR: nRecs > 1 not yet supported.\n");      printf("ERROR: nrecs > 1 not yet supported.\n");
     return 1;     return 1;
   }   }
  
   
   inrec=inrecset->records[0];   inrec=inrecset->records[0];
   char *inchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "T_STEP"};   char *inchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "T_STEP"};
  
Line 480  int DoIt(void)
Line 375  int DoIt(void)
     }     }
   }   }
  
     int mflipin = drms_getkey_int(inrec, "MFLIPPED", &status);
     if (status != DRMS_SUCCESS) mflipin=0;
   
     if (mflipout != mflipin)
       flipm=1;
     else
       flipm=0;
   
   status=drms_stage_records(inrecset, 1, 0);   status=drms_stage_records(inrecset, 1, 0);
   if (status != DRMS_SUCCESS)   if (status != DRMS_SUCCESS)
   {   {
Line 500  int irec=0;
Line 403  int irec=0;
   }   }
   lmode=lmin;   lmode=lmin;
  
   
   tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL);   tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL);
   tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL);   tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL);
   tstep=drms_getkey_float(inrec, "T_STEP", NULL);  //  tstep=drms_getkey_float(inrec, "T_STEP", NULL);
  
   n_in=(tstop-tstart)/tstep;   n_in=(tstop-tstart)/tstep;
   if (n_in != gapsize)   if (n_in != gapsize)
Line 514  int irec=0;
Line 416  int irec=0;
  
   tsample=tstep;   tsample=tstep;
  
   
   
   // changed n_in to gapsize in following call   // changed n_in to gapsize in following call
   if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))   if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
   {   {
Line 523  int irec=0;
Line 423  int irec=0;
     return 0;     return 0;
   }   }
  
   /* Adjust detrend parameters according to the number of  
      available points. */  
   if (n_sampled_in<detrend_length)  
   {  
     detrend_length = n_sampled_in;  
     detrend_skip = detrend_length;  
   }  
   
   int idtf;  
   if (detrend_first > 0)  
   {  
     for (idtf=0;idtf<detrend_first;idtf++)  
       gaps[idtf]=0;  
   }  
   
   if (n_sampled_out == -1)   if (n_sampled_out == -1)
     n_sampled_out = n_sampled_in;     n_sampled_out = n_sampled_in;
   df1 = tsamplex*n_sampled_out;   df1 = tsamplex*n_sampled_out;
  
   if (fftout || fft1out || powout || arpowout || mavgout)    if (fftout || fft1out || powout || mavgout)
     plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data,     plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data,
                             (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);                             (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);
  
Line 550  int irec=0;
Line 435  int irec=0;
   if (mavgout)   if (mavgout)
     msum = (float *)calloc(n_sampled_out/2,sizeof(float));     msum = (float *)calloc(n_sampled_out/2,sizeof(float));
  
   
   tstartout=tstart+ifirst*tsample;  
   tstopout=tstartout+df1;  
   sprint_time(tstartstr, tstartout, "TAI", 0);  
   
   
   /* Create output data set. */   /* Create output data set. */
   if (tsout || fftout)   if (tsout || fftout)
   {   {
Line 576  int irec=0;
Line 455  int irec=0;
       length[0]=2*npow;       length[0]=2*npow;
     else     else
       length[0]=npow;       length[0]=npow;
     if (powout || arpowout || fft1out)      if (powout || fft1out)
       length[1]=2*lmode+1;       length[1]=2*lmode+1;
     else     else
       length[1]=1;       length[1]=1;
Line 590  int irec=0;
Line 469  int irec=0;
      were nan's not removed in the input window function.  the present version      were nan's not removed in the input window function.  the present version
      assumes the window function is correct. */      assumes the window function is correct. */
  
   startind[0]=2*ifirst;  
   endind[0]=2*(ifirst + n_sampled_out) - 1;  
   startind[1]=0;   startind[1]=0;
   endind[1]=lmode;   endind[1]=lmode;
  
Line 599  int irec=0;
Line 476  int irec=0;
     segin = drms_segment_lookup(inrec, segnamein);     segin = drms_segment_lookup(inrec, segnamein);
   else   else
     segin = drms_segment_lookupnum(inrec, 0);     segin = drms_segment_lookupnum(inrec, 0);
   if (segin)    if (segin == NULL)
 //    inarr = drms_segment_read(segin, usetype, &status);    {
       fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld\n", inrec->recnum);
       return 0;
     }
   
     double tstart0 = tstart+ifirst*tsample;
     int iset=0;
     for (iset=0; iset < ntimechunks; iset++)
     {
       tstartout=tstart0 + iset * chunksecs;
       tstopout=tstartout+chunksecs;
       sprint_time(tstartstr, tstartout, "TAI", 0);
   
       startind[0]=2*(ifirst+iset*n_sampled_out);
       endind[0]=2*(ifirst+iset*n_sampled_out + n_sampled_out) - 1;
   
     inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);     inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);
   if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL)      if (status != DRMS_SUCCESS || inarr == NULL)
   {   {
     fprintf(stderr, "ERROR: problem reading input segment: irec = %d, status = %d\n", irec, status);        fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld, status = %d\n", inrec->recnum, status);
     return 0;     return 0;
   }   }
   datarr=(float *)(inarr->data);   datarr=(float *)(inarr->data);
  
   
   outrec = drms_create_record(drms_env,   outrec = drms_create_record(drms_env,
                               outseries,                               outseries,
                               lifetime,                               lifetime,
Line 645  int irec=0;
Line 536  int irec=0;
   /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/   /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
   for (m=0; m<=lmode; m++)   for (m=0; m<=lmode; m++)
   {   {
 //    (*history)("### processing m = %d ###\n",m);  
     /* Read chunk of data corresponding to next m from the data file. */  
 /*  
     start_index[1] = m;  
     if (!(tile_sds = vds_select(in_vds,0,0,1,1,start_index,counts,intervals)))  
     {  
       fprintf(stderr,"vds_select failed");  
       DIE(9);  
     }  
     if (!sds_data(tile_sds))  
     {  
       fprintf(stderr,"*** lnu vds_select returned no data.\n");  
       sds_free(&tile_sds);  
       DIE(10);  
     }  
     in_data = (float *)(sds_data(tile_sds));  
     if (sds_dim0(tile_sds)<0)  
     {  
       sds_free(&tile_sds);  
       DIE(soi_errno);  
     }  
 */  
  
 //    in_data=datarr + m*2*gapsize;  
     in_data=datarr + m*2*n_sampled_out;     in_data=datarr + m*2*n_sampled_out;
 // changed this expression after using drms_segment_readslice above  
   
  
     /* Extracts data copy */     /* Extracts data copy */
 //    extract_data(n_sampled_in, gaps, in_data, data);  
    extract_data(n_sampled_out, gaps+gapoffset, in_data, data);    extract_data(n_sampled_out, gaps+gapoffset, in_data, data);
 // likewise  
   
     /* Detrend entire time series if required */  
     if (detrend_const != 0)  
     {  
       detrend_order = floor(detrend_length/detrend_const)+2;  
   
 //      for (i=0;i<num_sections; i++)  
 //        printf("[%d:%d]\n",sections[2*i],sections[2*i+1]);  
       cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps+gapoffset,  
                          detrend_order, detrend_length, detrend_skip,  
                          num_sections, sections, detrend_first);  
     }  
   
     /* Fill gaps in entire time series if required */  
     memcpy(gaps2, gaps, n_sampled_in*sizeof(int));  
     if (ifill != 0 && max_ar_order > 0)  
     {  
 //      (*history)("Filling gaps, ");  
       cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2+gapoffset,  
                       fu_min_percent, max_ar_order, fu_iter, fu_fit_edges,  
                       &ar_order, (_Complex float *)ar_coeff);  
 /*  
       {  
         FILE *fh;  
         unsigned char c;  
         fh = fopen("gaps.bin","w");  
         for(i=0; i<n_sampled_in; i++)  
         {  
           c = (unsigned char) gaps2[i];  
           fwrite(&c,1,1,fh);  
         }  
         fclose(fh);  
       }  
 */  
 //      (*history)("AR model order used = %d.\n",ar_order);  
     }  
   
  
     /* pad with zeros if the last point output (n_sampled_out+ifirst)     /* pad with zeros if the last point output (n_sampled_out+ifirst)
        is after the last data points read in (nread). */        is after the last data points read in (nread). */
     if (arpowout)  
     {  
       memmove(data, ar_coeff, 2*(ar_order+1)*sizeof(float));  
       memset(&data[2*(ar_order+1)],0,2*(n_sampled_out-ar_order-1)*sizeof(float));  
     }  
     else  
     {  
       memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));       memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
       if ((ifirst+n_sampled_out) >= n_sampled_in)       if ((ifirst+n_sampled_out) >= n_sampled_in)
         memset(&data[2*n_sampled_in], 0,          memset(&data[2*n_sampled_in], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));
                2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));  
     }  
     if (fftout || fft1out || powout || arpowout || mavgout)  
     {  
 //      (*history)("Computing FFT.\n");  
       fftwf_execute(plan);  
     }  
   
  
 //    (*history)("Preparing output.\n");        if (fftout || fft1out || powout || mavgout)
           fftwf_execute(plan);
  
     if (tsout || fftout)     if (tsout || fftout)
       memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));       memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));
Line 759  int irec=0;
Line 572  int irec=0;
     else     else
     {     {
       /* Compute power spectrum from complex FFT. */       /* Compute power spectrum from complex FFT. */
       if (arpowout)  
       {  
         for (i = 0; i < n_sampled_out; i++)  
           pow[i] = 1.0/(data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]);  
       }  
       else  
       {  
         for (i = 0; i < n_sampled_out; i++)         for (i = 0; i < n_sampled_out; i++)
           pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];           pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];
       }          if (powout)
       if (powout || arpowout)  
       {       {
         /* Do positive m */         /* Do positive m */
         memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));         memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));
Line 797  int irec=0;
Line 602  int irec=0;
         else         else
           for (i=0; i<n_sampled_out/2+mshift; i++)           for (i=0; i<n_sampled_out/2+mshift; i++)
             msum[i-mshift] += pow[i];             msum[i-mshift] += pow[i];
         if (m > 0) {            if (m > 0)
             {
           /* Do negative m */           /* Do negative m */
           mshift = floor(df1*splitting(lmode,-m)+0.5f);           mshift = floor(df1*splitting(lmode,-m)+0.5f);
           if (mshift >=0)           if (mshift >=0)
Line 812  int irec=0;
Line 618  int irec=0;
  
   } /* End of m loop */   } /* End of m loop */
  
 /*  
   hiter_new(&outKey_list, &outrec->keywords);  
   while ( (outKey=(DRMS_Keyword_t *)hiter_getnext(&outKey_list)) )  
      {  
      char *keyName = outKey->info->name;  
      DRMS_Value_t inValue = {DRMS_TYPE_STRING, NULL};  
      inValue = drms_getkey_p(inrec, keyName, &status);  
      if (status) fprintf(stderr,"*** COPY drms_getkey_p %s status=%d\n",keyName,status);  
      drms_setkey_p(outrec, keyName, &inValue);  
      if ((inValue.type == DRMS_TYPE_STRING) && inValue.value.string_val)  
      free(inValue.value.string_val);  
      inValue.value.string_val = NULL;  
      }  
 */  
   drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);   drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
   drms_setkey_int(outrec, "QUALITY", 0);   drms_setkey_int(outrec, "QUALITY", 0);
   drms_setkey_time(outrec, "T_START", tstartout);   drms_setkey_time(outrec, "T_START", tstartout);
Line 834  int irec=0;
Line 626  int irec=0;
   drms_setkey_float(outrec, "T_STEP", tsamplex);   drms_setkey_float(outrec, "T_STEP", tsamplex);
   drms_setkey_int(outrec, "NDT", n_sampled_out);   drms_setkey_int(outrec, "NDT", n_sampled_out);
  
   
   
   tnow = (double)time(NULL);   tnow = (double)time(NULL);
   tnow += UNIX_epoch;   tnow += UNIX_epoch;
   drms_setkey_time(outrec, "DATE", tnow);   drms_setkey_time(outrec, "DATE", tnow);
  
   drms_free_array(inarr);  
   if (inrecset)  
   {  
     drms_close_records(inrecset, DRMS_FREE_RECORD);  
   }  
   
   if (mavgout)   if (mavgout)
   {   {
     c=1.0f/(2*lmode+1);     c=1.0f/(2*lmode+1);
Line 853  int irec=0;
Line 637  int irec=0;
       out_data[i] = msum[imin+i] * c;       out_data[i] = msum[imin+i] * c;
   }   }
  
   if(ar_coeff)  
     free(ar_coeff);  
   
   status=drms_segment_write(segout, outarr, 0);   status=drms_segment_write(segout, outarr, 0);
   if (status != DRMS_SUCCESS)   if (status != DRMS_SUCCESS)
   {   {
Line 865  int irec=0;
Line 646  int irec=0;
  
   drms_close_record(outrec, DRMS_INSERT_RECORD);   drms_close_record(outrec, DRMS_INSERT_RECORD);
  
   if (fftout || fft1out || powout || arpowout || mavgout)    }
   
     drms_free_array(inarr);
     if (inrecset)
     {
       drms_close_records(inrecset, DRMS_FREE_RECORD);
     }
   
   
     if (fftout || fft1out || powout || mavgout)
     fftwf_destroy_plan(plan);     fftwf_destroy_plan(plan);
  
   wt=getwalltime();   wt=getwalltime();
Line 875  int irec=0;
Line 665  int irec=0;
 //    printf("number of records created = %d\n", nsuccess); //    printf("number of records created = %d\n", nsuccess);
     fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n",     fprintf(stdout, "total time spent: %.2f ms wall time, %.2f ms cpu time\n",
             wt-wt0, ct-ct0);             wt-wt0, ct-ct0);
     printf("module %s successful completion\n", cmdparams.argv[0]);  
   }   }
  
 //  xmem_leakreport();    printf("module %s successful completion\n", cmdparams.argv[0]);
   return 0;   return 0;
 } }
  
Line 1138  static int extract_gaps(int n_in, int *g
Line 927  static int extract_gaps(int n_in, int *g
   return 0;   return 0;
 } }
  
   
 /*  
    Read a list of sections [start_sample:end_sample] that should  
    contain continuous data. When detrending, jumps are allow to  
    occur between sections. The sections file should be a text file  
    of the form:  
   
    num  
    start_1  end_1  
    start_2  end_2  
    ...  
    start_num  end_num  
 */  
 static int read_section_file(char *filename, int *num_sections, int **sections)  
 {  
   FILE *fh;  
   int i;  
   
   fh = fopen(filename,"r");  
   if (fh!=NULL)  
   {  
     fscanf(fh,"%d",num_sections);  
     *sections = (int *)malloc(2*(*num_sections)*sizeof(int));  
   
     for (i=0 ;i<2*(*num_sections); i+=2)  
     {  
       fscanf(fh,"%d",&(*sections)[i]);  
       fscanf(fh,"%d",&(*sections)[i+1]);  
   
       if ((*sections)[i]>(*sections)[i+1] ||  
           (i>0 && (*sections)[i]<=(*sections)[i-1]))  
       {  
         fprintf(stderr,"Invalid sections file, sections overlap.\n");  
         return 2;  
       }  
     }  
     fclose(fh);  
     return 0;  
   }  
   else  
     return 1;  
 }  
   


Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

Karen Tian
Powered by
ViewCVS 0.9.4