/* this JSOC module is a port of an SOI module written by Rasmus Larsen.
 * ported by Tim Larson.
 * tim doubts this works correctly with nskip or navg
 */

/* 
 *  Detrending/gapfilling/differencing module
 *  ts_fiddle_rml upgrades ts_fiddle
 */

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <math.h>
#include <assert.h>
#include <fftw3.h>

#include "jsoc_main.h"
#include "fitsio.h"

#define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
#define kNOTSPECIFIED "not specified"

char *module_name = "jtsfiddle";
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 $";

/* Command line arguments: */
ModuleArgs_t module_args[] =
{
  {ARG_STRING,   "in", "", "input data records"},
  {ARG_STRING,   "tsout", kNOTSPECIFIED, "output dataseries for timeseries"},
  {ARG_STRING,   "powout", kNOTSPECIFIED, "output dataseries for power spectra"},
  {ARG_STRING,   "fftout", kNOTSPECIFIED, "output dataseries for fft's"},
  {ARG_STRING,   "fft1out", kNOTSPECIFIED, "output dataseries for fft's reordered"},
  {ARG_STRING,   "arpowout", kNOTSPECIFIED, "output dataseries for AR power"},
  {ARG_STRING,   "mavgout", kNOTSPECIFIED, "output dataseries for m-averaged power spectra"},
//  {ARG_STRING,   "out", "", "output data series"},
  {ARG_STRING,   "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"},
  {ARG_STRING,   "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"},
  {ARG_STRING,   "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
  {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,      "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,   "TAG", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, usually used as an additional primekey"},
  {ARG_STRING,   "VERSION", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, useful for selecting obsolete records"},

  {ARG_STRING,   "LOGFILE", "stdout", ""},
  {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_FLOAT,    "TSAMPLE", "60", "", ""},
  {ARG_INT,      "IFILL", "1", "", ""},
  {ARG_INT,      "IFIRST", "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", "50.0", "", ""},
  {ARG_INT,      "DETREND_LENGTH", "1600", "", ""},
  {ARG_INT,      "DETREND_SKIP", "1440", "", ""},
  {ARG_INT,      "DETREND_FIRST", "0", "", ""},
/*
  {ARG_INT,      "TSOUT", "0", "", ""},
  {ARG_INT,      "FFTOUT", "0", "", ""},
  {ARG_INT,      "FFT1OUT", "0", "", ""},
  {ARG_INT,      "POWOUT", "0", "", ""},
  {ARG_INT,      "ARPOWOUT", "0", "", ""},
  {ARG_INT,      "MAVGOUT", "0", "", ""},
*/
  {ARG_INT,      "NAVG", "0", "", ""},
  {ARG_INT,      "NSKIP", "0", "", ""},
  {ARG_INT,      "NOFF", "0", "", ""},
  {ARG_INT,      "IMIN", "0", "", ""},
  {ARG_INT,      "IMAX", "-1", "", ""},
  {ARG_END},
};

#include "saveparm.c"
#include "timing.c"
#include "set_history.c"

extern void cdetrend_discontig( int n, _Complex float *data, int *isgood, 
				int degree, int length, int skip, 
				int m, int *sect_last, int detrend_first);

int cfahlman_ulrych(int n, _Complex float *data, int *isgood, 
		    int minpercentage, int maxorder, int iterations, 
		    int padends, int *order, _Complex float *ar_coeff);

//char *getdetrendversion(void);
//char *getgapfillversion(void);

/* global variables holding the values of the command line variables. */
static char *logfile, *gapfile, *sectionfile;
static int lmode, n_sampled_out, ifill, flipm;
static int ar_order, max_ar_order, fu_iter, fu_fit_edges;
static int fu_min_percent;
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;


/* Prototypes for local functions. */
static double splitting(int l, int m);

/* Prototypes for external functions */
extern void cffti_(int *, float *);
extern void cfftb_(int *, float *, float *);
static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
		  float tsample_in, float *tsample_out,
		  int *num_sections, int *sections);
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 **************/
int DoIt(void)
{
  int i;
  int m, n_in, n_sampled_in;
  TIME epoch, step, start, stop;
  int gapsize, istart, istop, data_dim, detrend_order;
  int npow, mshift;
  float tsamplex, df1, c;
  int *agaps;
  int *gaps, *gaps2;
  float *msum, *pow;
  float *data, *out_data, *in_data;
  fftwf_plan plan;
  int num_sections, *sections;
  float *ar_coeff=NULL;
  float *datarr, *datptr, *outptr;
  float *datahold;
  char tstartstr[100];
  int lmin, lmax;

  int fstatus = 0;
  fitsfile *fitsptr;
  long fdims[1];
  long fpixel[]={1};
  int *anynul=0;

  int instart[2], inend[2];
  int tslength[2], tsstart[2], tsend[2], tstotal[2];
  int fft1length[2], fft1start[2], fft1end[2], fft1total[2];
  int powlength[2], powstart[2], powend[2], powtotal[2];
  int status=DRMS_SUCCESS;
  int newstat=0;

  DRMS_Segment_t *segin = NULL;
  DRMS_Segment_t *segout = NULL;
  DRMS_Array_t *inarr = NULL;
  DRMS_Array_t *outarr = NULL;
  DRMS_Record_t *inrec = NULL;
  DRMS_Record_t *outrec = NULL;
  DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
  DRMS_RecLifetime_t lifetime;
  long long histrecnum = -1;
  DRMS_Segment_t *gapseg = NULL;
  DRMS_Array_t *gaparr = NULL;

  HIterator_t outKey_list;
  DRMS_Keyword_t *outKey;
  TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */

  int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
  int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);

  double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout;
  double wt0, wt1, wt2, wt3, wt;
  double ut0, ut1, ut2, ut3, ut;
  double st0, st1, st2, st3, st;
  double ct0, ct1, ct2, ct3, ct;

  wt0=getwalltime();
  ct0=getcputime(&ut0, &st0);

  /* Read input parameters. */
  logfile     = (char *)cmdparams_save_str    (&cmdparams, "LOGFILE", &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);
//  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);
//  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);
  fftout      = cmdparams_save_int    (&cmdparams, "FFTOUT", &newstat);
  fft1out     = cmdparams_save_int    (&cmdparams, "FFT1OUT", &newstat);
  powout      = cmdparams_save_int    (&cmdparams, "POWOUT", &newstat);
  arpowout    = cmdparams_save_int    (&cmdparams, "ARPOWOUT", &newstat);
  mavgout     = cmdparams_save_int    (&cmdparams, "MAVGOUT", &newstat);
*/
  navg        = cmdparams_save_int    (&cmdparams, "NAVG", &newstat);
  nskip       = cmdparams_save_int    (&cmdparams, "NSKIP", &newstat);
  noff        = cmdparams_save_int    (&cmdparams, "NOFF", &newstat);
  imin        = cmdparams_save_int    (&cmdparams, "IMIN", &newstat);
  imax        = cmdparams_save_int    (&cmdparams, "IMAX", &newstat);

  char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);
