/* this JSOC module is a port of an SOI module.
 * ported by Tim Larson.
 */

/* 
 *  peakbagging program - pkbgn24.c 
 */

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

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

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

char *module_name = "jpkbgn";
char *cvsinfo_jpkbgn = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jpkbgn.c,v 1.12 2015/12/30 00:18:25 tplarson Exp $";

ModuleArgs_t module_args[] =
{
   {ARG_STRING, "in", "", "input data records"},
   {ARG_STRING, "out", "", "output data series"},
   {ARG_STRING, "ffttmp", "", ""},
   {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,   "LOGFILE", "stdout", ""},
   {ARG_INT,      "READFFT", "0", "", ""},
   {ARG_STRING,   "GAPFILE", "", "record or file containing window function"},
   {ARG_INT,      "IPAR", "1", "", ""},
   {ARG_STRING,   "PARAMFILE", "nofile", ""},
   {ARG_STRING,   "PAR1FILE", "", ""},
   {ARG_STRING,   "RESFILE", "result", ""},
   {ARG_STRING,   "CROSSFILE", "", ""},
   {ARG_INT,      "NDELTAL", "6", "", ""},
   {ARG_INT,      "NOISETYPE", "1", "", ""},
   {ARG_STRING,   "CROSSNFILE", "/dev/null", ""},
   {ARG_INT,      "NNOIB", "1", "", ""},

   {ARG_STRING,   "NOIB", "40000 43000", ""},

   {ARG_INT,      "NNOIM", "0", "", ""},
   {ARG_STRING,   "NOIM", " ", ""},
   {ARG_FLOAT,    "CSUBM", "0.5", "", ""},
   {ARG_INT,      "IFOLLOW", "1", "", ""},
   {ARG_FLOAT,    "ANOISE", "0.0", "", ""},
   {ARG_FLOAT,    "AOTHER", "0.0", "", ""},
//   {ARG_INT,      "LMODE", "", "", ""},
   {ARG_INT,      "IODD", "1", "", ""},
   {ARG_INT,      "NMIN", "0", "", ""},
   {ARG_INT,      "NMAX", "50", "", ""},
   {ARG_FLOAT,    "FMIN", "1000", "", ""},
   {ARG_FLOAT,    "FMAX", "4000", "", ""},
   {ARG_INT,      "IMFREQ", "0", "", ""},
   {ARG_INT,      "ICASE", "3", "", ""},
   {ARG_INT,      "NACOEFF", "5", "", ""},
   {ARG_INT,      "NDT", "51840", "", ""},
//   {ARG_FLOAT,    "TSAMPLE", "60", "", ""},
   {ARG_INT,      "IDF", "1", "", ""},
   {ARG_INT,      "NDF0", "20", "", ""},
   {ARG_INT,      "NDF1", "50", "", ""},
   {ARG_FLOAT,    "CDF", "5", "", ""},
   {ARG_INT,      "NDM", "10", "", ""},
   {ARG_INT,      "NDL", "6", "", ""},
   {ARG_FLOAT,    "XSAFE", "10.0", "", ""},
   {ARG_FLOAT,    "XSAFE1", "1.0", "", ""},
   {ARG_FLOAT,    "DELMIN", "0.01", "", ""},
   {ARG_INT,      "MAXCHI", "50", "", ""},
   {ARG_INT,      "MAXGRAD", "30", "", ""},
   {ARG_INT,      "MAXCOF", "4", "", ""},
   {ARG_INT,      "IFIX", "0", "", ""},
   {ARG_INT,      "IFIRST", "0", "", ""},
   {ARG_INT,      "FLIPM", "0", "", ""},
   {ARG_INT,      "NLOBES", "", "", ""},
   {ARG_INT,      "NFOUR", "0", "", ""},
   {ARG_INT,      "GONGFLAG", "0", "", ""},
   {ARG_INT,      "PRINTFIT", "0", "", ""},
   {ARG_FLOAT,    "CDETREND", "50.0", "", ""},
   {ARG_INT,      "IFILL", "1", "", ""},
   {ARG_INT,      "FDIFF", "0", "", ""},
   {ARG_INT,      "IASYM", "0", "", ""},
   {ARG_INT,      "ICT", "0", "", ""},
   {ARG_INT,      "IWOOD", "0", "", ""},
   {ARG_FLOAT,    "WOODB1", "0.0", "", ""},
   {ARG_FLOAT,    "WOODB2", "0.0", "", ""},
   {ARG_INT,      "CTX", "1.0", "", ""},

   {ARG_END}
};

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

