proj/globalhs/apps/jtsslice.c

/*
this module creates power spectra out of slices of timeseries
based on jtsfiddle, but no detrending or gapfilling
*/

/* 
 *  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 = "jtsslice";
char *cvsinfo_jtsslice = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsslice.c,v 1.12 2013/04/28 08:05:21 tplarson Exp $";

/* Command line arguments: */
ModuleArgs_t module_args[] =
{
  {ARG_STRING,   "in", "", "input data records"},
  {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,   "GAPFILE", "", "record or file containing window function"},

  {ARG_INT,      "NDT", "-1", "", ""},
  {ARG_STRING,   "TTOTAL",   NULL,          "total length of time in output"},

  {ARG_INT,      "IFIRST", "0", "", ""},
//  {ARG_INT,      "FLIPM", "0", "", ""},

  {ARG_INT,      "TSOUT", "0", "", ""},
  {ARG_INT,      "FFTOUT", "0", "", ""},
  {ARG_INT,      "FFT1OUT", "0", "", ""},
  {ARG_INT,      "POWOUT", "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"

#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. */
static char *gapfile;
static int lmode, n_sampled_out, flipm;
static float tsample;
static int ifirst, tsout,fftout, fft1out, powout, 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);

/************ Here begins the main module **************/
int DoIt(void)
{
  int i;
  int start_index[2], counts[2], intervals[2];
  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;
  float *msum, *pow;
  float *data, *out_data, *in_data;
  fftwf_plan plan;
  int num_sections, *sections;
  char tstartstr[100];

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

  int length[2], startind[2], endind[2];
  float *datarr;
  int status=DRMS_SUCCESS;
  int newstat=0;
  char *ttotal;
  double nseconds;

  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. */
  gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
  n_sampled_out = cmdparams_save_int  (&cmdparams, "NDT", &newstat);
  ifirst      = cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
//  flipm       = cmdparams_save_int    (&cmdparams, "FLIPM", &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);
  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;

  ttotal=(char *)cmdparams_save_str(&cmdparams, "TTOTAL", &newstat);
  status=drms_names_parseduration(&ttotal, &nseconds, 1);
  if (status != DRMS_SUCCESS)
  {
    fprintf(stderr, "ERROR: problem parsing TTOTAL, = %s\n", ttotal);
    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;
  }


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

  /* Default is to just output the timeseries. */
  if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0)) 
    tsout=1;
  /* Only one type of output allowed at a time */
  if ((tsout+fftout+fft1out+powout+mavgout) !=1) 
  {
    fprintf(stderr, "ERROR: only one type of output allowed at a time\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;
  }


  DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
                           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);
    return 1;
  }

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