//  char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
  char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);
  char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);
  int seginflag = strcmp(kNOTSPECIFIED, segnamein);
  int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
  char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
  char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
  int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
  int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
  if (permflag)
    lifetime = DRMS_PERMANENT;
  else
    lifetime = DRMS_TRANSIENT;

  char *tag = (char *)cmdparams_save_str(&cmdparams, "TAG", &newstat);
  int tagflag = strcmp(kNOTSPECIFIED, tag);
  char *version = (char *)cmdparams_save_str(&cmdparams, "VERSION", &newstat);
  int verflag = strcmp(kNOTSPECIFIED, version);

  char *tsseries = (char *)cmdparams_save_str(&cmdparams, "tsout", &newstat);
  tsout = strcmp(kNOTSPECIFIED, tsseries);
  char *powseries = (char *)cmdparams_save_str(&cmdparams, "powout", &newstat);
  powout = strcmp(kNOTSPECIFIED, powseries);
  char *fftseries = (char *)cmdparams_save_str(&cmdparams, "fftout", &newstat);
  fftout = strcmp(kNOTSPECIFIED, fftseries);
  char *fft1series = (char *)cmdparams_save_str(&cmdparams, "fft1out", &newstat);
  fft1out = strcmp(kNOTSPECIFIED, fft1series);
  char *arpowseries = (char *)cmdparams_save_str(&cmdparams, "arpowout", &newstat);
  arpowout = strcmp(kNOTSPECIFIED, arpowseries);
  char *mavgseries = (char *)cmdparams_save_str(&cmdparams, "mavgout", &newstat);
  mavgout = strcmp(kNOTSPECIFIED, mavgseries);

  char *serieslist[6];
  enum seriesindex  {TSOUT, FFTOUT, FFT1OUT, POWOUT, ARPOWOUT, MAVGOUT};
  int  flagarr[6] = {tsout, fftout, fft1out, powout, arpowout, mavgout};
  int mfliparr[6] = {0, 0, 0, 0, 0, 0};
  int msignarr[6] = {0, 0, 0, 0, 0, 0};
  serieslist[TSOUT]=tsseries;
  serieslist[FFTOUT]=fftseries;
  serieslist[FFT1OUT]=fft1series;
  serieslist[POWOUT]=powseries;
  serieslist[ARPOWOUT]=arpowseries;
  serieslist[MAVGOUT]=mavgseries;
  long long histrecnumarr[6]={-1, -1, -1, -1, -1, -1};
  DRMS_RecordSet_t *outrecsetarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};
  DRMS_Record_t       *outrecarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};
  DRMS_Segment_t      *segoutarr[6]={NULL, NULL, NULL, NULL, NULL, NULL};

  /**** Sanity check of input parameters. ****/

  if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0)) 
  {
    fprintf(stderr, "ERROR: must specify an output dataseries\n");
    return 1;
  }

  if (navg <= 0) 
    navg=1;
  if (nskip <= 0) 
    nskip=1;
  if ((navg != 1) && (nskip != 1))
  {
    fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");
    return 1;
  }
  if (noff < 0 || noff >= nskip) 
  {
    fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
    return 1;
  }
  if (arpowout && !ifill)
  {
    fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");
    return 1;
  }
  if (fu_min_percent <= 0 || fu_min_percent > 100)
  {
    fprintf(stderr, "ERROR: FU_MIN_PERCENT must be > 0 and <= 100\n");
    return 1;
  }

  if (newstat) 
  {
    fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
    cpsave_decode_error(newstat);
    return 1;
  }  
  else if (savestrlen != strlen(savestr)) 
  {
    fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
    return 1;
  }


 // these must be present in the output dataseries and variable, not links or constants
 // 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...
  char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
  int is, ishold, itest, mflip;
  char *holdseries="";
