JSOC/proj/globalhs/apps/jtsfiddle.c - diff

Return to jtsfiddle.c CVS log

Up to [Development] / JSOC / proj / globalhs / apps

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

version 1.2, 2010/04/23 05:47:51

version 1.3, 2011/08/16 02:29:10

Line 1

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

Line 33 char *module_name = "jtsfiddle";

Line 34 char *module_name = "jtsfiddle";

ModuleArgs_t module_args[] =

{

{ARG_STRING, "in", "", "input data records"},

{ARG_STRING, "out", "", "output data series"},

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

Line 48 ModuleArgs_t module_args[] =

Line 55 ModuleArgs_t module_args[] =

// {ARG_FLOAT, "TSAMPLE", "60", "", ""},

{ARG_INT, "IFILL", "1", "", ""},

{ARG_INT, "IFIRST", "0", "", ""},

{ARG_INT, "FLIPM", "1", "", ""},

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

{ARG_INT, "MAX_AR_ORDER", "360", "", ""},

{ARG_INT, "FU_FIT_EDGES", "1", "", ""},

{ARG_INT, "FU_ITER", "3", "", ""},

Line 56 ModuleArgs_t module_args[] =

Line 63 ModuleArgs_t module_args[] =

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

Line 72 ModuleArgs_t module_args[] =

Line 82 ModuleArgs_t module_args[] =

#include "saveparm.c"

#include "timing.c"

#include "set_history.c"

#define DIE(code) { fprintf(stderr,"jtsfiddle died with error code %d\n",(code)); return 0;}

Line 80 static char *logfile, *gapfile, *section

Line 91 static char *logfile, *gapfile, *section

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;

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;

Line 102 static int read_section_file(char *filen

Line 113 static int read_section_file(char *filen

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;

Line 115 int DoIt(void)

Line 125 int DoIt(void)

fftwf_plan plan;

int num_sections, *sections;

float *ar_coeff=NULL;

int ndt;

float *datarr, *datptr, *outptr;

float *datahold;

char tstartstr[100];

int lmin, lmax;

int fstatus = 0;

fitsfile *fitsptr;

Line 122 int DoIt(void)

Line 137 int DoIt(void)

long fpixel[]={1};

int *anynul=0;

int length[2];

int instart[2], inend[2];

float *datarr;

int tslength[2], tsstart[2], tsend[2];

int fft1length[2], fft1start[2], fft1end[2];

int powlength[2], powstart[2], powend[2];

int status=DRMS_SUCCESS;

int newstat=0;

Line 156 int DoIt(void)

Line 173 int DoIt(void)

ct0=getcputime(&ut0, &st0);

/* Read input parameters. */

logfile = cmdparams_save_str (&cmdparams, "LOGFILE", &newstat);

logfile = (char *)cmdparams_save_str (&cmdparams, "LOGFILE", &newstat);

gapfile = cmdparams_save_str (&cmdparams, "GAPFILE", &newstat);

gapfile = (char *)cmdparams_save_str (&cmdparams, "GAPFILE", &newstat);

sectionfile = cmdparams_save_str (&cmdparams, "SECTIONFILE", &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);

Line 167 int DoIt(void)

Line 184 int DoIt(void)

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

// 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 = cmdparams_save_str(&cmdparams, "in", &newstat);

char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat);

char *outSeries = cmdparams_save_str(&cmdparams, "out", &newstat);

// char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);

char *segnamein = cmdparams_save_str(&cmdparams, "segin", &newstat);

char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat);

char *segnameout = cmdparams_save_str(&cmdparams, "segout", &newstat);

char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat);

int seginflag = strcmp(kNOTSPECIFIED, segnamein);

int segoutflag = strcmp(kNOTSPECIFIED, segnameout);

char *histlinkname = cmdparams_save_str(&cmdparams, "histlink", &newstat);

char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);

char *srclinkname = cmdparams_save_str(&cmdparams, "srclink", &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)

Line 198 int DoIt(void)

Line 218 int DoIt(void)

else

lifetime = DRMS_TRANSIENT;

#include "histinclude.c"

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_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. ****/

/* Default is to just output the timeseries. */

if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0))

tsout=1;

/* Only one type of output allowed at a time */

if ((tsout+fftout+fft1out+powout+mavgout+arpowout) !=1)

{

fprintf(stderr, "ERROR: only one type of output allowed at a time\n");

fprintf(stderr, "ERROR: must specify an output dataseries\n");

return 0;

return 1;

}

if (navg <= 0)

navg=1;