void detrend(float *data, float *gaps, int length, float cdetrend);
void fill_gaps(float *data, float *gaps, float *newgaps, int length);
//char *getpkbgapsversion(void);

void sub24_(char *logfile,char *fftfile, float *rgaps, int *ipar, char *paramfile, char *par1file, char *resfile,
       char *crossfile, int *ndeltal, int *noisetype, char *crossnfile, int *nnoib, int *noib, int *nnoim, int *noim,
       float *csubm, int *ifol, float *anoise, float *aother, int *lin, int *ioddin, int *iasymin, int *ictin, float *ctxin,
       int *iwoodin, float *woodb1, float *woodb2,  int *nmin, int *nmax, float *fmin, float *fmax,
       int *imfin, int *icasein, int *nain, int *idtin, int *ndtin, float *tsin, int *idf, int *ndf0, int *ndf1,
       float *cdf, int *ndmin, int *ndl,
       float *xsin, float *xs1in, float *delx, int *maxchi, int *maxgrad, int *maxcof, int *recin,
       int *nlobesin, int *nfourin, int *gongin, int *printflag,
       int, int, int, int, int, int, int);

void cffti_(int *, float *);
void cfftb_(int *, float *, float *);

//void getpkbgnversion_(char *, int);


int DoIt(void)
{
  int status = DRMS_SUCCESS;
  int newstat = 0;
  int noib[10], noim[10];

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

  int zero=0;
  int i, j;
  float c,sum;
  char *ptr;

  FILE *fileptr;

  int readfft, savefft;
  int in_nsets;

  int fdiff,flipm,m,navail,nread,nx1;
  float isign,*data;
  TIME epoch, step, start, stop;
  int gapsize, istart, istop;
  char *ffttmp, *fftfull, *infft;
  float cdetrend;
  char *syscall;

  char *gapfile, *logfile, *paramfile, *par1file, *resfile, *crossfile, *crossnfile;
  int loglen, fftlen, paramlen, par1len, reslen, crosslen, crossnlen;
  int ipar, ndeltal, noisetype, nnoib, nnoim;
  char *noibstr, *noimstr;
  float csubm, anoise, aother, tsample, cdf, xsafe, xsafe1, delmin, fmin, fmax;
  int ifollow, lmode, iodd, nmin, nmax, imfreq, icase, nacoeff;
  int ndt, idf, ndf0, ndf1, ndm, ndl, maxchi, maxgrad, maxcof;

  int ifix, ifirst, ifill;
  int nlobes,nfour,gongflag,printfit;
  int iasym, ict, iwood;
  float ctx;
  float woodb1, woodb2;

  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;
  DRMS_RecordSet_t *ffttmprecset = NULL;
  DRMS_Record_t *ffttmprec = NULL;
  DRMS_Segment_t *ffttmpseg = NULL;
  FILE *ffttmpfile = NULL;
  char tstartstr[100];
  char fftfile[DRMS_MAXPATHLEN+1];
  char dirstr[DRMS_MAXPATHLEN+1];

  TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
  double tstart, tstartsave, tstop, tstopsave, tstep;

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

/*#ifdef __linux__*/
/* Used to have to use this:
  int regsz = 4;
After upgrade have to use this: */
#ifdef __ia64
  int regsz = 4;
#else
  int regsz = 1;
#endif

  float *rgaps, *igaps;
  float *gaps, *gaps1, *gaps2, *gaps3;
  float *help;
  float *datr, *dati, *data3;
  float *datarr;

  ffttmp      = (char *)cmdparams_save_str    (&cmdparams, "ffttmp", &newstat);
  logfile     = (char *)cmdparams_save_str    (&cmdparams, "LOGFILE", &newstat);
  readfft     =         cmdparams_save_int    (&cmdparams, "READFFT", &newstat);
//  savefft     =         cmdparams_save_int    (&cmdparams, "SAVEFFT", &newstat);
  gapfile     = (char *)cmdparams_save_str    (&cmdparams, "GAPFILE", &newstat);
  ipar        =         cmdparams_save_int    (&cmdparams, "IPAR", &newstat);
  paramfile   = (char *)cmdparams_save_str    (&cmdparams, "PARAMFILE", &newstat);
  par1file    = (char *)cmdparams_save_str    (&cmdparams, "PAR1FILE", &newstat);
  resfile     = (char *)cmdparams_save_str    (&cmdparams, "RESFILE", &newstat);
  crossfile   = (char *)cmdparams_save_str    (&cmdparams, "CROSSFILE", &newstat);
  ndeltal     =         cmdparams_save_int    (&cmdparams, "NDELTAL", &newstat);
  noisetype   =         cmdparams_save_int    (&cmdparams, "NOISETYPE", &newstat);
  crossnfile  = (char *)cmdparams_save_str    (&cmdparams, "CROSSNFILE", &newstat);
  nnoib       =         cmdparams_save_int    (&cmdparams, "NNOIB", &newstat);
  noibstr     = (char *)cmdparams_save_str    (&cmdparams, "NOIB", &newstat);
  nnoim       =         cmdparams_save_int    (&cmdparams, "NNOIM", &newstat);
  noimstr     = (char *)cmdparams_save_str    (&cmdparams, "NOIM", &newstat);
  csubm       =         cmdparams_save_float  (&cmdparams, "CSUBM", &newstat);
  ifollow     =         cmdparams_save_int    (&cmdparams, "IFOLLOW", &newstat);
  anoise      =         cmdparams_save_float  (&cmdparams, "ANOISE", &newstat);
  aother      =         cmdparams_save_float  (&cmdparams, "AOTHER", &newstat);
//  lmode       = cmdparams_save_int    (&cmdparams, "LMODE", &newstat);
  iodd        =         cmdparams_save_int    (&cmdparams, "IODD", &newstat);
  nmin        =         cmdparams_save_int    (&cmdparams, "NMIN", &newstat);
  nmax        =         cmdparams_save_int    (&cmdparams, "NMAX", &newstat);
  fmin        =         cmdparams_save_float  (&cmdparams, "FMIN", &newstat);
  fmax        =         cmdparams_save_float  (&cmdparams, "FMAX", &newstat);
  imfreq      =         cmdparams_save_int    (&cmdparams, "IMFREQ", &newstat);
  icase       =         cmdparams_save_int    (&cmdparams, "ICASE", &newstat);
  nacoeff     =         cmdparams_save_int    (&cmdparams, "NACOEFF", &newstat);
  ndt         =         cmdparams_save_int    (&cmdparams, "NDT", &newstat);
//  tsample     = cmdparams_save_float  (&cmdparams, "TSAMPLE", &newstat);
  idf         =         cmdparams_save_int    (&cmdparams, "IDF", &newstat);
  ndf0        =         cmdparams_save_int    (&cmdparams, "NDF0", &newstat);
  ndf1        =         cmdparams_save_int    (&cmdparams, "NDF1", &newstat);
  cdf         =         cmdparams_save_float  (&cmdparams, "CDF", &newstat);
  ndm         =         cmdparams_save_int    (&cmdparams, "NDM", &newstat);
  ndl         =         cmdparams_save_int    (&cmdparams, "NDL", &newstat);
  xsafe       =         cmdparams_save_float  (&cmdparams, "XSAFE", &newstat);
  xsafe1      =         cmdparams_save_float  (&cmdparams, "XSAFE1", &newstat);
  delmin      =         cmdparams_save_float  (&cmdparams, "DELMIN", &newstat);
  maxchi      =         cmdparams_save_int    (&cmdparams, "MAXCHI", &newstat);
  maxgrad     =         cmdparams_save_int    (&cmdparams, "MAXGRAD", &newstat);
  maxcof      =         cmdparams_save_int    (&cmdparams, "MAXCOF", &newstat);
  ifix        =         cmdparams_save_int    (&cmdparams, "IFIX", &newstat);
  ifirst      =         cmdparams_save_int    (&cmdparams, "IFIRST", &newstat);
  flipm       =         cmdparams_save_int    (&cmdparams, "FLIPM", &newstat);
  nlobes      =         cmdparams_save_int    (&cmdparams, "NLOBES", &newstat);
  nfour       =         cmdparams_save_int    (&cmdparams, "NFOUR", &newstat);
  gongflag    =         cmdparams_save_int    (&cmdparams, "GONGFLAG", &newstat);
  printfit    =         cmdparams_save_int    (&cmdparams, "PRINTFIT", &newstat);
  cdetrend    =         cmdparams_save_float  (&cmdparams, "CDETREND", &newstat);
  ifill       =         cmdparams_save_int    (&cmdparams, "IFILL", &newstat);
  fdiff       =         cmdparams_save_int    (&cmdparams, "FDIFF", &newstat);
  iasym       =         cmdparams_save_int    (&cmdparams, "IASYM", &newstat);
  ict         =         cmdparams_save_int    (&cmdparams, "ICT", &newstat);
  iwood       =         cmdparams_save_int    (&cmdparams, "IWOOD", &newstat);
  woodb1      =         cmdparams_save_float  (&cmdparams, "WOODB1", &newstat);
  woodb2      =         cmdparams_save_float  (&cmdparams, "WOODB2", &newstat);
  ctx         =         cmdparams_save_int    (&cmdparams, "CTX", &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;
*/

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


  int histflag = strncasecmp("none", outseries, 4);
  if (histflag)
  {
    long long histrecnum;
    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);
    if (histlink != NULL)
    {
      histrecnum=set_history(histlink);
      if (histrecnum < 0)
      {
        fprintf(stderr, "ERROR: problem creating record in history dataseries for output dataseries %s\n", outseries);
        drms_close_record(tempoutrec, DRMS_FREE_RECORD);
        return 1;
      }
      printf("histrecnum=%ld\n",histrecnum);
    }
    else
    {
      fprintf(stderr,"ERROR: could not find history link in output dataseries %s\n", outseries);
      return 1;
    }

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

  }


/* Set PARAMFILE to PAR1FILE if not set */

  if (strcmp(paramfile,"nofile") == 0)
    paramfile=par1file;

  char *crosspath;
  char *dir, *sub, *hold;
  struct stat crosstat;
  if (fstatus=stat(crossfile, &crosstat))
  {
    DRMS_RecordSet_t *crossrecset = drms_open_records(drms_env, crossfile, &status);
    if (status != DRMS_SUCCESS || crossrecset == NULL) 
    {
      fprintf(stderr,"ERROR: problem finding leakage matrix: file status = %d, DRMS status = %d\n",fstatus,status);
      return 1;
    }
    if (crossrecset->n == 0)
    {
      fprintf(stderr,"ERROR: crossfile recordset contains no records\n");
      return 1;
    }
    if (crossrecset->n > 1)
    {
      fprintf(stderr,"ERROR: crossfile recordset with more than 1 record not supported.\n");
      return 1;
    } 
    DRMS_Segment_t *crosseg = drms_segment_lookupnum(crossrecset->records[0], 0);
    crosspath = (char *)malloc(DRMS_MAXPATHLEN+1);
    if (crosseg != NULL && crosspath != NULL)
      status=drms_record_directory(crossrecset->records[0],crosspath,0); 	
    if (crosseg == NULL || status != DRMS_SUCCESS || crosspath[0] == '\0')
    {
      fprintf(stderr, "ERROR: problem finding leakage matrix segment: status = %d\n", status);
      return 1;
    }
    dir=drms_getkey_string(crossrecset->records[0], "dirname", &status);
    sub=strtok(dir,"/");
    while ((hold=strtok(NULL,"/")) != NULL)
      sub=hold;
    strcat(crosspath,"/");
    strcat(crosspath,sub);
    strcat(crosspath,"/");
    if (verbflag)
      printf("crosspath = %s\n", crosspath);
  }
  else
  {
    crosspath=strdup(crossfile);
  }

/* Read window function */
  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_FLOAT, &status);
    if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)
    {
      fprintf(stderr, "ERROR: problem reading gaps segment: status = %d\n", status);
      return 1;
    }
    rgaps=(float *)(gaparr->data);
    gapsize=gaparr->axis[0];
  }
  else
  {
    fits_get_img_size(fitsptr, 1, fdims, &fstatus);
    gapsize=fdims[0];
    rgaps=(float *)(malloc(gapsize*sizeof(float)));
    fits_read_pix(fitsptr, TFLOAT, fpixel, gapsize, NULL, rgaps, 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);

  igaps=(float *)(malloc(gapsize*sizeof(float)));
  gaps=(float *)(malloc(gapsize*sizeof(float)));
  gaps1=(float *)(malloc(gapsize*sizeof(float)));
  gaps2=(float *)(malloc(gapsize*sizeof(float)));
  datr=(float *)(malloc(gapsize*sizeof(float)));
  dati=(float *)(malloc(gapsize*sizeof(float)));
/*
help=(float *)(malloc((4*gapsize+30)*sizeof(float)));
Changed 20061022 */
  help=(float *)(malloc((4*ndt+30)*sizeof(float)));
  gaps3=(float *)(malloc(ndt*sizeof(float)));
  data3=(float *)(malloc(2*ndt*sizeof(float)));

  isign = 1.0;
  if (flipm != 0)
    isign=-1.0;
  sum=0.0;
  for (i=0; i<gapsize; i++)
  {
    sum=sum+rgaps[i];
    igaps[i] = isign*rgaps[i];
    if (rgaps[i] == 0.0)
      gaps[i]=0.0;
    else
      gaps[i]=1.0;
  }
  if (verbflag)
    printf("%f\n",sum);

  cffti_(&ndt, help);


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

  if (inrecset->n > 1)
  {
    fprintf(stderr,"ERROR: nrecs > 1 not yet supported.\n");
    return 1;
  }


  inrec=inrecset->records[0];
  int itest;
  char *inchecklist[] = {"T_START", "T_STOP", "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;
    }
  }

  tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL);
  tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL);
  tstep=drms_getkey_float(inrec, "T_STEP", NULL);

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

  nread=navail;
  tsample=tstep;

  if ((ifirst+ndt) < navail)
    nread=ifirst+ndt;

  if (verbflag)
    printf("%d %d\n",navail,nread);

  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;

  char *tag = drms_getkey_string(inrec, "TAG", &status);
  if (status != DRMS_SUCCESS)
    tag=strdup("none");
  if (readfft == 0)
  {
    ffttmprec = drms_create_record(drms_env, ffttmp, DRMS_PERMANENT, &status);

    if (status != DRMS_SUCCESS) 
    {
      fprintf(stderr,"ERROR: couldn't create a record in ffttmp dataseries %s, status = %d\n", ffttmp, status);
      return 0;
    }

    drms_copykeys(ffttmprec, inrec, 0, kDRMS_KeyClass_Explicit);
    drms_setkey_int(ffttmprec, "LMODE", lmode);
    drms_setkey_int(ffttmprec, "NDT", ndt);
    tnow = (double)time(NULL);
    tnow += UNIX_epoch;
    drms_setkey_time(ffttmprec, "DATE", tnow);

    ffttmpseg = drms_segment_lookupnum(ffttmprec, 0);

    sprintf(fftfile, "fft%d", lmode);
    ffttmpfile =  drms_segment_fopen(ffttmpseg, fftfile, 0, &status);
    if (status != DRMS_SUCCESS) 
    {
      fprintf(stderr,"ERROR: couldn't open a file to write ffttmp record, status = %d\n", status);
      return 0;
    }

    drms_segment_filename(ffttmpseg,fftfile);

  }
  else
  {
    char *ffttmpquery = malloc(100);
    sprint_time(tstartstr, tstart, "TAI", 0); 
    sprintf(ffttmpquery, "%s[%s][%d][%d][%s]", ffttmp,tstartstr,lmode,ndt,tag);
//printf("ffttmpquery = %s\n", ffttmpquery);
    ffttmprecset = drms_open_records(drms_env, ffttmpquery, &status);
    if (status != DRMS_SUCCESS || ffttmprecset == NULL || ffttmprecset->n == 0) 
    {
      fprintf(stderr,"ERROR: couldn't open ffttmp record %s, status = %d\n", ffttmpquery, status);
      return 0;
    }

    ffttmprec = ffttmprecset->records[0];
    ffttmpseg = drms_segment_lookupnum(ffttmprec, 0);

    drms_segment_filename(ffttmpseg,fftfile);

  }

  if (verbflag)
    printf("READFFT=%d,fftfull=%s\n",readfft,fftfile);