/*
  char *cvsinfo;
  cvsinfo = (char *)malloc(1024);
  strcpy(cvsinfo,"$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.15 2013/07/02 20:33:16 tplarson Exp $");
  strcat(cvsinfo,"\n");
  strcat(cvsinfo,getdetrendversion());
  strcat(cvsinfo,"\n");
  strcat(cvsinfo,getgapfillversion());
*/
  for (is=0;is<6;is++)
  {

    if (!flagarr[is])
      continue;

    DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
                                                   serieslist[is],
                                                   DRMS_TRANSIENT, 
                                                   &status);

    if (status != DRMS_SUCCESS) 
    {
     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
     return 1;
    }

    DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);


    // set up ancillary dataseries for processing metadata
    if (histlink != NULL) 
    {
      if (strcmp(holdseries, histlink->info->target_series))
      {
        ishold=is;
        holdseries=strdup(histlink->info->target_series);
        histrecnumarr[is]=set_history(histlink);
        if (histrecnumarr[is] < 0)
        {
          fprintf(stderr, "ERROR: problem creating record in history dataseries for output dataseries %s\n", serieslist[is]);
          drms_close_record(tempoutrec, DRMS_FREE_RECORD);
          return 1;
        }
      }
      else
      {
        histrecnumarr[is]=histrecnumarr[ishold];
      }
    }
    else
    {
      fprintf(stderr,"WARNING: could not find history link in output dataseries %s\n", serieslist[is]);
    }

    for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
    {
      DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
      if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1)
      {
        fprintf(stderr, "ERROR: output keyword %s in series %s is either missing, constant, or a link\n", outchecklist[itest], serieslist[is]);
        drms_close_record(tempoutrec, DRMS_FREE_RECORD);
        return 1;
      }
    }

    mflip = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
    if (status == DRMS_SUCCESS)
      mfliparr[is]=mflip;    

    drms_close_record(tempoutrec, DRMS_FREE_RECORD);
  }

// for efficiency require that tsseries and fftseries have the same value of MFLIPPED.  this restriction may be lifted as noted below.
  if ((tsout && fftout) && (mfliparr[TSOUT] != mfliparr[FFTOUT]))
  {
    fprintf(stderr, "ERROR: series %s and %s must have the same value of MFLIPPED\n", tsseries, fftseries);
    return 1;
  }

  if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
  {
    DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status);
    if (status != DRMS_SUCCESS || gaprecset == NULL) 
    {
      fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
      return 1;
    }
    if (gaprecset->n == 0)
    {
      fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
      return 1;
    }
    if (gaprecset->n > 1)
    {
      fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
      return 1;
    } 
    gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);
    if (gapseg != NULL)
      gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status);
    if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
    {
      fprintf(stderr, "problem reading gaps segment: status = %d\n", status);
      return 1;
    }
    agaps=(int *)(gaparr->data);
    gapsize=gaparr->axis[0];
  }
  else
  {
    fits_get_img_size(fitsptr, 1, fdims, &fstatus);
    gapsize=fdims[0];
    agaps=(int *)(malloc(gapsize*sizeof(int)));
    fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus); 
    fits_close_file(fitsptr, &fstatus);
    if (fstatus) 
    {
      fits_report_error(stderr, fstatus);
      return 1;
    }
  }

  if (verbflag)
    printf("gapfile read, gapsize = %d\n",gapsize);

  data_dim=gapsize;
  if (n_sampled_out>gapsize) 
    data_dim=n_sampled_out;

  /* allocate temporary storage */
  gaps=(int *)(malloc(gapsize*sizeof(int)));
  gaps2=(int *)(malloc(gapsize*sizeof(int)));
  data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));
  ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float));