if (nskip <= 0)

Line 218 int DoIt(void)

Line 261 int DoIt(void)

if ((navg != 1) && (nskip != 1))

{

fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n");

return 0;

return 1;

}

if (noff < 0 || noff >= nskip)

{

fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");

return 0;

return 1;

}

if(arpowout && !ifill)

{

fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");

return 0;

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

int is, itest, mflip;

char *holdseries="";

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

{

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[is-1];

}

else

{

fprintf(stderr,"WARNING: could not find history link in output dataseries %s\n", serieslist[is]);

}

int itest;

for (itest=0; itest < ARRLENGTH(outchecklist); itest++)

{

DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);

if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)

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

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

return 1;

}

drms_close_record(tempoutrec, DRMS_FREE_RECORD);

mflip = drms_getkey_int(tempoutrec, "MFLIPPED", &status);

if (status == DRMS_SUCCESS)

mfliparr[is]=mflip;

// xmem_off();

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

{

Line 258 int DoIt(void)

Line 372 int DoIt(void)

if (status != DRMS_SUCCESS || gaprecset == NULL)

{

fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);

return 0;

return 1;

}

if (gaprecset->n == 0)

{

fprintf(stderr,"ERROR: gapfile recordset contains no records\n");

return 0;

return 1;

}

if (gaprecset->n > 1)

{

fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");

return 0;

return 1;

}

gapseg = drms_segment_lookupnum(gaprecset->records[0], 0);

if (gapseg != NULL)

Line 276 int DoIt(void)

Line 390 int DoIt(void)

if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL)

{

fprintf(stderr, "problem reading gaps segment: status = %d\n", status);

return 0;

return 1;

}

agaps=(int *)(gaparr->data);

gapsize=gaparr->axis[0];

Line 291 int DoIt(void)

Line 405 int DoIt(void)

if (fstatus)

{

fits_report_error(stderr, fstatus);

return 0;

return 1;

}

Line 301 int DoIt(void)

Line 415 int DoIt(void)

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

{

Line 320 int DoIt(void)

Line 443 int DoIt(void)

if (status != DRMS_SUCCESS || secrecset == NULL)

{

fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status);

return 0;

return 1;

}

if (secrecset->n == 0)

{

fprintf(stderr,"ERROR: sectionfile recordset contains no records\n");

return 0;

return 1;

}

if (secrecset->n > 1)

{

fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n");

return 0;

return 1;

}

num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL);

if (num_sections < 1)

{

fprintf(stderr,"ERROR: Invalid NSECS keywords\n");

return 0;

return 1;

}

sections = (int *)malloc(2*(num_sections)*sizeof(int));

char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL);

Line 351 int DoIt(void)

Line 474 int DoIt(void)

if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1]))

{

fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n");

return 0;

return 1;

}

free(sectkey);

Line 360 int DoIt(void)

Line 483 int DoIt(void)

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

/* allocate temporary storage */

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

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

DRMS_RecordSet_t *inRecSet = drms_open_records(drms_env, inRecQuery, &status);

if (status != DRMS_SUCCESS || inRecSet == NULL)

if (status != DRMS_SUCCESS || inrecset == NULL)

{

fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);

return 0;

return 1;

}

if (inrecset->n == 0)

if (inRecSet->n == 0)

{

printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]);

printf("ERROR: input recordset contains no records\n");

return 0;

return 1;

}

if (inrecset->n > 1)

if (inRecSet->n > 1)

{

printf("ERROR: nRecs > 1 not yet supported.\n");

printf("ERROR: nrecs > 1 not yet supported.\n");

return 0;

return 1;

}

inrec=inrecset->records[0];

inrec=inRecSet->records[0];

char *inchecklist[] = {"T_START", "T_STOP", "QUALITY", "LMIN", "LMAX", "T_STEP"};

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)

if (inkeytest == NULL)

{

fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);

drms_close_records(inRecSet, DRMS_FREE_RECORD);

drms_close_records(inrecset, DRMS_FREE_RECORD);

return 0;

return 1;

}

ndt=drms_getkey_int(inrec, "NDT", NULL);

tstep=drms_getkey_float(inrec, "T_STEP", NULL); //assume constant across all input records

tsample=tstep;

// insert loop over input here

int mflipin=drms_getkey_int(inrec, "MFLIPPED", &status);

int iRec=0;

if (status != DRMS_SUCCESS) mflipin=0;

inrec = inRecSet->records[iRec];

int lmin=drms_getkey_int(inrec, "LMIN", NULL);