//  char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsslice.c,v 1.12 2013/04/28 08:05:21 tplarson Exp $");
  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   
  char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};

  int itest;
  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 is either missing, constant, or a link\n", outchecklist[itest]);
      drms_close_record(tempoutrec, DRMS_FREE_RECORD);
      return 1;
    }
  }

  cadence=drms_getkey_float(tempoutrec, "T_STEP", &status);
  double chunksecs = n_sampled_out*cadence;
  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);
    return 1;
  }
  int ntimechunks = nseconds/chunksecs;
  tsample=cadence;

  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))
  {
    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, "ERROR: 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) 
    {
      fprintf(stderr, "ERROR: ");
      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;

  num_sections=1;
  sections = malloc(2*sizeof(int));
  sections[0] = 0;
  sections[1] = data_dim-1;

  /* allocate temporary storage */
  gaps=(int *)(malloc(gapsize*sizeof(int)));
  data=(float *)(fftwf_malloc(2*data_dim*sizeof(float)));

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

  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;

  if (mflipout != mflipin)
    flipm=1;
  else
    flipm=0;

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

  if (n_sampled_out == -1) 
    n_sampled_out = n_sampled_in;
  df1 = tsamplex*n_sampled_out;

  if (fftout || fft1out || powout || 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 *)calloc(n_sampled_out/2,sizeof(float));

  if (fft1out || powout || mavgout)
  {
    /* 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)));
  }


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

  int nrecsout = ntimechunks * inrecset->n;
  DRMS_RecordSet_t *outrecset = drms_create_records(drms_env, nrecsout, outseries, DRMS_PERMANENT, &status);
  if (status != DRMS_SUCCESS || outrecset == NULL) 
  {
   fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
   return 1;
  }

  int gapoffset, gapoffset0 = ifirst;
  double tstart0 = tstartsave+ifirst*tsample;
  int irec;
  for (irec=0;irec<inrecset->n;irec++)
  {

    inrec = inrecset->records[irec];

    int lmin=drms_getkey_int(inrec, "LMIN", NULL);
    int 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)
      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 time-series, gapsize=%d, n_in=%d\n", gapsize, n_in);
      return 0;
    }

    /* Create output data set. */
    if (tsout || fftout) 
    {
      length[0]=2*n_sampled_out;
      length[1]=lmode+1;
    }
    else 
    {
      if (fft1out) 
        length[0]=2*npow;
      else 
        length[0]=npow;
      if (powout || fft1out) 
        length[1]=2*lmode+1;
      else 
        length[1]=1;
    }

    startind[1]=0;
    endind[1]=lmode;

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

    int iset;
    for (iset=0; iset < ntimechunks; iset++)
    {
      tstartout=tstart0 + iset * chunksecs;
      tstopout=tstartout+chunksecs;
      sprint_time(tstartstr, tstartout, "TAI", 0);
      gapoffset = gapoffset0 + iset*n_sampled_out;

      if (verbflag>1)
      {
        printf("  starting tstart = %s, ",tstartstr);
        if (irec == 0)
        {
          int ig, gtotal=0;
          for (ig=0;ig<n_sampled_out;ig++)
            gtotal+=gaps[gapoffset+ig];
          float dutycycle = (float)gtotal/n_sampled_out;
          printf("dutycycle = %f\n", dutycycle);
        }
      }

// set ifirst to gapoffset for reading slice
      ifirst=gapoffset;
      startind[0]=2*(ifirst);
      endind[0]=2*(ifirst + n_sampled_out) - 1;
// set ifirst=0 because data is read starting at ifirst by drms_segment_readslice
      ifirst=0;

      inarr = drms_segment_readslice(segin, usetype, startind, endind, &status);
      if (status != DRMS_SUCCESS || inarr == NULL)
      {
        fprintf(stderr, "ERROR: problem reading input segment: recnum = %ld, status = %d\n", inrec->recnum, status);
        return 0;
      }
      datarr=(float *)(inarr->data);

/*
      outrec = drms_create_record(drms_env, outseries, lifetime, &status);
      if (status != DRMS_SUCCESS || outrec==NULL)
      {
         fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);
         return 0;
      }
*/
//      outrec = outrecset->records[ntimechunks*irec+iset];
      outrec = outrecset->records[iset*inrecset->n + irec];

      if (histlink != NULL)
        drms_setlink_static(outrec, histlinkname,  histrecnum);
      if (srclink != NULL)
        drms_setlink_static(outrec, srclinkname,  inrec->recnum);

      if (segoutflag)
        segout = drms_segment_lookup(outrec, segnameout);
      else
        segout = drms_segment_lookupnum(outrec, 0);
      outarr = drms_array_create(usetype, 2, length, NULL, &status);
      if (status != DRMS_SUCCESS || outarr == NULL || segout == NULL)
      {
        fprintf(stderr,"ERROR: problem creating output array or segment: length = [%d, %d], status = %d", length[0], length[1], status);
        return 0; 
      }
      out_data = outarr->data;

  /***** 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);

        in_data=datarr + m*2*n_sampled_out;

        /* Extracts data copy */
        extract_data(n_sampled_out, gaps+gapoffset, in_data, data);

        /* pad with zeros if the last point output (n_sampled_out+ifirst) 
           is after the last data points read in (nread). */
        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], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));      

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

        if (tsout || fftout) 
          memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));

        else if (fft1out) 
        {
        /* Do positive m */
          memcpy(&out_data[2*(lmode+m)*npow], &data[2*imin], 2*npow*sizeof(float));
          if (m > 0) 
          {
          /* Do negative m */
        for (i=0; i<npow; i++) 
        {
           /* Conjugate for negative m */
          out_data[2*((lmode-m)*npow+i)] = data[2*(n_sampled_out-imin-i)];
          out_data[2*((lmode-m)*npow+i)+1] = -data[2*(n_sampled_out-imin-i)+1];
        }
          }
        }
        else 
        {
        /* Compute power spectrum from complex FFT. */      
          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) 
          {
          /* Do positive m */
            memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));
        if (m > 0) 
            {
            /* Do negative m */
          if (imin)
            out_data[(lmode-m)*npow] = pow[n_sampled_out-imin];
          else
            out_data[(lmode-m)*npow] = pow[0];
          for (i=1; i<npow;i++) 
            out_data[(lmode-m)*npow+i] = pow[n_sampled_out-imin-i];
            }
          }
          else 
          {
          /* m-averaged output */
          /* Sum all freqs, not only those in output. Should be fixed. */
          /* Maybe should allow for wrapping around Nyquist. */
          /* Do positive m */
            mshift = floor(df1*splitting(lmode,m)+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)+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];
            }
          }
        }

      } /* End of m loop */

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

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

      if (mavgout) 
      {
        c=1.0f/(2*lmode+1);
        for (i=0; i<npow; i++) 
          out_data[i] = msum[imin+i] * c;
      }

      outarr->bzero=segout->bzero;
      outarr->bscale=segout->bscale;
      status=drms_segment_write(segout, outarr, 0);
      if (status != DRMS_SUCCESS)
      {
        fprintf(stderr, "ERROR: problem writing output segment: status = %d, T_START = %s, LMODE = %d, histrecnum = %lld\n", status, tstartstr, lmode, histrecnum);
        return 0;
      }
      drms_free_array(outarr);
      drms_free_array(inarr);
//    drms_close_record(outrec, DRMS_INSERT_RECORD);

    }

  }

  drms_close_records(inrecset,  DRMS_FREE_RECORD);
  drms_close_records(outrecset, DRMS_INSERT_RECORD);

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

  wt=getwalltime();
  ct=getcputime(&ut, &st);
  if (verbflag) 
  {
    fprintf(stdout, "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;
}

---