// now done below
//  if (fft1out || powout || mavgout || arpowout)
//    pow=(float *)(malloc(n_sampled_out*sizeof(float)));

  /* Read location of jumps. */
  if (!strcmp(sectionfile,kNOTSPECIFIED))
  {
    /* No sections file specified. Assume that the whole 
       time series in one section. */
    num_sections=1;
    sections = malloc(2*sizeof(int));
    sections[0] = 0;
    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);
    }
  }

  if (verbflag)
    printf("num_sections = %d\n",num_sections);

  DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);

  if (status != DRMS_SUCCESS || inrecset == NULL) 
  {
    fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
    return 1;
  }
  if (inrecset->n == 0)
  {
    fprintf(stderr, "ERROR: input recordset contains no records\n");
    return 1;
  }

  inrec=inrecset->records[0];
  char *inchecklist[] = {"T_START", "T_STOP", "QUALITY", "LMIN", "LMAX", "T_STEP"};
  for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
  {
    DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
    if (inkeytest == NULL)
    {
      fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
      drms_close_records(inrecset, DRMS_FREE_RECORD);
      return 1;
    }
  }

  cadence=drms_getkey_float(inrec, "T_STEP", NULL); //assume constant across all input records
  tsample=cadence;
  tstartsave=drms_getkey_time(inrec, "T_START", NULL); //must be constant across all input records unless GAPFILE is implemented differently

  int mflipin=drms_getkey_int(inrec, "MFLIPPED", &status);
  if (status != DRMS_SUCCESS) mflipin=0;

// already required mfliparr[TSOUT] == mfliparr[FFTOUT] above
  if (tsout && mfliparr[TSOUT] != mflipin)
    flipm=1;
  else if (fftout && mfliparr[FFTOUT] != mflipin)
    flipm=1;
  else
    flipm=0;

  for (is=2;is<5;is++)
  {
    if (!flagarr[is])
      continue;
    if (mfliparr[is] != (mflipin || flipm))
      msignarr[is]=-1;
    else
      msignarr[is]=1;
  }

  if (mflipin)
    msignarr[MAVGOUT]=1;
  else
    msignarr[MAVGOUT]=-1;

  // changed n_in to gapsize in following call
  if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections))
  {
    fprintf(stderr, "ERROR: problem in extract_gaps\n");
    return 1;
  }

  /* 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) 
    n_sampled_out = n_sampled_in;
  df1 = tsamplex*n_sampled_out;

  if (fftout || fft1out || powout || arpowout || mavgout)     
    plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, 
			    (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE);

  /* Set sum to zero before starting */
  if (mavgout) 
    msum = (float *)malloc((n_sampled_out/2)*sizeof(float));

  if (fft1out || powout || mavgout || arpowout)
  {
    /* Set first and last output index if not set as a parameter. */
    if (imax < imin) 
    {
      imin=0;
      imax=n_sampled_out/2-1;
    }
    npow=imax-imin+1;
    pow=(float *)(malloc(n_sampled_out*sizeof(float)));
    datahold=(float *)malloc(2*npow*sizeof(float));
  }