int lmax=drms_getkey_int(inrec, "LMAX", NULL);

if (lmin != lmax)

{

printf("ERROR: lmin != lmax not yet supported.\n");

return 0;

}

lmode=lmin;

tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL);

// already required mfliparr[TSOUT] == mfliparr[FFTOUT] above

tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL);

if (tsout && mfliparr[TSOUT] != mflipin)

tstep=drms_getkey_float(inrec, "T_STEP", NULL);

flipm=1;

else if (fftout && mfliparr[FFTOUT] != mflipin)

flipm=1;

else

flipm=0;

n_in=(tstop-tstart)/tstep;

for (is=2;is<5;is++)

if (n_in != gapsize)

{

fprintf(stderr, "ERROR: gaps seem incompatible with time-series, gapsize=%d, n_in=%d\n", gapsize, n_in);

if (!flagarr[is])

return 0;

continue;

if (mfliparr[is] != (mflipin || flipm))

msignarr[is]=-1;

else

msignarr[is]=1;

}

tsample=tstep;

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

return 1;

}

/* Adjust detrend parameters according to the number of

Line 448 int iRec=0;

Line 559 int iRec=0;

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;

Line 460 int iRec=0;

Line 579 int iRec=0;

if (mavgout)

msum = (float *)calloc(n_sampled_out/2,sizeof(float));

if (fft1out || powout || mavgout || arpowout)

tstartout=tstart+ifirst*tsample;

tstopout=tstartout+df1;

/* Create output data set. */

if (tsout || fftout)

{

length[0]=2*n_sampled_out;

length[1]=lmode+1;

}

else

{

/* Set first and last output index if not set as a parameter. */

if (imax < imin)

Line 481 int iRec=0;

Line 589 int iRec=0;

}

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]=2*n_sampled_out;

tslength[1]=1;

tsstart[0]=0;

tsend[0]=2*n_sampled_out-1;

}

if (fft1out)

length[0]=2*npow;

{

else

fft1length[0]=2*npow;

length[0]=npow;

fft1length[1]=1;

if (powout || arpowout || fft1out)

fft1start[0]=0;

length[1]=2*lmode+1;

fft1end[0]=2*npow-1;

else

}

length[1]=1;

if (powout || arpowout || mavgout)

{

powlength[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;

}

/* Loops through all the data and update the gaps to reflect missing

data. */

/* the above comment is from the original code. it refers to an ifdef

that would read in all the data and run extract_data on it in case there

were nan's not removed in the input window function. the present version

assumes the window function is correct. */

// insert loop over input here

int irec=0;

inrec = inrecset->records[irec];

lmin=drms_getkey_int(inrec, "LMIN", NULL);

lmax=drms_getkey_int(inrec, "LMAX", NULL);

if (lmin != lmax)

{

printf("ERROR: lmin != lmax not yet supported.\n");

return 0;

}

lmode=lmin;

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

n_in=(tstop-tstart)/tstep;

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)

if (segin == NULL)

inarr = drms_segment_read(segin, usetype, &status);

if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL)

{

fprintf(stderr, "ERROR: problem reading input segment: iRec = %d, status = %d\n", iRec, status);

fprintf(stderr, "ERROR: problem looking up input segment: recnum = %lld\n", inrec->recnum);

return 0;

}

datarr=(float *)(inarr->data);

for (is=0; is<6; is++)

{

if (!flagarr[is])

continue;

outrec = drms_create_record(drms_env,

outrecarr[is] = drms_create_record(drms_env,

outSeries,

serieslist[is],

lifetime,

DRMS_PERMANENT,

&status);

if (status != DRMS_SUCCESS || outrec==NULL)

if (status != DRMS_SUCCESS)

{

fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status);

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

return 0;

}

if (histlink)

if (histrecnumarr[is] > 0)

drms_setlink_static(outrec, histlinkname, histrecnum);

drms_setlink_static(outrecarr[is], histlinkname, histrecnumarr[is]);

if (srclink)

drms_setlink_static(outrecarr[is], srclinkname, inrec->recnum);

drms_setlink_static(outrec, srclinkname, inrec->recnum);

if (segoutflag)

segout = drms_segment_lookup(outrec, segnameout);

segoutarr[is] = drms_segment_lookup(outrecarr[is], segnameout);

else

segout = drms_segment_lookupnum(outrec, 0);

segoutarr[is] = drms_segment_lookupnum(outrecarr[is], 0);

outarr = drms_array_create(usetype, 2, length, NULL, &status);

if (segoutarr[is] == NULL)

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

fprintf(stderr, "ERROR: problem looking up outputput segment in series %s\n", serieslist[is]);

return 0;

}