//printf("cdetrend=%f,ifill=%d\n",cdetrend,ifill);

/* Loop over m and do FFT's */
  if (readfft == 0)
  {
    if (seginflag)
      segin = drms_segment_lookup(inrec, segnamein);
    else
      segin = drms_segment_lookupnum(inrec, 0);
    if (segin)
      inarr = drms_segment_read(segin, usetype, &status);
    if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL)
    {
      fprintf(stderr, "problem reading input segment: status = %d\n", status);
      return 0;
    }

    datarr=(float *)(inarr->data);
 
    for (m=0; m<= lmode; m++)
    {

//      printf("%d %d\n",ndt,m);

      data=datarr + m*2*gapsize;

      for (j=0; j<navail; j++)
      {
        if (isnan(data[2*j])) 
          datr[j]=0.0;
        else
          datr[j]=rgaps[j]*data[2*j];
        if (isnan(data[2*j+1])) 
          dati[j]=0.0;
        else
          dati[j]=igaps[j]*data[2*j+1];
      }

      if (cdetrend != 0.0)
      {
        detrend(datr,gaps,navail,cdetrend);
        detrend(dati,gaps,navail,cdetrend);
      }

      if (ifill == 0)
      {
        for (i=0; i<navail; i++)
        {
          gaps1[i]=gaps[i];
        }
      }
      else
      {
        fill_gaps(datr,gaps,gaps1,navail);
        if (m !=0 )
        {
          fill_gaps(dati,gaps,gaps1,navail);
        }
      }

      if (fdiff!=0)
      {
        for (j=0; j<(navail-1); j++)
        {
          if ((gaps1[j]==1) & (gaps1[j+1]==1))
          {
            gaps2[j]=1.0;
            datr[j]=datr[j+1]-datr[j];
            dati[j]=dati[j+1]-dati[j];
          }
          else
          {
            gaps2[j]=0.0;
            datr[j]=0.0;
            dati[j]=0.0;
          }
        }
        gaps2[navail-1]=0.0;
        datr[navail-1]=0.0;
        dati[navail-1]=0.0;
      }
      else
      {
        for (j=0; j<navail; j++)
        {
          gaps2[j]=gaps1[j];
        }
      }

/* Make complex array.
Now use data3 instead of reusing data in order to allow ndt>navail. */
      nx1=ndt;
      if (nread < (ifirst+ndt))
        nx1=nread-ifirst;

      for (j=0; j<nx1; j++)
      {
        gaps3[j]=gaps2[j+ifirst];
        data3[2*j]=datr[j+ifirst];
        data3[2*j+1]=dati[j+ifirst];
      }
      for (j = nx1; j < ndt; j++)
      {
        gaps3[j] = 0;
        data3[2*j] = 0;
        data3[2*j+1] = 0;
      }

      if ((ifix == 1) && (m == 0))
      {
        c = sqrt(2.0);
        for (j = 0; j < ndt; j++)
        {
          data3[2*j] *= c;
          data3[2*j+1] *= c;
        }
      }

      cfftb_ (&ndt, data3, help);

      fwrite ( (char*)data3, 8*ndt, 1, ffttmpfile);


    } /* End of m loop */

  /* Pad the file for buffered reading in Fortran when NDT is small */
  fwrite ( (char*)data3, 8*ndt, 1, ffttmpfile);
  drms_segment_fclose(ffttmpfile);
  if (verbflag)
    printf("fft file written\n");
  }
  else
  {
/* Fill dummy time-series to make gaps */
    for (j=0; j<navail; j++)
    {
      datr[j]=j;
    }
    if (cdetrend != 0.0)
    {
      detrend(datr,gaps,navail,cdetrend);
    }
    if (ifill == 0)
    {
      for (i=0; i<navail; i++)
      {
        gaps1[i]=gaps[i];
      }
    }
    else
    {
      fill_gaps(datr,gaps,gaps1,navail);
    }
    if (fdiff!=0)
    {
      for (j=0; j<(navail-1); j++)
      {
        if ((gaps1[j]==1) & (gaps1[j+1]==1))
        {
          gaps2[j]=1.0;
        }
        else
        {
          gaps2[j]=0.0;
        }
      }
      gaps2[navail-1]=0.0;
    }
    else
    {
      for (j=0; j<navail; j++)
      {
        gaps2[j]=gaps1[j];
      }
    }
    nx1=ndt;
    if (nread < (ifirst+ndt))
      nx1=nread-ifirst;
    for (j=0; j<nx1; j++)
    {
      gaps3[j]=gaps2[j+ifirst];
    }
    for (j = nx1; j < ndt; j++)
    {
      gaps3[j] = 0;
    }
  }

  drms_free_array(inarr);
  if (inrecset)
  {
    drms_close_records(inrecset, DRMS_FREE_RECORD);
  }
  if (verbflag)
    printf("input recordset closed\n");

  i = -1;
  for ( ptr = noibstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) 
    noib[++i] = atoi(ptr);

  i = -1;
  for ( ptr = noimstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) 
    noim[++i] = atoi(ptr);

  loglen = strlen(logfile);
  fftlen = strlen(fftfile);
  paramlen = strlen(paramfile);
  par1len = strlen(par1file);
  reslen = strlen(resfile);
  crosslen = strlen(crosspath);
  crossnlen = strlen(crossnfile);