// needed to implement mfliparr[TSOUT] != mfliparr[FFTOUT] 
//  float *datatemp=(float *)malloc(2*n_sampled_out*sizeof(float));
  float *tmpptr=data;

  if (tsout || fftout) 
  {
    tslength[0]=tstotal[0]=2*n_sampled_out;
    tslength[1]=1;
    tsstart[0]=0;
    tsend[0]=2*n_sampled_out-1;
  }
  if (fft1out)
  {
    fft1length[0]=fft1total[0]=2*npow;
    fft1length[1]=1;
    fft1start[0]=0;
    fft1end[0]=2*npow-1;
  }
  if (powout || arpowout || mavgout)
  {
    powlength[0]=powtotal[0]=npow;
    powlength[1]=1;
    powstart[0]=0;
    powend[0]=npow-1;
  }
  instart[0]=0;
  inend[0]=2*gapsize - 1;

  status=drms_stage_records(inrecset, 1, 0);
  if (status != DRMS_SUCCESS)
  {
    fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status);
    return 1;
  }

  for (is=0; is<6; is++)
  {
    if (!flagarr[is])
      continue;

    outrecsetarr[is] = drms_create_records(drms_env, inrecset->n, serieslist[is], DRMS_PERMANENT, &status);

    if (status != DRMS_SUCCESS || outrecsetarr[is] == NULL) 
    {
     fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
     return 1;
    }
  }

  int irec;
  for (irec=0;irec<inrecset->n;irec++)
  {

    inrec = inrecset->records[irec];

    lmin=drms_getkey_int(inrec, "LMIN", NULL);
    lmax=drms_getkey_int(inrec, "LMAX", NULL);
    if (lmin != lmax)
    {
      fprintf(stderr, "ERROR: lmin != lmax not yet supported.\n");
      return 0;
    }
    lmode=lmin;

    if (verbflag > 1) 
    {
      printf("starting l=%d\n", lmode);
    }

    tstart=drms_getkey_time(inrec, "T_START", NULL);
    tstop =drms_getkey_time(inrec, "T_STOP", NULL);
    cadence=drms_getkey_float(inrec, "T_STEP", NULL);
    if (tstart != tstartsave)
    {
// to lift this restriction must be able to specify multiple gap files/records as input
      sprint_time(tstartstr, tstart, "TAI", 0);
      fprintf(stderr, "ERROR: all input records must have same T_START, record %d has %s\n", irec, tstartstr);
      return 0;
    }

    n_in=(tstop-tstart)/cadence;
    if (n_in != gapsize)
    {
      fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, n_in=%d\n", gapsize, n_in);
      return 0;
    }

    tstartout=tstart+ifirst*tsample;
    tstopout=tstartout+df1;
    sprint_time(tstartstr, tstartout, "TAI", 0);

    if (seginflag)
      segin = drms_segment_lookup(inrec, segnamein);
    else
      segin = drms_segment_lookupnum(inrec, 0);
    if (segin == NULL)
    {
      fprintf(stderr, "ERROR: problem looking up input segment: recnum = %lld\n", inrec->recnum);
      return 0;
    }

    for (is=0; is<6; is++)
    {
      if (!flagarr[is])
        continue;
/*
      outrecarr[is] = drms_create_record(drms_env, serieslist[is], DRMS_PERMANENT, &status);
      if (status != DRMS_SUCCESS) 
      {
       fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status);
       return 0;
      }
*/
      outrec=outrecsetarr[is]->records[irec];
      if (histrecnumarr[is] > 0)
        drms_setlink_static(outrec, histlinkname,  histrecnumarr[is]);
      drms_setlink_static(outrec, srclinkname,  inrec->recnum);

      if (segoutflag)
        segoutarr[is] = drms_segment_lookup(outrec, segnameout);
      else
        segoutarr[is] = drms_segment_lookupnum(outrec, 0);
      if (segoutarr[is] == NULL)
      {
        fprintf(stderr, "ERROR: problem looking up outputput segment in series %s\n", serieslist[is]);
        return 0;
      }

// the following is not needed if at least one of the first N-1 dimensions are 0 in the jsd, 
// or if the first N-1 dimensions in the jsd are always the dimensions desired.
// it *is* needed any time one must override the jsd defaults
/*
      switch(is)
      {
        case TSOUT:
        case FFTOUT:
          segoutarr[is]->axis[0]=2*n_sampled_out;
          segoutarr[is]->axis[1]=lmode+1;
          break;
        case FFT1OUT:
          segoutarr[is]->axis[0]=2*npow;
          segoutarr[is]->axis[1]=2*lmode+1;
          break;
        case POWOUT:
        case ARPOWOUT:
          segoutarr[is]->axis[0]=npow;
          segoutarr[is]->axis[1]=2*lmode+1;
          break;
        case MAVGOUT:
          segoutarr[is]->axis[0]=npow;
          segoutarr[is]->axis[1]=1;
          break;
        default:
          ;
      }
*/
      switch(is)
      {
        case TSOUT:
        case FFTOUT:
          tstotal[1]=lmode+1;
          break;
        case FFT1OUT:
          fft1total[1]=2*lmode+1;
          break;
        case POWOUT:
        case ARPOWOUT:
          powtotal[1]=2*lmode+1;
          break;
        case MAVGOUT:
        default:
          ;
      }

    }

    if (mavgout)
      for (i=0;i<n_sampled_out/2;i++) 
        msum[i] = 0.0;

  /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
    for (m=0; m<=lmode; m++) 
    {

      if (verbflag > 2) 
      {
        printf("starting m=%d\n", m);
      }

      instart[1]=m;
      inend[1]=m;
      inarr = drms_segment_readslice(segin, usetype, instart, inend, &status);
      if (status != DRMS_SUCCESS || inarr == NULL )
      {
        fprintf(stderr, "ERROR: problem reading input segment: status = %d, recnum = %lld\n", status, inrec->recnum);
        return 0;
      }
      in_data=(float *)(inarr->data);

    /* Extracts data copy */
      extract_data(n_sampled_in, gaps, in_data, data);

    /* Detrend entire time series if required */
      if (detrend_const != 0) 
      {
        detrend_order = floor(detrend_length/detrend_const)+2;
        cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps, 
			 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) 
      {
        cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2,
		      fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, 
		      &ar_order, (_Complex float *)ar_coeff);
      }

      drms_free_array(inarr);

      memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));
      if ((ifirst+n_sampled_out) >= n_sampled_in) 
        memset(&data[2*(n_sampled_in-ifirst)], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));      

      if (tsout)
      {
// code to flip m on this output separately.  in the present code this has already been accomplished by setting flipm.
/*
        if (mfliparr[TSOUT] != mflipin)
        {
          for (i = 0; i < n_sampled_out; i++)
          {
            datatemp[2*i]=data[2*i];
            datatemp[2*i+1]=-data[2*i+1];
          }
          tmpptr=datatemp;
        }
        else
          tmpptr=data;
*/
        tsstart[1]=m;
        tsend[1]=m;
        outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
        outarr->bzero=segoutarr[TSOUT]->bzero;
        outarr->bscale=segoutarr[TSOUT]->bscale;
        status=drms_segment_writeslice_ext(segoutarr[TSOUT], outarr, tsstart, tsend, tstotal, 0);
        free(outarr);
        if (status != DRMS_SUCCESS)
        {
          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]);
          return 0;
        }
      }

      if (fftout || fft1out || powout || mavgout)
        fftwf_execute(plan);

      if (fftout)
      {
// code to flip m on this output separately.  in the present code this has already been accomplished by setting flipm.
/*
        if (mfliparr[FFTOUT] != mflipin)
        {
          datatemp[0]=data[0];
          datatemp[1]=-data[1];
          for (i = 1; i < n_sampled_out; i++)
          {
            datatemp[2*i]=data[2*(n_sampled_out-i)];
            datatemp[2*i+1]=-data[2*(n_sampled_out-i)+1];
          }
          tmpptr=datatemp;
        }
        else
          tmpptr=data;
*/
        tsstart[1]=m;
        tsend[1]=m;
        outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
        outarr->bzero=segoutarr[FFTOUT]->bzero;
        outarr->bscale=segoutarr[FFTOUT]->bscale;
        status=drms_segment_writeslice_ext(segoutarr[FFTOUT], outarr, tsstart, tsend, tstotal, 0);
        free(outarr);
        if (status != DRMS_SUCCESS)
        {
          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]);
          return 0;
        }
      }

      if (fft1out)
      {

        fft1start[1]=lmode+m*msignarr[FFT1OUT];
        fft1end[1]=lmode+m*msignarr[FFT1OUT];
        outarr = drms_array_create(usetype, 2, fft1length, data+2*imin, &status);
        outarr->bzero=segoutarr[FFT1OUT]->bzero;
        outarr->bscale=segoutarr[FFT1OUT]->bscale;
        status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
        free(outarr);
        if (status != DRMS_SUCCESS)
        {
          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]);
          return 0;
        }

        if (m > 0) 
        {
          /* Do negative m */
          for (i=0; i<npow; i++) 
          {
            /* Conjugate for negative m */
            datahold[2*i] = data[2*(n_sampled_out-imin-i)];
            datahold[2*i+1] = -data[2*(n_sampled_out-imin-i)+1];
          }
          fft1start[1]=lmode-m*msignarr[FFT1OUT];
          fft1end[1]=lmode-m*msignarr[FFT1OUT];
          outarr = drms_array_create(usetype, 2, fft1length, datahold, &status);
          outarr->bzero=segoutarr[FFT1OUT]->bzero;
          outarr->bscale=segoutarr[FFT1OUT]->bscale;
          status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
          free(outarr);
          if (status != DRMS_SUCCESS)
          {
            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]);
            return 0;
          }
        }
      }

      if (powout || mavgout)
        for (i = 0; i < n_sampled_out; i++) 
          pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1];

      if (powout)
      {
        powstart[1]=lmode+m*msignarr[POWOUT];
        powend[1]=lmode+m*msignarr[POWOUT];
        outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status);
        outarr->bzero=segoutarr[POWOUT]->bzero;
        outarr->bscale=segoutarr[POWOUT]->bscale;
        status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
        free(outarr);
        if (status != DRMS_SUCCESS)
        {
          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]);
          return 0;
        }

        if (m > 0) 
        {
          /* Do negative m */
          if (imin)
            datahold[0] = pow[n_sampled_out-imin];
          else
            datahold[0] = pow[0];
          for (i=1; i<npow;i++) 
            datahold[i] = pow[n_sampled_out-imin-i];

          powstart[1]=lmode-m*msignarr[POWOUT];
          powend[1]=lmode-m*msignarr[POWOUT];
          outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
          outarr->bzero=segoutarr[POWOUT]->bzero;
          outarr->bscale=segoutarr[POWOUT]->bscale;
          status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
          free(outarr);
          if (status != DRMS_SUCCESS)
          {
            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]);
            return 0;
          }
        }
      }

      if (mavgout)
      {
        mshift = floor(df1*splitting(lmode,m*msignarr[MAVGOUT])+0.5f);
        if (mshift >= 0) 
          for (i=0; i<n_sampled_out/2-mshift; i++)
            msum[i] += pow[mshift+i];
        else
          for (i=0; i<n_sampled_out/2+mshift; i++)
            msum[i-mshift] += pow[i];
        if (m > 0) 
        {
          /* Do negative m */
          mshift = floor(df1*splitting(lmode,-m*msignarr[MAVGOUT])+0.5f);
          if (mshift >=0)
            for (i=1; i<n_sampled_out/2-mshift;i++)
              msum[i] += pow[n_sampled_out-mshift-i];
          else
            for (i=1; i<n_sampled_out/2+mshift; i++)
              msum[i-mshift] += pow[n_sampled_out-i];
        }
      }

      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));
        fftwf_execute(plan);
        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]);

        powstart[1]=lmode+m*msignarr[ARPOWOUT];
        powend[1]=lmode+m*msignarr[ARPOWOUT];
        outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status);
        outarr->bzero=segoutarr[ARPOWOUT]->bzero;
        outarr->bscale=segoutarr[ARPOWOUT]->bscale;
        status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
        free(outarr);
        if (status != DRMS_SUCCESS)
        {
          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]);
          return 0;
        }

        if (m > 0) 
        {
          /* Do negative m */
          if (imin)
            datahold[0] = pow[n_sampled_out-imin];
          else
            datahold[0] = pow[0];
          for (i=1; i<npow;i++) 
            datahold[i] = pow[n_sampled_out-imin-i];

          powstart[1]=lmode-m*msignarr[ARPOWOUT];
          powend[1]=lmode-m*msignarr[ARPOWOUT];
          outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
          outarr->bzero=segoutarr[ARPOWOUT]->bzero;
          outarr->bscale=segoutarr[ARPOWOUT]->bscale;
          status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
          free(outarr);
          if (status != DRMS_SUCCESS)
          {
            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]);
             return 0;
          }
        }
      }

    } //end of loop over m

    if (mavgout)
    {
      c=1.0f/(2*lmode+1);
      for (i=0; i<npow; i++) 
        datahold[i] = msum[imin+i] * c;
      outarr = drms_array_create(usetype, 2, powlength, datahold, &status);
      outarr->bzero=segoutarr[MAVGOUT]->bzero;
      outarr->bscale=segoutarr[MAVGOUT]->bscale;
      status=drms_segment_write(segoutarr[MAVGOUT], outarr, 0);
      free(outarr);
      if (status != DRMS_SUCCESS)
      {
        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);
        return 0;
      }
    }

    for (is=0;is<6;is++)
    {
      if (!flagarr[is])
        continue;

      outrec=outrecsetarr[is]->records[irec];
      drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);
      drms_setkey_int(outrec, "QUALITY", 0);
      drms_setkey_time(outrec, "T_START", tstartout);
      drms_setkey_time(outrec, "T_STOP", tstopout);
      drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2);
      drms_setkey_float(outrec, "T_STEP", tsamplex);
      drms_setkey_int(outrec, "NDT", n_sampled_out);
      if (tagflag)
        drms_setkey_string(outrec, "TAG", tag);
      if (verflag)
         drms_setkey_string(outrec, "VERSION", version);

      tnow = (double)time(NULL);
      tnow += UNIX_epoch;
      drms_setkey_time(outrec, "DATE", tnow);