out_data = outarr->data;

// 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:

;

}

/***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/

for (m=0; m<=lmode; m++)

{

// (*history)("### processing m = %d ###\n",m);

/* Read chunk of data corresponding to next m from the data file. */

if (verbflag > 1)

start_index[1] = m;

if (!(tile_sds = vds_select(in_vds,0,0,1,1,start_index,counts,intervals)))

{

fprintf(stderr,"vds_select failed");

DIE(9);

}

if (!sds_data(tile_sds))

{

fprintf(stderr,"*** lnu vds_select returned no data.\n");

fprintf(stdout, "starting m=%d\n", m);

sds_free(&tile_sds);

DIE(10);

}

in_data = (float *)(sds_data(tile_sds));

if (sds_dim0(tile_sds)<0)

instart[1]=m;

inend[1]=m;

inarr = drms_segment_readslice(segin, usetype, instart, inend, &status);

if (status != DRMS_SUCCESS || inarr == NULL )

{

sds_free(&tile_sds);

fprintf(stderr, "ERROR: problem reading input segment: status = %d, recnum = %lld\n", status, inrec->recnum);

DIE(soi_errno);

return 0;

}

in_data=(float *)(inarr->data);

in_data=datarr + m*2*gapsize;

/* Extracts data copy */

extract_data(n_sampled_in, gaps, in_data, data);

Line 578 int iRec=0;

Line 754 int iRec=0;

if (detrend_const != 0)

{

detrend_order = floor(detrend_length/detrend_const)+2;

// for (i=0;i<num_sections; i++)

// printf("[%d:%d]\n",sections[2*i],sections[2*i+1]);

cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps,

detrend_order, detrend_length, detrend_skip,

num_sections,sections);

num_sections, sections, detrend_first);

}

/* Fill gaps in entire time series if required */

memcpy(gaps2, gaps, n_sampled_in*sizeof(int));

if (ifill != 0 && max_ar_order > 0)

{

// (*history)("Filling gaps, ");

cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2,

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)

{

FILE *fh;

for (i = 0; i < n_sampled_out; i++)

unsigned char c;

fh = fopen("gaps.bin","w");

for(i=0; i<n_sampled_in; i++)

{

c = (unsigned char) gaps2[i];

datatemp[2*i]=data[2*i];

fwrite(&c,1,1,fh);

datatemp[2*i+1]=-data[2*i+1];

}

fclose(fh);

tmpptr=datatemp;

}

else

tmpptr=data;

// (*history)("AR model order used = %d.\n",ar_order);

tsstart[1]=m;

tsend[1]=m;

outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);

status=drms_segment_writeslice(segoutarr[TSOUT], outarr, tsstart, tsend, 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);

/* pad with zeros if the last point output (n_sampled_out+ifirst)

if (fftout)

is after the last data points read in (nread). */

if (arpowout)

{

memmove(data, ar_coeff, 2*(ar_order+1)*sizeof(float));

// code to flip m on this output separately. in the present code this has already been accomplished by setting flipm.

memset(&data[2*(ar_order+1)],0,2*(n_sampled_out-ar_order-1)*sizeof(float));

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

status=drms_segment_writeslice(segoutarr[FFTOUT], outarr, tsstart, tsend, 0);

free(outarr);

if (status != DRMS_SUCCESS)

{

memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float));

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

if ((ifirst+n_sampled_out) >= n_sampled_in)

return 0;

memset(&data[2*n_sampled_in], 0,

2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float));

}

if (fftout || fft1out || powout || arpowout || mavgout)

{

// (*history)("Computing FFT.\n");

fftwf_execute(plan);

}

if (fft1out)