printf("crosspath=%s, crosslen=%d\n",crosspath,crosslen);

/* Call Fortran program */
  sub24_(logfile,fftfile, gaps3, &ipar,
       paramfile, par1file, resfile,
       crosspath, &ndeltal, &noisetype, crossnfile, &nnoib, noib, &nnoim, noim,
       &csubm, &ifollow, &anoise, &aother, &lmode, &iodd, &iasym, &ict, &ctx,
       &iwood, &woodb1, &woodb2, &nmin, &nmax, &fmin, &fmax,
       &imfreq, &icase, &nacoeff, &zero, &ndt, &tsample, &idf, &ndf0, &ndf1,
       &cdf, &ndm, &ndl,
       &xsafe, &xsafe1, &delmin, &maxchi, &maxgrad, &maxcof, &regsz,
       &nlobes, &nfour,
       &gongflag, &printfit,
       loglen, fftlen, paramlen, par1len, reslen, crosslen, crossnlen);

  if (readfft == 0)
    drms_close_record(ffttmprec,DRMS_INSERT_RECORD);
  else
    drms_close_records(ffttmprecset,DRMS_FREE_RECORD);

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


void detrend(
float *data,
float *gaps,
int length,
float cdetrend
)
{
#define chunksz 1440
#define maxpols 50

extern float sdot_(int *, float *, int *, float *, int *);
extern float saxpy_(int *, float *, float *, int *, float *, int *);
extern float sscal_(int *, float *, float *, int *);

int i,j,k,l,n_good,ndegree;
int goodlist[chunksz];
float pols[maxpols][chunksz];
float a,cg,rms2,sum,x[chunksz];
float *help;
int one=1;

for (i=0; i<length ; i=i+chunksz) {
  help=data+i;
  n_good=0;
  cg=0.0;
  for (j=i; j<i+chunksz; j++) {
    if (gaps[j] != 0) {
      goodlist[n_good]=j-i;
      cg=cg+j-i;
      n_good++;
    }
  }
  if (n_good != 0) {
    cg=cg/n_good;
    rms2=0.0;
    for (j=0; j<n_good; j++) {rms2=rms2+(goodlist[j]-cg)*(goodlist[j]-cg);}
    ndegree = floor((goodlist[n_good-1]-goodlist[0])/cdetrend)+2;
    for (k=0; k<n_good; k++) {
      x[k]=(goodlist[k]-cg)/sqrt(rms2);
    }
    for (k=0; k<n_good; k++) {pols[0][k]=1.0/sqrt(n_good);}
    for (k=0; k<n_good; k++) {pols[1][k]=x[k];}
    for (j=2; j<ndegree; j++) {
      for (k=0; k<n_good; k++) {
        pols[j][k]=(2.0*j)/j*x[k]*pols[j-1][k]-(j-1.0)/j*pols[j-2][k];
      }
      for (l=0; l<j; l++) {
        a=-sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one);
        saxpy_(&n_good,&a,&pols[l][0],&one,&pols[j][0],&one);
      }
      a=1.0/sqrt(sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one));
      sscal_(&n_good,&a,&pols[j][0],&one);
    }
    for (j=0; j<ndegree; j++) {
      sum=0.0;
      for (k=0; k<n_good; k++) {
        sum=sum+pols[j][k]*help[goodlist[k]];
      }
      for (k=0; k<n_good; k++) {
        help[goodlist[k]] -= sum*pols[j][k];
      }
    }
  } /* if n_good */

} /* for i */

}