//      drms_close_record(outrec, DRMS_INSERT_RECORD);

    }

  }

  drms_close_records(inrecset, DRMS_FREE_RECORD);
  for (is=0;is<6;is++)
  {
    if (!flagarr[is])
      continue;
    drms_close_records(outrecsetarr[is], DRMS_INSERT_RECORD);
  }

  if(ar_coeff != NULL)
    free(ar_coeff);

  if (fftout || fft1out || powout || arpowout || mavgout)   
    fftwf_destroy_plan(plan); 

  wt=getwalltime();
  ct=getcputime(&ut, &st);
  if (verbflag) 
  {
    printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", 
            wt-wt0, ct-ct0);
  }

  printf("module %s successful completion\n", cmdparams.argv[0]);
  return 0;
}


static double splitting(int l, int m)
{
  double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0};
  double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8};
  double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5};
  /* f0 is the frequency used for generating the above a-coeff */
  double f0 = 1500.0;
  /* fcol is the frequency for which to optimize the collaps */
  double fcol = 800.0;

  double ll,x,x3,x5,a1,a2,a3,a5,frac,lx;
  int ix;

  if (l == 0) 
    ll = 1; 
  else 
    ll = sqrt(l*(l+1.));
  x = m/ll;
  x3 = x*x*x;
  x5 = x3*x*x;
  lx = 5*log10(l*f0/fcol)-4;
  ix = floor(lx);
  frac = lx-ix;
  if (lx <= 0) {
    ix = 0;
    frac = 0.0;
  }
  if (lx >=8) {
    ix = 7;
    frac = 1.0;
  }
  a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1];
  a2 = 0.0;
  a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1];
  a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1];

  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;
}