{

// (*history)("Preparing output.\n");

fft1start[1]=lmode+m*msignarr[FFT1OUT];

fft1end[1]=lmode+m*msignarr[FFT1OUT];

if (tsout || fftout)

outarr = drms_array_create(usetype, 2, fft1length, data+2*imin, &status);

// memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float));

status=drms_segment_writeslice(segoutarr[FFT1OUT], outarr, fft1start, fft1end, 0);

// the above line is above, already moved shifted ifirst to start of data, call should be this:

free(outarr);

memcpy(&out_data[2*m*n_sampled_out], data, 2*n_sampled_out*sizeof(float));

if (status != DRMS_SUCCESS)

else if (fft1out)

{

/* Do positive m */

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

memcpy(&out_data[2*(lmode+m)*npow], &data[2*imin], 2*npow*sizeof(float));

return 0;

}

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

datahold[2*i] = data[2*(n_sampled_out-imin-i)];

out_data[2*((lmode-m)*npow+i)+1] = -data[2*(n_sampled_out-imin-i)+1];

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

status=drms_segment_writeslice(segoutarr[FFT1OUT], outarr, fft1start, fft1end, 0);

free(outarr);

if (status != DRMS_SUCCESS)

{

return 0;

}

else

{

/* Compute power spectrum from complex FFT. */

if (arpowout)

{

for (i = 0; i < n_sampled_out; i++)

pow[i] = 1.0/(data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]);

}

else

{

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

if (powout)

{

/* Do positive m */

powstart[1]=lmode+m*msignarr[POWOUT];

memcpy(&out_data[(lmode+m)*npow], &pow[imin], npow*sizeof(float));

powend[1]=lmode+m*msignarr[POWOUT];

outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status);

status=drms_segment_writeslice(segoutarr[POWOUT], outarr, powstart, powend, 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)

out_data[(lmode-m)*npow] = pow[n_sampled_out-imin];

datahold[0] = pow[n_sampled_out-imin];

else

out_data[(lmode-m)*npow] = pow[0];

datahold[0] = pow[0];

for (i=1; i<npow;i++)

out_data[(lmode-m)*npow+i] = pow[n_sampled_out-imin-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);

status=drms_segment_writeslice(segoutarr[POWOUT], outarr, powstart, powend, 0);

free(outarr);

if (status != DRMS_SUCCESS)

{

return 0;

}

else

}

if (mavgout)

{

/* m-averaged output */

mshift = floor(df1*splitting(lmode,m*msignarr[MAVGOUT])+0.5f);

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

if (m > 0)

{

/* Do negative m */

mshift = floor(df1*splitting(lmode,-m)+0.5f);

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

Line 706 int iRec=0;

Line 932 int iRec=0;

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

status=drms_segment_writeslice(segoutarr[ARPOWOUT], outarr, powstart, powend, 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 m loop */

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

hiter_new(&outKey_list, &outrec->keywords);

powend[1]=lmode-m*msignarr[ARPOWOUT];

while ( (outKey=(DRMS_Keyword_t *)hiter_getnext(&outKey_list)) )

outarr = drms_array_create(usetype, 2, powlength, datahold, &status);

status=drms_segment_writeslice(segoutarr[ARPOWOUT], outarr, powstart, powend, 0);

free(outarr);

if (status != DRMS_SUCCESS)

{

char *keyName = outKey->info->name;

DRMS_Value_t inValue = {DRMS_TYPE_STRING, NULL};

return 0;

inValue = drms_getkey_p(inrec, keyName, &status);

if (status) fprintf(stderr,"*** COPY drms_getkey_p %s status=%d\n",keyName,status);

drms_setkey_p(outrec, keyName, &inValue);

if ((inValue.type == DRMS_TYPE_STRING) && inValue.value.string_val)

free(inValue.value.string_val);

inValue.value.string_val = NULL;

}

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

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=outrecarr[is];

drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit);

drms_setkey_int(outrec, "QUALITY", 0);

drms_setkey_time(outrec, "T_START", tstartout);

Line 732 int iRec=0;

Line 1006 int iRec=0;

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

drms_free_array(inarr);

drms_close_record(outrec, DRMS_INSERT_RECORD);

if (inRecSet)

{

drms_close_records(inRecSet, DRMS_FREE_RECORD);

}

if (mavgout)

if (inrecset != NULL)

{

c=1.0f/(2*lmode+1);

drms_close_records(inrecset, DRMS_FREE_RECORD);

for (i=0; i<npow; i++)

out_data[i] = msum[imin+i] * c;

}

if(ar_coeff)

free(ar_coeff);

drms_segment_write(segout, outarr, 0);

drms_close_record(outrec, DRMS_INSERT_RECORD);

if (fftout || fft1out || powout || arpowout || mavgout)

fftwf_destroy_plan(plan);

Line 767 int iRec=0;

Line 1033 int iRec=0;

// printf("number of records created = %d\n", nsuccess);

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

}

// xmem_leakreport();

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

return 0;

}

Line 821 static double splitting(int l, int m)

Line 1086 static double splitting(int l, int m)

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

{

Line 922 static int extract_gaps(int n_in, int *g

Line 1189 static int extract_gaps(int n_in, int *g

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;

Legend:

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

Karen Tian