version 1.1, 2010/01/29 21:13:02
|
version 1.2, 2010/04/23 05:47:51
|
Line 41 ModuleArgs_t module_args[] = |
|
Line 41 ModuleArgs_t module_args[] = |
|
{ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"}, | {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"}, |
{ARG_INT, "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL}, | {ARG_INT, "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL}, |
{ARG_STRING, "LOGFILE", "stdout", ""}, | {ARG_STRING, "LOGFILE", "stdout", ""}, |
{ARG_STRING, "GAPFILE", "", ""}, |
{ARG_STRING, "GAPFILE", "", "record or file containing window function"}, |
{ARG_STRING, "SECTIONFILE", "none", ""}, |
{ARG_STRING, "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"}, |
// {ARG_INT, "LMODE", "", "", ""}, |
|
{ARG_INT, "NDT", "-1", "", ""}, | {ARG_INT, "NDT", "-1", "", ""}, |
{ARG_FLOAT, "TSAMPLE", "60", "", ""}, |
// {ARG_FLOAT, "TSAMPLE", "60", "", ""}, |
{ARG_INT, "IFILL", "1", "", ""}, | {ARG_INT, "IFILL", "1", "", ""}, |
{ARG_INT, "IFIRST", "0", "", ""}, | {ARG_INT, "IFIRST", "0", "", ""}, |
{ARG_INT, "FLIPM", "1", "", ""}, | {ARG_INT, "FLIPM", "1", "", ""}, |
|
|
DRMS_Type_t usetype = DRMS_TYPE_FLOAT; | DRMS_Type_t usetype = DRMS_TYPE_FLOAT; |
DRMS_RecLifetime_t lifetime; | DRMS_RecLifetime_t lifetime; |
long long histrecnum = -1; | long long histrecnum = -1; |
|
DRMS_Segment_t *gapseg = NULL; |
|
DRMS_Array_t *gaparr = NULL; |
| |
HIterator_t outKey_list; | HIterator_t outKey_list; |
DRMS_Keyword_t *outKey; | DRMS_Keyword_t *outKey; |
|
|
int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); | int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); |
int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); | int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); |
| |
double tstart, tstartsave, tstop, tstopsave, tstep; |
double tstart, tstartsave, tstop, tstopsave, tstep, tstartout, tstopout; |
double wt0, wt1, wt2, wt3, wt; | double wt0, wt1, wt2, wt3, wt; |
double ut0, ut1, ut2, ut3, ut; | double ut0, ut1, ut2, ut3, ut; |
double st0, st1, st2, st3, st; | double st0, st1, st2, st3, st; |
|
|
logfile = cmdparams_save_str (&cmdparams, "LOGFILE", &newstat); | logfile = cmdparams_save_str (&cmdparams, "LOGFILE", &newstat); |
gapfile = cmdparams_save_str (&cmdparams, "GAPFILE", &newstat); | gapfile = cmdparams_save_str (&cmdparams, "GAPFILE", &newstat); |
sectionfile = cmdparams_save_str (&cmdparams, "SECTIONFILE", &newstat); | sectionfile = cmdparams_save_str (&cmdparams, "SECTIONFILE", &newstat); |
// lmode = cmdparams_save_int (&cmdparams, "LMODE", &newstat); |
|
n_sampled_out = cmdparams_save_int (&cmdparams, "NDT", &newstat); | n_sampled_out = cmdparams_save_int (&cmdparams, "NDT", &newstat); |
tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat); |
// tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat); |
ifill = cmdparams_save_int (&cmdparams, "IFILL", &newstat); | ifill = cmdparams_save_int (&cmdparams, "IFILL", &newstat); |
max_ar_order= cmdparams_save_int (&cmdparams, "MAX_AR_ORDER", &newstat); | max_ar_order= cmdparams_save_int (&cmdparams, "MAX_AR_ORDER", &newstat); |
fu_fit_edges= cmdparams_save_int (&cmdparams, "FU_FIT_EDGES", &newstat); | fu_fit_edges= cmdparams_save_int (&cmdparams, "FU_FIT_EDGES", &newstat); |
|
|
return 0; | return 0; |
} | } |
| |
// insert checks on output keywords here |
|
|
// these must be present in the output dataseries and variable, not links or constants |
|
char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX"}; |
|
|
|
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) |
|
{ |
|
fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]); |
|
drms_close_record(tempoutrec, DRMS_FREE_RECORD); |
|
return 0; |
|
} |
|
} |
|
|
drms_close_record(tempoutrec, DRMS_FREE_RECORD); | drms_close_record(tempoutrec, DRMS_FREE_RECORD); |
| |
// xmem_off(); | // xmem_off(); |
| |
| |
DRMS_RecordSet_t *inRecSet = drms_open_records(drms_env, inRecQuery, &status); |
if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus)) |
|
|
if (status != DRMS_SUCCESS || inRecSet == NULL) |
|
{ | { |
fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status); |
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 0; | return 0; |
} | } |
|
if (gaprecset->n == 0) |
if (inRecSet->n == 0) |
|
{ | { |
printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]); |
fprintf(stderr,"ERROR: gapfile recordset contains no records\n"); |
return 0; | return 0; |
} | } |
|
if (gaprecset->n > 1) |
if (inRecSet->n > 1) |
|
{ | { |
printf("nRecs > 1 not yet supported.\n"); |
fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n"); |
return 0; | return 0; |
} | } |
|
gapseg = drms_segment_lookupnum(gaprecset->records[0], 0); |
|
if (gapseg != NULL) |
inrec=inRecSet->records[0]; |
gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status); |
// insert checks on input keywords here |
if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL) |
|
|
int lmin=drms_getkey_int(inrec, "LMIN", NULL); |
|
int lmax=drms_getkey_int(inrec, "LMAX", NULL); |
|
if (lmin != lmax) |
|
{ | { |
printf("lmin != lmax not yet supported.\n"); |
fprintf(stderr, "problem reading gaps segment: status = %d\n", status); |
return 0; | return 0; |
} | } |
lmode=lmin; |
agaps=(int *)(gaparr->data); |
|
gapsize=gaparr->axis[0]; |
fits_open_file(&fitsptr, gapfile, READONLY, &fstatus); |
} |
fits_report_error(stderr, fstatus); |
else |
|
{ |
fits_get_img_size(fitsptr, 1, fdims, &fstatus); | fits_get_img_size(fitsptr, 1, fdims, &fstatus); |
gapsize=fdims[0]; | gapsize=fdims[0]; |
agaps=(int *)(malloc(gapsize*sizeof(int))); | agaps=(int *)(malloc(gapsize*sizeof(int))); |
|
|
fits_report_error(stderr, fstatus); | fits_report_error(stderr, fstatus); |
return 0; | return 0; |
} | } |
|
} |
| |
printf("gapfile read, gapsize = %ld\n",gapsize); |
printf("gapfile read, gapsize = %d\n",gapsize); |
| |
data_dim=gapsize; | data_dim=gapsize; |
if (n_sampled_out>gapsize) | if (n_sampled_out>gapsize) |
data_dim=n_sampled_out; | data_dim=n_sampled_out; |
| |
/* Read location of jumps. */ | /* Read location of jumps. */ |
if (!strncmp(sectionfile,"none",4)) |
if (!strcmp(sectionfile,kNOTSPECIFIED)) |
{ | { |
/* No sections file specified. Assume that the whole | /* No sections file specified. Assume that the whole |
time series in one section. */ | time series in one section. */ |
|
|
} | } |
else | else |
{ | { |
if (read_section_file(sectionfile, &num_sections, §ions)) |
int secstat; |
|
if (secstat=read_section_file(sectionfile, &num_sections, §ions)) |
|
{ |
|
DRMS_RecordSet_t *secrecset = drms_open_records(drms_env, sectionfile, &status); |
|
if (status != DRMS_SUCCESS || secrecset == NULL) |
|
{ |
|
fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status); |
|
return 0; |
|
} |
|
if (secrecset->n == 0) |
|
{ |
|
fprintf(stderr,"ERROR: sectionfile recordset contains no records\n"); |
|
return 0; |
|
} |
|
if (secrecset->n > 1) |
|
{ |
|
fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n"); |
|
return 0; |
|
} |
|
num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL); |
|
if (num_sections < 1) |
|
{ |
|
fprintf(stderr,"ERROR: Invalid NSECS keywords\n"); |
|
return 0; |
|
} |
|
sections = (int *)malloc(2*(num_sections)*sizeof(int)); |
|
char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL); |
|
sscanf(strtok(sectkey," \n"),"%d",&(sections[0])); |
|
sscanf(strtok(NULL," \n"),"%d",&(sections[1])); |
|
int i; |
|
for (i=2 ;i<2*(num_sections); i+=2) |
|
{ |
|
sscanf(strtok(NULL," \n"), "%d", &(sections)[i]); |
|
sscanf(strtok(NULL," \n"), "%d", &(sections)[i+1]); |
|
|
|
if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1])) |
{ | { |
fprintf(stderr, "ERROR: problem with section file\n"); |
fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n"); |
return 0; | return 0; |
} | } |
} | } |
|
free(sectkey); |
|
} |
|
} |
| |
printf("num_sections= %d.\n",num_sections); |
printf("num_sections = %d\n",num_sections); |
| |
/* allocate temporary storage */ | /* allocate temporary storage */ |
gaps=(int *)(malloc(gapsize*sizeof(int))); | gaps=(int *)(malloc(gapsize*sizeof(int))); |
|
|
data=(float *)(fftwf_malloc(2*data_dim*sizeof(float))); | data=(float *)(fftwf_malloc(2*data_dim*sizeof(float))); |
ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float)); | 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) |
|
{ |
|
fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status); |
|
return 0; |
|
} |
|
|
|
if (inRecSet->n == 0) |
|
{ |
|
printf("WARNING: input recordset contains no records\nmodule %s successful completion\n", cmdparams.argv[0]); |
|
return 0; |
|
} |
|
|
|
if (inRecSet->n > 1) |
|
{ |
|
printf("ERROR: nRecs > 1 not yet supported.\n"); |
|
return 0; |
|
} |
|
|
|
|
|
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) |
|
{ |
|
fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]); |
|
drms_close_records(inRecSet, DRMS_FREE_RECORD); |
|
return 0; |
|
} |
|
} |
|
|
|
|
|
// insert loop over input here |
|
int iRec=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); | tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL); |
tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL); | tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL); |
tstep=drms_getkey_float(inrec, "T_STEP", NULL); | tstep=drms_getkey_float(inrec, "T_STEP", NULL); |
|
|
return 0; | return 0; |
} | } |
| |
if (extract_gaps(n_in, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections)) |
tsample=tstep; |
|
|
|
|
|
|
|
// 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"); | fprintf(stderr, "ERROR: problem in extract_gaps\n"); |
return 0; | return 0; |
|
|
n_sampled_out = n_sampled_in; | n_sampled_out = n_sampled_in; |
df1 = tsamplex*n_sampled_out; | df1 = tsamplex*n_sampled_out; |
| |
double tstartout=tstart+ifirst*tsample; |
|
double tstopout=tstartout+df1; |
|
|
|
if (fftout || fft1out || powout || arpowout || mavgout) | if (fftout || fft1out || powout || arpowout || mavgout) |
plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, | plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, |
(fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE); | (fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE); |
| |
|
/* Set sum to zero before starting */ |
|
if (mavgout) |
|
msum = (float *)calloc(n_sampled_out/2,sizeof(float)); |
|
|
|
|
|
tstartout=tstart+ifirst*tsample; |
|
tstopout=tstartout+df1; |
|
|
|
|
/* Create output data set. */ | /* Create output data set. */ |
if (tsout || fftout) | if (tsout || fftout) |
{ | { |
Line 378 double tstopout=tstartout+df1; |
|
Line 492 double tstopout=tstartout+df1; |
|
} | } |
| |
| |
/* Set sum to zero before starting */ |
|
if (mavgout) |
|
msum = (float *)calloc(n_sampled_out/2,sizeof(float)); |
|
|
|
/* Loops through all the data and update the gaps to reflect missing | /* Loops through all the data and update the gaps to reflect missing |
data. */ | data. */ |
start_index[0] = 0; |
/* the above comment is from the original code. it refers to an ifdef |
start_index[1] = 0; |
that would read in all the data and run extract_data on it in case there |
counts[0] = 2*n_in; |
were nan's not removed in the input window function. the present version |
counts[1] = 1; |
assumes the window function is correct. */ |
intervals[0] = 0; |
|
intervals[1] = 0; |
|
| |
| |
|
|
//insert loop over input here |
|
int iRec=0; |
|
inrec = inRecSet->records[iRec]; |
|
if (seginflag) | if (seginflag) |
segin = drms_segment_lookup(inrec, segnamein); | segin = drms_segment_lookup(inrec, segnamein); |
else | else |
|
|
inarr = drms_segment_read(segin, usetype, &status); | inarr = drms_segment_read(segin, usetype, &status); |
if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL) | if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL) |
{ | { |
fprintf(stderr, "problem reading input segment: iRec = %d, status = %d\n", iRec, status); |
fprintf(stderr, "ERROR: problem reading input segment: iRec = %d, status = %d\n", iRec, status); |
// continue; |
|
return 0; | return 0; |
} | } |
datarr=(float *)(inarr->data); | datarr=(float *)(inarr->data); |
|
|
if (status != DRMS_SUCCESS || outrec==NULL) | if (status != DRMS_SUCCESS || outrec==NULL) |
{ | { |
fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status); | fprintf(stderr,"ERROR: unable to open record in output dataseries, status = %d\n", status); |
// error=2; |
|
// break; |
|
return 0; | return 0; |
} | } |
| |
|
|
cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2, | cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2, |
fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, | fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, |
&ar_order, (_Complex float *)ar_coeff); | &ar_order, (_Complex float *)ar_coeff); |
|
/* |
{ | { |
FILE *fh; | FILE *fh; |
unsigned char c; | unsigned char c; |
|
|
} | } |
fclose(fh); | fclose(fh); |
} | } |
|
*/ |
// (*history)("AR model order used = %d.\n",ar_order); | // (*history)("AR model order used = %d.\n",ar_order); |
} | } |
| |
|
|
// (*history)("Preparing output.\n"); | // (*history)("Preparing output.\n"); |
| |
if (tsout || fftout) | if (tsout || fftout) |
memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], |
// memcpy(&out_data[2*m*n_sampled_out], &data[2*ifirst], 2*n_sampled_out*sizeof(float)); |
2*n_sampled_out*sizeof(float)); |
// the above line is above, already moved shifted ifirst to start of data, call should be this: |
|
memcpy(&out_data[2*m*n_sampled_out], data, 2*n_sampled_out*sizeof(float)); |
|
|
else if (fft1out) | else if (fft1out) |
{ | { |
/* Do positive m */ | /* Do positive m */ |
Line 625 drms_setkey_time(outrec, "T_START", tsta |
|
Line 730 drms_setkey_time(outrec, "T_START", tsta |
|
drms_setkey_time(outrec, "T_STOP", tstopout); | drms_setkey_time(outrec, "T_STOP", tstopout); |
drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2); | drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2); |
drms_setkey_float(outrec, "T_STEP", tsamplex); | drms_setkey_float(outrec, "T_STEP", tsamplex); |
|
drms_setkey_int(outrec, "NDT", n_sampled_out); |
|
|
| |
| |
tnow = (double)time(NULL); | tnow = (double)time(NULL); |
Line 952 static int read_section_file(char *filen |
|
Line 1059 static int read_section_file(char *filen |
|
(i>0 && (*sections)[i]<=(*sections)[i-1])) | (i>0 && (*sections)[i]<=(*sections)[i-1])) |
{ | { |
fprintf(stderr,"Invalid sections file, sections overlap.\n"); | fprintf(stderr,"Invalid sections file, sections overlap.\n"); |
return 1; |
return 2; |
} | } |
} | } |
fclose(fh); | fclose(fh); |