/* Extract the data samples actually used by skipping or averaging
   data. Replace missing data that are not marked as gaps with zero. */
void extract_data(int n_in, int *gaps, float *data_in, float *data_out)
{
  int i,j, nmiss = 0;
// the following check will be skipped in production executables since they will be built with NDEBUG defined.
// it doesn't matter because the check is already done in DoIt().
  assert(nskip==1 || navg==1);
  if ((navg == 1) && (nskip == 1)) 
  {
    for (i=0; i<n_in; i++) 
    {
      if (gaps[i]==0 )
      {
	data_out[2*i] = 0.0f;
	data_out[2*i+1] = 0.0f;
      }
      else
      {
	if (isnan(data_in[2*i]) || isnan(data_in[2*i+1]))
	{
	  data_out[2*i] = 0.0f;
	  data_out[2*i+1] = 0.0f;
	  gaps[i] = 0;
	  nmiss++;
	}
	else
	{
	  data_out[2*i] = data_in[2*i];
	  data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1];
	}
      }
    }
  }
  else if (nskip != 1) 
  {
    for (i = 0; i<n_in; i++) 
    {
      if (gaps[i]==0 )
      {
	data_out[2*i] = 0.0f;
	data_out[2*i+1] = 0.0f;
      }
      else
      {
	int ix = nskip*i+noff;
	if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
	{
	  data_out[2*i] = 0.0f;
	  data_out[2*i+1] = 0.0f;
	  gaps[i] = 0;
	  nmiss++;
	}
	else
	{
	  data_out[2*i] = data_in[2*ix];
	  data_out[2*i+1] =  flipm ? -data_in[2*ix+1] : data_in[2*ix+1];
	}
      }
    }
  }
  else if (navg != 1) 
  {
    for (i = 0; i<n_in; i++) 
    {
      if (gaps[i]==0 )
      {
	data_out[2*i] = 0.0f;
	data_out[2*i+1] = 0.0f;
      }
      else
      {
	float avgr = 0.0f;
	float avgi = 0.0f;
	for (j=0; j<navg; j++) 
	{
	  int ix = navg*i+j;
	  if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1]))
	  {
	    data_out[2*i] = 0.0f;
	    data_out[2*i+1] = 0.0f;
	    gaps[i] = 0;
	    nmiss++;
	    break;
	  }
	  else 
	  {
	    avgr += data_in[2*ix];
	    avgi += data_in[2*ix+1];
	  }
	}
	if (j == navg) 
	{
	  data_out[2*i] = avgr/navg;
	  data_out[2*i+1] = avgi/navg;
	}
      }
    }  
  }
}