void memcof(float *, int, int, float *, float *);
void fixrts(float *, int);
void predic(float *, int, float *, int, float *, int);

void fill_gaps(
float *data,
float *gaps,
float *newgaps,
int length
)
{
#include "nr.h"
#include "nrutil.h"
#define chunksz 1440
#define maxpoles 10
#define maxgap 5

int i,j,k;
int n_good, gap_start, gap_end, gap_size, type;
float *help,good_points[chunksz];
int goodlist[chunksz];
int npoles;
float pm,cof[maxpoles];
float input[maxpoles];
float f_predic[maxgap],b_predic[maxgap];

npoles=6;
for (i=0; i<length ; i=i+chunksz) {
  help=data+i;
  n_good=0;
  for (j=i; j<i+chunksz; j++) {
    if (gaps[j] != 0) {
      goodlist[n_good]=j-i;
      good_points[n_good]=data[j];
      n_good++;
      newgaps[j]=gaps[j];
    }
  }
  if (n_good != 0) {
    memcof(good_points-1,n_good,npoles,&pm,cof-1);
    fixrts(cof-1,npoles);
    for (j=1; j<n_good; j++) {
      if ((goodlist[j]-goodlist[j-1]) > 1) {
        gap_start=goodlist[j-1]+1;
        gap_end=goodlist[j]-1;
        gap_size=gap_end-gap_start+1;
        if (gap_size <= maxgap) {
          type=0;
/* Forward prediction */
/* Check that there are npoles good points before gap */
          if ((j >= npoles) && (goodlist[j-npoles] == goodlist[j-1]-npoles+1)) {
            for (k=0; k<npoles; k++) input[k]=help[gap_start-npoles+k];
            predic(input-1,npoles,cof-1,npoles,f_predic-1,gap_size);
            type +=1;
          }
/* Backwards prediction */
/* Check that there are npoles good points after gap */
          if ((j+npoles <= n_good) && (goodlist[j+npoles-1] == goodlist[j]+npoles-1)) {
            for (k=0; k<npoles; k++) input[k]=help[gap_end+npoles-k];
/*          for (k=0; k<npoles; k++) printf("%d %d %f\n",gap_end,k,help[gap_end+npoles-k]); */
            predic(input-1,npoles,cof-1,npoles,b_predic-1,gap_size);
            type +=2;
          }
          if (type == 3) {
/* Average predictions */
            for (k=0; k<gap_size; k++)
              help[gap_start+k] = (f_predic[k]+b_predic[gap_size-k-1])/2.;
          }
          if (type == 1) {
/* Use forward prediction */
            for (k=0; k<gap_size; k++)
              help[gap_start+k] = f_predic[k];
          }
          if (type == 2) {
/* Use backwards prediction */
            for (k=0; k<gap_size; k++)
              help[gap_start+k] = b_predic[gap_size-k-1];
          }
          if (type !=0 ) {
            for (k=0; k<gap_size; k++)
              newgaps[gap_start+k+i]=1;
          }
        }
      }
    } /* for j */
  } /* if n_good */
} /* for i */

}