/* Extract the windows function for the samples actually used after
   subsampling or averaging. Modify the list of continous sections
   accordingly, and make sure we do not average across section boundaries. */
static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, 
		  float tsample_in, float *tsample_out, 
		  int *num_sections, int *sections)
{
// the following check will be skipped in production executables since they will be built with NDEBUG defined.
// it doesn't matter because the check is already done in DoIt().
  assert(nskip==1 || navg==1);
  int i,j,sect, start,stop;


  if (*num_sections<1)
  {
    fprintf(stderr,"Apparantly no sections of data available.");
    return 1;
  }
  /* Mask out data that in between sections if this hasn't already been
     done in gapfile. */
  for (i=0; i<sections[0] && i<n_in; i++)
    gaps_in[i] = 0;
  for (sect=1; sect<(*num_sections); sect++)
  {
    for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++)
      gaps_in[i] = 0;
  }
  for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++)
    gaps_in[i] = 0;
    

  if ((navg == 1) && (nskip == 1)) 
  {
    *n_out = n_in;
    *tsample_out = tsample_in;
    for (i=0; i<*n_out; i++) 
      gaps_out[i] = gaps_in[i];
  }
  else if (nskip != 1) 
  {
    *n_out = n_in/nskip;
    *tsample_out = nskip*tsample_in;
    for (i=0; i<*n_out; i++) 
    {
      gaps_out[i] = gaps_in[nskip*i+noff];
    }
    for (i=0; i<2*(*num_sections); i+=2)
    {
      start = (int)ceilf(((float)(sections[i]-noff))/nskip);
      stop = (int)floorf(((float)(sections[i+1]-noff))/nskip);
      if (start <= stop)
      {
	sections[i] = start;
	sections[i+1] = stop;
      }
      else 
      {
	/* This section was skipped entirely. */
	for (j=i; j< 2*(*num_sections-1); j+=2) 
	{
	  sections[j] = sections[j+2];
	  sections[j+1] = sections[j+3];
	}
	*num_sections -= 1;
      }
    }
  }
  else  if (navg != 1) 
  {
    sect = 0;
    *n_out = n_in/navg;
    *tsample_out = tsample*navg;
    for (i=0; i<*n_out; i++) 
    {      
      int igx = 1;
      while (sect < *num_sections && 
	     !(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1]))
	sect++;

      /* Don't allow averaging of data from different sections. */
      if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1])
      {
	for (j=0; j<navg; j++) 
	  igx *= gaps_in[navg*i+j];
	gaps_out[i] = igx;
      }
      else
	gaps_out[i] = 0;
    }
    for (i=0; i<2*(*num_sections); i+=2)
    {
      start = (int)ceilf(((float)sections[i])/navg);
      stop = (int)floorf(((float)sections[i+1])/navg);
      if (start <= stop)
      {
	sections[i] = start;
	sections[i+1] = stop;
      }
      else 
      {
	/* This section was skipped entirely. */
	for (j=i; j< 2*(*num_sections-1); j+=2) 
	{
	  sections[j] = sections[j+2];
	  sections[j+1] = sections[j+3];
	}
	*num_sections -= 1;
      }
    }
  }  
  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;
}