| 1 |
/* this JSOC module is a port of an SOI module written by Rasmus Larsen. |
| 2 |
* ported by Tim Larson. |
| 3 |
* tim doubts this works correctly with nskip or navg |
| 4 |
*/ |
| 5 |
|
| 6 |
/* |
| 7 |
* Detrending/gapfilling/differencing module |
| 8 |
* ts_fiddle_rml upgrades ts_fiddle |
| 9 |
*/ |
| 10 |
|
| 11 |
#include <stdio.h> |
| 12 |
#include <stdlib.h> |
| 13 |
#include <strings.h> |
| 14 |
#include <time.h> |
| 15 |
#include <sys/time.h> |
| 16 |
#include <sys/resource.h> |
| 17 |
#include <math.h> |
| 18 |
#include <assert.h> |
| 19 |
#include <fftw3.h> |
| 20 |
|
| 21 |
#include "jsoc_main.h" |
| 22 |
#include "fitsio.h" |
| 23 |
|
| 24 |
#define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0])) |
| 25 |
#define kNOTSPECIFIED "not specified" |
| 26 |
|
| 27 |
char *module_name = "jtsfiddle"; |
| 28 |
char *cvsinfo_jtsfiddle = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.14 2013/04/28 08:05:21 tplarson Exp $"; |
| 29 |
|
| 30 |
/* Command line arguments: */ |
| 31 |
ModuleArgs_t module_args[] = |
| 32 |
{ |
| 33 |
{ARG_STRING, "in", "", "input data records"}, |
| 34 |
{ARG_STRING, "tsout", kNOTSPECIFIED, "output dataseries for timeseries"}, |
| 35 |
{ARG_STRING, "powout", kNOTSPECIFIED, "output dataseries for power spectra"}, |
| 36 |
{ARG_STRING, "fftout", kNOTSPECIFIED, "output dataseries for fft's"}, |
| 37 |
{ARG_STRING, "fft1out", kNOTSPECIFIED, "output dataseries for fft's reordered"}, |
| 38 |
{ARG_STRING, "arpowout", kNOTSPECIFIED, "output dataseries for AR power"}, |
| 39 |
{ARG_STRING, "mavgout", kNOTSPECIFIED, "output dataseries for m-averaged power spectra"}, |
| 40 |
// {ARG_STRING, "out", "", "output data series"}, |
| 41 |
{ARG_STRING, "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"}, |
| 42 |
{ARG_STRING, "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"}, |
| 43 |
{ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"}, |
| 44 |
{ARG_STRING, "srclink", "SRCDATA", "name of link to source data"}, |
| 45 |
{ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"}, |
| 46 |
{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}, |
| 47 |
{ARG_STRING, "TAG", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, usually used as an additional primekey"}, |
| 48 |
{ARG_STRING, "VERSION", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, useful for selecting obsolete records"}, |
| 49 |
|
| 50 |
{ARG_STRING, "LOGFILE", "stdout", ""}, |
| 51 |
{ARG_STRING, "GAPFILE", "", "record or file containing window function"}, |
| 52 |
{ARG_STRING, "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"}, |
| 53 |
|
| 54 |
{ARG_INT, "NDT", "-1", "", ""}, |
| 55 |
// {ARG_FLOAT, "TSAMPLE", "60", "", ""}, |
| 56 |
{ARG_INT, "IFILL", "1", "", ""}, |
| 57 |
{ARG_INT, "IFIRST", "0", "", ""}, |
| 58 |
// {ARG_INT, "FLIPM", "0", "", ""}, |
| 59 |
{ARG_INT, "MAX_AR_ORDER", "360", "", ""}, |
| 60 |
{ARG_INT, "FU_FIT_EDGES", "1", "", ""}, |
| 61 |
{ARG_INT, "FU_ITER", "3", "", ""}, |
| 62 |
{ARG_INT, "FU_MIN_PERCENT", "90", "", ""}, |
| 63 |
{ARG_FLOAT, "CDETREND", "50.0", "", ""}, |
| 64 |
{ARG_INT, "DETREND_LENGTH", "1600", "", ""}, |
| 65 |
{ARG_INT, "DETREND_SKIP", "1440", "", ""}, |
| 66 |
{ARG_INT, "DETREND_FIRST", "0", "", ""}, |
| 67 |
/* |
| 68 |
{ARG_INT, "TSOUT", "0", "", ""}, |
| 69 |
{ARG_INT, "FFTOUT", "0", "", ""}, |
| 70 |
{ARG_INT, "FFT1OUT", "0", "", ""}, |
| 71 |
{ARG_INT, "POWOUT", "0", "", ""}, |
| 72 |
{ARG_INT, "ARPOWOUT", "0", "", ""}, |
| 73 |
{ARG_INT, "MAVGOUT", "0", "", ""}, |
| 74 |
*/ |
| 75 |
{ARG_INT, "NAVG", "0", "", ""}, |
| 76 |
{ARG_INT, "NSKIP", "0", "", ""}, |
| 77 |
{ARG_INT, "NOFF", "0", "", ""}, |
| 78 |
{ARG_INT, "IMIN", "0", "", ""}, |
| 79 |
{ARG_INT, "IMAX", "-1", "", ""}, |
| 80 |
{ARG_END}, |
| 81 |
}; |
| 82 |
|
| 83 |
#include "saveparm.c" |
| 84 |
#include "timing.c" |
| 85 |
#include "set_history.c" |
| 86 |
|
| 87 |
extern void cdetrend_discontig( int n, _Complex float *data, int *isgood, |
| 88 |
int degree, int length, int skip, |
| 89 |
int m, int *sect_last, int detrend_first); |
| 90 |
|
| 91 |
int cfahlman_ulrych(int n, _Complex float *data, int *isgood, |
| 92 |
int minpercentage, int maxorder, int iterations, |
| 93 |
int padends, int *order, _Complex float *ar_coeff); |
| 94 |
|
| 95 |
//char *getdetrendversion(void); |
| 96 |
//char *getgapfillversion(void); |
| 97 |
|
| 98 |
/* global variables holding the values of the command line variables. */ |
| 99 |
static char *logfile, *gapfile, *sectionfile; |
| 100 |
static int lmode, n_sampled_out, ifill, flipm; |
| 101 |
static int ar_order, max_ar_order, fu_iter, fu_fit_edges; |
| 102 |
static int fu_min_percent; |
| 103 |
static int detrend_length, detrend_skip, detrend_first; |
| 104 |
static float tsample, detrend_const; |
| 105 |
static int ifirst, tsout,fftout, fft1out, powout, arpowout, mavgout; |
| 106 |
static int navg, nskip, noff, imin, imax; |
| 107 |
|
| 108 |
|
| 109 |
/* Prototypes for local functions. */ |
| 110 |
static double splitting(int l, int m); |
| 111 |
|
| 112 |
/* Prototypes for external functions */ |
| 113 |
extern void cffti_(int *, float *); |
| 114 |
extern void cfftb_(int *, float *, float *); |
| 115 |
static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, |
| 116 |
float tsample_in, float *tsample_out, |
| 117 |
int *num_sections, int *sections); |
| 118 |
static void extract_data(int n_in, int *gaps, float *data_in, float *data_out); |
| 119 |
static int read_section_file(char *filename, int *num_section, int **section); |
| 120 |
|
| 121 |
/************ Here begins the main module **************/ |
| 122 |
int DoIt(void) |
| 123 |
{ |
| 124 |
int i; |
| 125 |
int m, n_in, n_sampled_in; |
| 126 |
TIME epoch, step, start, stop; |
| 127 |
int gapsize, istart, istop, data_dim, detrend_order; |
| 128 |
int npow, mshift; |
| 129 |
float tsamplex, df1, c; |
| 130 |
int *agaps; |
| 131 |
int *gaps, *gaps2; |
| 132 |
float *msum, *pow; |
| 133 |
float *data, *out_data, *in_data; |
| 134 |
fftwf_plan plan; |
| 135 |
int num_sections, *sections; |
| 136 |
float *ar_coeff=NULL; |
| 137 |
float *datarr, *datptr, *outptr; |
| 138 |
float *datahold; |
| 139 |
char tstartstr[100]; |
| 140 |
int lmin, lmax; |
| 141 |
|
| 142 |
int fstatus = 0; |
| 143 |
fitsfile *fitsptr; |
| 144 |
long fdims[1]; |
| 145 |
long fpixel[]={1}; |
| 146 |
int *anynul=0; |
| 147 |
|
| 148 |
int instart[2], inend[2]; |
| 149 |
int tslength[2], tsstart[2], tsend[2], tstotal[2]; |
| 150 |
int fft1length[2], fft1start[2], fft1end[2], fft1total[2]; |
| 151 |
int powlength[2], powstart[2], powend[2], powtotal[2]; |
| 152 |
int status=DRMS_SUCCESS; |
| 153 |
int newstat=0; |
| 154 |
|
| 155 |
DRMS_Segment_t *segin = NULL; |
| 156 |
DRMS_Segment_t *segout = NULL; |
| 157 |
DRMS_Array_t *inarr = NULL; |
| 158 |
DRMS_Array_t *outarr = NULL; |
| 159 |
DRMS_Record_t *inrec = NULL; |
| 160 |
DRMS_Record_t *outrec = NULL; |
| 161 |
DRMS_Type_t usetype = DRMS_TYPE_FLOAT; |
| 162 |
DRMS_RecLifetime_t lifetime; |
| 163 |
long long histrecnum = -1; |
| 164 |
DRMS_Segment_t *gapseg = NULL; |
| 165 |
DRMS_Array_t *gaparr = NULL; |
| 166 |
|
| 167 |
HIterator_t outKey_list; |
| 168 |
DRMS_Keyword_t *outKey; |
| 169 |
TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */ |
| 170 |
|
| 171 |
int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); |
| 172 |
int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); |
| 173 |
|
| 174 |
double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout; |
| 175 |
double wt0, wt1, wt2, wt3, wt; |
| 176 |
double ut0, ut1, ut2, ut3, ut; |
| 177 |
double st0, st1, st2, st3, st; |
| 178 |
double ct0, ct1, ct2, ct3, ct; |
| 179 |
|
| 180 |
wt0=getwalltime(); |
| 181 |
ct0=getcputime(&ut0, &st0); |
| 182 |
|
| 183 |
/* Read input parameters. */ |
| 184 |
logfile = (char *)cmdparams_save_str (&cmdparams, "LOGFILE", &newstat); |
| 185 |
gapfile = (char *)cmdparams_save_str (&cmdparams, "GAPFILE", &newstat); |
| 186 |
sectionfile = (char *)cmdparams_save_str (&cmdparams, "SECTIONFILE", &newstat); |
| 187 |
n_sampled_out = cmdparams_save_int (&cmdparams, "NDT", &newstat); |
| 188 |
// tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat); |
| 189 |
ifill = cmdparams_save_int (&cmdparams, "IFILL", &newstat); |
| 190 |
max_ar_order= cmdparams_save_int (&cmdparams, "MAX_AR_ORDER", &newstat); |
| 191 |
fu_fit_edges= cmdparams_save_int (&cmdparams, "FU_FIT_EDGES", &newstat); |
| 192 |
fu_iter = cmdparams_save_int (&cmdparams, "FU_ITER", &newstat); |
| 193 |
fu_min_percent = cmdparams_save_int (&cmdparams, "FU_MIN_PERCENT", &newstat); |
| 194 |
ifirst = cmdparams_save_int (&cmdparams, "IFIRST", &newstat); |
| 195 |
// flipm = cmdparams_save_int (&cmdparams, "FLIPM", &newstat); |
| 196 |
detrend_const = cmdparams_save_float(&cmdparams, "CDETREND", &newstat); |
| 197 |
detrend_length = cmdparams_save_int (&cmdparams, "DETREND_LENGTH", &newstat); |
| 198 |
detrend_skip = cmdparams_save_int (&cmdparams, "DETREND_SKIP", &newstat); |
| 199 |
detrend_first = cmdparams_save_int (&cmdparams, "DETREND_FIRST", &newstat); |
| 200 |
/* |
| 201 |
tsout = cmdparams_save_int (&cmdparams, "TSOUT", &newstat); |
| 202 |
fftout = cmdparams_save_int (&cmdparams, "FFTOUT", &newstat); |
| 203 |
fft1out = cmdparams_save_int (&cmdparams, "FFT1OUT", &newstat); |
| 204 |
powout = cmdparams_save_int (&cmdparams, "POWOUT", &newstat); |
| 205 |
arpowout = cmdparams_save_int (&cmdparams, "ARPOWOUT", &newstat); |
| 206 |
mavgout = cmdparams_save_int (&cmdparams, "MAVGOUT", &newstat); |
| 207 |
*/ |
| 208 |
navg = cmdparams_save_int (&cmdparams, "NAVG", &newstat); |
| 209 |
nskip = cmdparams_save_int (&cmdparams, "NSKIP", &newstat); |
| 210 |
noff = cmdparams_save_int (&cmdparams, "NOFF", &newstat); |
| 211 |
imin = cmdparams_save_int (&cmdparams, "IMIN", &newstat); |
| 212 |
imax = cmdparams_save_int (&cmdparams, "IMAX", &newstat); |
| 213 |
|
| 214 |
char *inrecquery = (char *)cmdparams_save_str(&cmdparams, "in", &newstat); |
| 215 |
// char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat); |
| 216 |
char *segnamein = (char *)cmdparams_save_str(&cmdparams, "segin", &newstat); |
| 217 |
char *segnameout = (char *)cmdparams_save_str(&cmdparams, "segout", &newstat); |
| 218 |
int seginflag = strcmp(kNOTSPECIFIED, segnamein); |
| 219 |
int segoutflag = strcmp(kNOTSPECIFIED, segnameout); |
| 220 |
char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat); |
| 221 |
char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat); |
| 222 |
int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat); |
| 223 |
int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat); |
| 224 |
if (permflag) |
| 225 |
lifetime = DRMS_PERMANENT; |
| 226 |
else |
| 227 |
lifetime = DRMS_TRANSIENT; |
| 228 |
|
| 229 |
char *tag = (char *)cmdparams_save_str(&cmdparams, "TAG", &newstat); |
| 230 |
int tagflag = strcmp(kNOTSPECIFIED, tag); |
| 231 |
char *version = (char *)cmdparams_save_str(&cmdparams, "VERSION", &newstat); |
| 232 |
int verflag = strcmp(kNOTSPECIFIED, version); |
| 233 |
|
| 234 |
char *tsseries = (char *)cmdparams_save_str(&cmdparams, "tsout", &newstat); |
| 235 |
tsout = strcmp(kNOTSPECIFIED, tsseries); |
| 236 |
char *powseries = (char *)cmdparams_save_str(&cmdparams, "powout", &newstat); |
| 237 |
powout = strcmp(kNOTSPECIFIED, powseries); |
| 238 |
char *fftseries = (char *)cmdparams_save_str(&cmdparams, "fftout", &newstat); |
| 239 |
fftout = strcmp(kNOTSPECIFIED, fftseries); |
| 240 |
char *fft1series = (char *)cmdparams_save_str(&cmdparams, "fft1out", &newstat); |
| 241 |
fft1out = strcmp(kNOTSPECIFIED, fft1series); |
| 242 |
char *arpowseries = (char *)cmdparams_save_str(&cmdparams, "arpowout", &newstat); |
| 243 |
arpowout = strcmp(kNOTSPECIFIED, arpowseries); |
| 244 |
char *mavgseries = (char *)cmdparams_save_str(&cmdparams, "mavgout", &newstat); |
| 245 |
mavgout = strcmp(kNOTSPECIFIED, mavgseries); |
| 246 |
|
| 247 |
char *serieslist[6]; |
| 248 |
enum seriesindex {TSOUT, FFTOUT, FFT1OUT, POWOUT, ARPOWOUT, MAVGOUT}; |
| 249 |
int flagarr[6] = {tsout, fftout, fft1out, powout, arpowout, mavgout}; |
| 250 |
int mfliparr[6] = {0, 0, 0, 0, 0, 0}; |
| 251 |
int msignarr[6] = {0, 0, 0, 0, 0, 0}; |
| 252 |
serieslist[TSOUT]=tsseries; |
| 253 |
serieslist[FFTOUT]=fftseries; |
| 254 |
serieslist[FFT1OUT]=fft1series; |
| 255 |
serieslist[POWOUT]=powseries; |
| 256 |
serieslist[ARPOWOUT]=arpowseries; |
| 257 |
serieslist[MAVGOUT]=mavgseries; |
| 258 |
long long histrecnumarr[6]={-1, -1, -1, -1, -1, -1}; |
| 259 |
DRMS_RecordSet_t *outrecsetarr[6]={NULL, NULL, NULL, NULL, NULL, NULL}; |
| 260 |
DRMS_Record_t *outrecarr[6]={NULL, NULL, NULL, NULL, NULL, NULL}; |
| 261 |
DRMS_Segment_t *segoutarr[6]={NULL, NULL, NULL, NULL, NULL, NULL}; |
| 262 |
|
| 263 |
/**** Sanity check of input parameters. ****/ |
| 264 |
|
| 265 |
if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0)) |
| 266 |
{ |
| 267 |
fprintf(stderr, "ERROR: must specify an output dataseries\n"); |
| 268 |
return 1; |
| 269 |
} |
| 270 |
|
| 271 |
if (navg <= 0) |
| 272 |
navg=1; |
| 273 |
if (nskip <= 0) |
| 274 |
nskip=1; |
| 275 |
if ((navg != 1) && (nskip != 1)) |
| 276 |
{ |
| 277 |
fprintf(stderr, "ERROR: one of navg and nskip must equal 1\n"); |
| 278 |
return 1; |
| 279 |
} |
| 280 |
if (noff < 0 || noff >= nskip) |
| 281 |
{ |
| 282 |
fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n"); |
| 283 |
return 1; |
| 284 |
} |
| 285 |
if (arpowout && !ifill) |
| 286 |
{ |
| 287 |
fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n"); |
| 288 |
return 1; |
| 289 |
} |
| 290 |
if (fu_min_percent <= 0 || fu_min_percent > 100) |
| 291 |
{ |
| 292 |
fprintf(stderr, "ERROR: FU_MIN_PERCENT must be > 0 and <= 100\n"); |
| 293 |
return 1; |
| 294 |
} |
| 295 |
|
| 296 |
if (newstat) |
| 297 |
{ |
| 298 |
fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat); |
| 299 |
cpsave_decode_error(newstat); |
| 300 |
return 1; |
| 301 |
} |
| 302 |
else if (savestrlen != strlen(savestr)) |
| 303 |
{ |
| 304 |
fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr)); |
| 305 |
return 1; |
| 306 |
} |
| 307 |
|
| 308 |
|
| 309 |
// these must be present in the output dataseries and variable, not links or constants |
| 310 |
// 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... |
| 311 |
char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"}; |
| 312 |
int is, ishold, itest, mflip; |
| 313 |
char *holdseries=""; |
| 314 |
/* |
| 315 |
char *cvsinfo; |
| 316 |
cvsinfo = (char *)malloc(1024); |
| 317 |
strcpy(cvsinfo,"$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.14 2013/04/28 08:05:21 tplarson Exp $"); |
| 318 |
strcat(cvsinfo,"\n"); |
| 319 |
strcat(cvsinfo,getdetrendversion()); |
| 320 |
strcat(cvsinfo,"\n"); |
| 321 |
strcat(cvsinfo,getgapfillversion()); |
| 322 |
*/ |
| 323 |
for (is=0;is<6;is++) |
| 324 |
{ |
| 325 |
|
| 326 |
if (!flagarr[is]) |
| 327 |
continue; |
| 328 |
|
| 329 |
DRMS_Record_t *tempoutrec = drms_create_record(drms_env, |
| 330 |
serieslist[is], |
| 331 |
DRMS_TRANSIENT, |
| 332 |
&status); |
| 333 |
|
| 334 |
if (status != DRMS_SUCCESS) |
| 335 |
{ |
| 336 |
fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status); |
| 337 |
return 1; |
| 338 |
} |
| 339 |
|
| 340 |
DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname); |
| 341 |
|
| 342 |
|
| 343 |
// set up ancillary dataseries for processing metadata |
| 344 |
if (histlink != NULL) |
| 345 |
{ |
| 346 |
if (strcmp(holdseries, histlink->info->target_series)) |
| 347 |
{ |
| 348 |
ishold=is; |
| 349 |
holdseries=strdup(histlink->info->target_series); |
| 350 |
histrecnumarr[is]=set_history(histlink); |
| 351 |
if (histrecnumarr[is] < 0) |
| 352 |
{ |
| 353 |
fprintf(stderr, "ERROR: problem creating record in history dataseries for output dataseries %s\n", serieslist[is]); |
| 354 |
drms_close_record(tempoutrec, DRMS_FREE_RECORD); |
| 355 |
return 1; |
| 356 |
} |
| 357 |
} |
| 358 |
else |
| 359 |
{ |
| 360 |
histrecnumarr[is]=histrecnumarr[ishold]; |
| 361 |
} |
| 362 |
} |
| 363 |
else |
| 364 |
{ |
| 365 |
fprintf(stderr,"WARNING: could not find history link in output dataseries %s\n", serieslist[is]); |
| 366 |
} |
| 367 |
|
| 368 |
for (itest=0; itest < ARRLENGTH(outchecklist); itest++) |
| 369 |
{ |
| 370 |
DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]); |
| 371 |
if (outkeytest == NULL || outkeytest->info->islink || outkeytest->info->recscope == 1) |
| 372 |
{ |
| 373 |
fprintf(stderr, "ERROR: output keyword %s in series %s is either missing, constant, or a link\n", outchecklist[itest], serieslist[is]); |
| 374 |
drms_close_record(tempoutrec, DRMS_FREE_RECORD); |
| 375 |
return 1; |
| 376 |
} |
| 377 |
} |
| 378 |
|
| 379 |
mflip = drms_getkey_int(tempoutrec, "MFLIPPED", &status); |
| 380 |
if (status == DRMS_SUCCESS) |
| 381 |
mfliparr[is]=mflip; |
| 382 |
|
| 383 |
drms_close_record(tempoutrec, DRMS_FREE_RECORD); |
| 384 |
} |
| 385 |
|
| 386 |
// for efficiency require that tsseries and fftseries have the same value of MFLIPPED. this restriction may be lifted as noted below. |
| 387 |
if ((tsout && fftout) && (mfliparr[TSOUT] != mfliparr[FFTOUT])) |
| 388 |
{ |
| 389 |
fprintf(stderr, "ERROR: series %s and %s must have the same value of MFLIPPED\n", tsseries, fftseries); |
| 390 |
return 1; |
| 391 |
} |
| 392 |
|
| 393 |
if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus)) |
| 394 |
{ |
| 395 |
DRMS_RecordSet_t *gaprecset = drms_open_records(drms_env, gapfile, &status); |
| 396 |
if (status != DRMS_SUCCESS || gaprecset == NULL) |
| 397 |
{ |
| 398 |
fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status); |
| 399 |
return 1; |
| 400 |
} |
| 401 |
if (gaprecset->n == 0) |
| 402 |
{ |
| 403 |
fprintf(stderr,"ERROR: gapfile recordset contains no records\n"); |
| 404 |
return 1; |
| 405 |
} |
| 406 |
if (gaprecset->n > 1) |
| 407 |
{ |
| 408 |
fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n"); |
| 409 |
return 1; |
| 410 |
} |
| 411 |
gapseg = drms_segment_lookupnum(gaprecset->records[0], 0); |
| 412 |
if (gapseg != NULL) |
| 413 |
gaparr = drms_segment_read(gapseg, DRMS_TYPE_INT, &status); |
| 414 |
if (status != DRMS_SUCCESS || gaparr == NULL || gapseg == NULL) |
| 415 |
{ |
| 416 |
fprintf(stderr, "problem reading gaps segment: status = %d\n", status); |
| 417 |
return 1; |
| 418 |
} |
| 419 |
agaps=(int *)(gaparr->data); |
| 420 |
gapsize=gaparr->axis[0]; |
| 421 |
} |
| 422 |
else |
| 423 |
{ |
| 424 |
fits_get_img_size(fitsptr, 1, fdims, &fstatus); |
| 425 |
gapsize=fdims[0]; |
| 426 |
agaps=(int *)(malloc(gapsize*sizeof(int))); |
| 427 |
fits_read_pix(fitsptr, TINT, fpixel, gapsize, NULL, agaps, anynul, &fstatus); |
| 428 |
fits_close_file(fitsptr, &fstatus); |
| 429 |
if (fstatus) |
| 430 |
{ |
| 431 |
fits_report_error(stderr, fstatus); |
| 432 |
return 1; |
| 433 |
} |
| 434 |
} |
| 435 |
|
| 436 |
if (verbflag) |
| 437 |
printf("gapfile read, gapsize = %d\n",gapsize); |
| 438 |
|
| 439 |
data_dim=gapsize; |
| 440 |
if (n_sampled_out>gapsize) |
| 441 |
data_dim=n_sampled_out; |
| 442 |
|
| 443 |
/* allocate temporary storage */ |
| 444 |
gaps=(int *)(malloc(gapsize*sizeof(int))); |
| 445 |
gaps2=(int *)(malloc(gapsize*sizeof(int))); |
| 446 |
data=(float *)(fftwf_malloc(2*data_dim*sizeof(float))); |
| 447 |
ar_coeff = (float *)malloc((max_ar_order+1)*2*sizeof(float)); |
| 448 |
// now done below |
| 449 |
// if (fft1out || powout || mavgout || arpowout) |
| 450 |
// pow=(float *)(malloc(n_sampled_out*sizeof(float))); |
| 451 |
|
| 452 |
/* Read location of jumps. */ |
| 453 |
if (!strcmp(sectionfile,kNOTSPECIFIED)) |
| 454 |
{ |
| 455 |
/* No sections file specified. Assume that the whole |
| 456 |
time series in one section. */ |
| 457 |
num_sections=1; |
| 458 |
sections = malloc(2*sizeof(int)); |
| 459 |
sections[0] = 0; |
| 460 |
sections[1] = data_dim-1; |
| 461 |
} |
| 462 |
else |
| 463 |
{ |
| 464 |
int secstat; |
| 465 |
if (secstat=read_section_file(sectionfile, &num_sections, §ions)) |
| 466 |
{ |
| 467 |
DRMS_RecordSet_t *secrecset = drms_open_records(drms_env, sectionfile, &status); |
| 468 |
if (status != DRMS_SUCCESS || secrecset == NULL) |
| 469 |
{ |
| 470 |
fprintf(stderr,"ERROR: problem reading sections: file status = %d, DRMS status = %d\n",secstat,status); |
| 471 |
return 1; |
| 472 |
} |
| 473 |
if (secrecset->n == 0) |
| 474 |
{ |
| 475 |
fprintf(stderr,"ERROR: sectionfile recordset contains no records\n"); |
| 476 |
return 1; |
| 477 |
} |
| 478 |
if (secrecset->n > 1) |
| 479 |
{ |
| 480 |
fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n"); |
| 481 |
return 1; |
| 482 |
} |
| 483 |
num_sections=drms_getkey_int(secrecset->records[0], "NSECS", NULL); |
| 484 |
if (num_sections < 1) |
| 485 |
{ |
| 486 |
fprintf(stderr,"ERROR: Invalid NSECS keywords\n"); |
| 487 |
return 1; |
| 488 |
} |
| 489 |
sections = (int *)malloc(2*(num_sections)*sizeof(int)); |
| 490 |
char *sectkey=drms_getkey_string(secrecset->records[0], "SECS", NULL); |
| 491 |
sscanf(strtok(sectkey," \n"),"%d",&(sections[0])); |
| 492 |
sscanf(strtok(NULL," \n"),"%d",&(sections[1])); |
| 493 |
int i; |
| 494 |
for (i=2 ;i<2*(num_sections); i+=2) |
| 495 |
{ |
| 496 |
sscanf(strtok(NULL," \n"), "%d", &(sections)[i]); |
| 497 |
sscanf(strtok(NULL," \n"), "%d", &(sections)[i+1]); |
| 498 |
|
| 499 |
if (sections[i]>sections[i+1] || (i>0 && sections[i]<=sections[i-1])) |
| 500 |
{ |
| 501 |
fprintf(stderr,"ERROR: Invalid SECS keyword, sections overlap.\n"); |
| 502 |
return 1; |
| 503 |
} |
| 504 |
} |
| 505 |
free(sectkey); |
| 506 |
} |
| 507 |
} |
| 508 |
|
| 509 |
if (verbflag) |
| 510 |
printf("num_sections = %d\n",num_sections); |
| 511 |
|
| 512 |
DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status); |
| 513 |
|
| 514 |
if (status != DRMS_SUCCESS || inrecset == NULL) |
| 515 |
{ |
| 516 |
fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status); |
| 517 |
return 1; |
| 518 |
} |
| 519 |
if (inrecset->n == 0) |
| 520 |
{ |
| 521 |
fprintf(stderr, "ERROR: input recordset contains no records\n"); |
| 522 |
return 1; |
| 523 |
} |
| 524 |
|
| 525 |
inrec=inrecset->records[0]; |
| 526 |
char *inchecklist[] = {"T_START", "T_STOP", "QUALITY", "LMIN", "LMAX", "T_STEP"}; |
| 527 |
for (itest=0; itest < ARRLENGTH(inchecklist); itest++) |
| 528 |
{ |
| 529 |
DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]); |
| 530 |
if (inkeytest == NULL) |
| 531 |
{ |
| 532 |
fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]); |
| 533 |
drms_close_records(inrecset, DRMS_FREE_RECORD); |
| 534 |
return 1; |
| 535 |
} |
| 536 |
} |
| 537 |
|
| 538 |
cadence=drms_getkey_float(inrec, "T_STEP", NULL); //assume constant across all input records |
| 539 |
tsample=cadence; |
| 540 |
tstartsave=drms_getkey_time(inrec, "T_START", NULL); //must be constant across all input records unless GAPFILE is implemented differently |
| 541 |
|
| 542 |
int mflipin=drms_getkey_int(inrec, "MFLIPPED", &status); |
| 543 |
if (status != DRMS_SUCCESS) mflipin=0; |
| 544 |
|
| 545 |
// already required mfliparr[TSOUT] == mfliparr[FFTOUT] above |
| 546 |
if (tsout && mfliparr[TSOUT] != mflipin) |
| 547 |
flipm=1; |
| 548 |
else if (fftout && mfliparr[FFTOUT] != mflipin) |
| 549 |
flipm=1; |
| 550 |
else |
| 551 |
flipm=0; |
| 552 |
|
| 553 |
for (is=2;is<5;is++) |
| 554 |
{ |
| 555 |
if (!flagarr[is]) |
| 556 |
continue; |
| 557 |
if (mfliparr[is] != (mflipin || flipm)) |
| 558 |
msignarr[is]=-1; |
| 559 |
else |
| 560 |
msignarr[is]=1; |
| 561 |
} |
| 562 |
|
| 563 |
if (mflipin) |
| 564 |
msignarr[MAVGOUT]=1; |
| 565 |
else |
| 566 |
msignarr[MAVGOUT]=-1; |
| 567 |
|
| 568 |
// changed n_in to gapsize in following call |
| 569 |
if (extract_gaps(gapsize, agaps, &n_sampled_in, gaps, tsample, &tsamplex, &num_sections, sections)) |
| 570 |
{ |
| 571 |
fprintf(stderr, "ERROR: problem in extract_gaps\n"); |
| 572 |
return 1; |
| 573 |
} |
| 574 |
|
| 575 |
/* Adjust detrend parameters according to the number of |
| 576 |
available points. */ |
| 577 |
if (n_sampled_in<detrend_length) |
| 578 |
{ |
| 579 |
detrend_length = n_sampled_in; |
| 580 |
detrend_skip = detrend_length; |
| 581 |
} |
| 582 |
|
| 583 |
int idtf; |
| 584 |
if (detrend_first > 0) |
| 585 |
{ |
| 586 |
for (idtf=0;idtf<detrend_first;idtf++) |
| 587 |
gaps[idtf]=0; |
| 588 |
} |
| 589 |
|
| 590 |
|
| 591 |
if (n_sampled_out == -1) |
| 592 |
n_sampled_out = n_sampled_in; |
| 593 |
df1 = tsamplex*n_sampled_out; |
| 594 |
|
| 595 |
if (fftout || fft1out || powout || arpowout || mavgout) |
| 596 |
plan = fftwf_plan_dft_1d(n_sampled_out, (fftwf_complex *)data, |
| 597 |
(fftwf_complex *)data, FFTW_BACKWARD, FFTW_MEASURE); |
| 598 |
|
| 599 |
/* Set sum to zero before starting */ |
| 600 |
if (mavgout) |
| 601 |
msum = (float *)malloc((n_sampled_out/2)*sizeof(float)); |
| 602 |
|
| 603 |
if (fft1out || powout || mavgout || arpowout) |
| 604 |
{ |
| 605 |
/* Set first and last output index if not set as a parameter. */ |
| 606 |
if (imax < imin) |
| 607 |
{ |
| 608 |
imin=0; |
| 609 |
imax=n_sampled_out/2-1; |
| 610 |
} |
| 611 |
npow=imax-imin+1; |
| 612 |
pow=(float *)(malloc(n_sampled_out*sizeof(float))); |
| 613 |
datahold=(float *)malloc(2*npow*sizeof(float)); |
| 614 |
} |
| 615 |
|
| 616 |
// needed to implement mfliparr[TSOUT] != mfliparr[FFTOUT] |
| 617 |
// float *datatemp=(float *)malloc(2*n_sampled_out*sizeof(float)); |
| 618 |
float *tmpptr=data; |
| 619 |
|
| 620 |
if (tsout || fftout) |
| 621 |
{ |
| 622 |
tslength[0]=tstotal[0]=2*n_sampled_out; |
| 623 |
tslength[1]=1; |
| 624 |
tsstart[0]=0; |
| 625 |
tsend[0]=2*n_sampled_out-1; |
| 626 |
} |
| 627 |
if (fft1out) |
| 628 |
{ |
| 629 |
fft1length[0]=fft1total[0]=2*npow; |
| 630 |
fft1length[1]=1; |
| 631 |
fft1start[0]=0; |
| 632 |
fft1end[0]=2*npow-1; |
| 633 |
} |
| 634 |
if (powout || arpowout || mavgout) |
| 635 |
{ |
| 636 |
powlength[0]=powtotal[0]=npow; |
| 637 |
powlength[1]=1; |
| 638 |
powstart[0]=0; |
| 639 |
powend[0]=npow-1; |
| 640 |
} |
| 641 |
instart[0]=0; |
| 642 |
inend[0]=2*gapsize - 1; |
| 643 |
|
| 644 |
status=drms_stage_records(inrecset, 1, 0); |
| 645 |
if (status != DRMS_SUCCESS) |
| 646 |
{ |
| 647 |
fprintf(stderr, "ERROR: drms_stage_records returned status = %d\n", status); |
| 648 |
return 1; |
| 649 |
} |
| 650 |
|
| 651 |
for (is=0; is<6; is++) |
| 652 |
{ |
| 653 |
if (!flagarr[is]) |
| 654 |
continue; |
| 655 |
|
| 656 |
outrecsetarr[is] = drms_create_records(drms_env, inrecset->n, serieslist[is], DRMS_PERMANENT, &status); |
| 657 |
|
| 658 |
if (status != DRMS_SUCCESS || outrecsetarr[is] == NULL) |
| 659 |
{ |
| 660 |
fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status); |
| 661 |
return 1; |
| 662 |
} |
| 663 |
} |
| 664 |
|
| 665 |
int irec; |
| 666 |
for (irec=0;irec<inrecset->n;irec++) |
| 667 |
{ |
| 668 |
|
| 669 |
inrec = inrecset->records[irec]; |
| 670 |
|
| 671 |
lmin=drms_getkey_int(inrec, "LMIN", NULL); |
| 672 |
lmax=drms_getkey_int(inrec, "LMAX", NULL); |
| 673 |
if (lmin != lmax) |
| 674 |
{ |
| 675 |
fprintf(stderr, "ERROR: lmin != lmax not yet supported.\n"); |
| 676 |
return 0; |
| 677 |
} |
| 678 |
lmode=lmin; |
| 679 |
|
| 680 |
if (verbflag > 1) |
| 681 |
{ |
| 682 |
printf("starting l=%d\n", lmode); |
| 683 |
} |
| 684 |
|
| 685 |
tstart=drms_getkey_time(inrec, "T_START", NULL); |
| 686 |
tstop =drms_getkey_time(inrec, "T_STOP", NULL); |
| 687 |
cadence=drms_getkey_float(inrec, "T_STEP", NULL); |
| 688 |
if (tstart != tstartsave) |
| 689 |
{ |
| 690 |
// to lift this restriction must be able to specify multiple gap files/records as input |
| 691 |
sprint_time(tstartstr, tstart, "TAI", 0); |
| 692 |
fprintf(stderr, "ERROR: all input records must have same T_START, record %d has %s\n", irec, tstartstr); |
| 693 |
return 0; |
| 694 |
} |
| 695 |
|
| 696 |
n_in=(tstop-tstart)/cadence; |
| 697 |
if (n_in != gapsize) |
| 698 |
{ |
| 699 |
fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, n_in=%d\n", gapsize, n_in); |
| 700 |
return 0; |
| 701 |
} |
| 702 |
|
| 703 |
tstartout=tstart+ifirst*tsample; |
| 704 |
tstopout=tstartout+df1; |
| 705 |
sprint_time(tstartstr, tstartout, "TAI", 0); |
| 706 |
|
| 707 |
if (seginflag) |
| 708 |
segin = drms_segment_lookup(inrec, segnamein); |
| 709 |
else |
| 710 |
segin = drms_segment_lookupnum(inrec, 0); |
| 711 |
if (segin == NULL) |
| 712 |
{ |
| 713 |
fprintf(stderr, "ERROR: problem looking up input segment: recnum = %lld\n", inrec->recnum); |
| 714 |
return 0; |
| 715 |
} |
| 716 |
|
| 717 |
for (is=0; is<6; is++) |
| 718 |
{ |
| 719 |
if (!flagarr[is]) |
| 720 |
continue; |
| 721 |
/* |
| 722 |
outrecarr[is] = drms_create_record(drms_env, serieslist[is], DRMS_PERMANENT, &status); |
| 723 |
if (status != DRMS_SUCCESS) |
| 724 |
{ |
| 725 |
fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", serieslist[is], status); |
| 726 |
return 0; |
| 727 |
} |
| 728 |
*/ |
| 729 |
outrec=outrecsetarr[is]->records[irec]; |
| 730 |
if (histrecnumarr[is] > 0) |
| 731 |
drms_setlink_static(outrec, histlinkname, histrecnumarr[is]); |
| 732 |
drms_setlink_static(outrec, srclinkname, inrec->recnum); |
| 733 |
|
| 734 |
if (segoutflag) |
| 735 |
segoutarr[is] = drms_segment_lookup(outrec, segnameout); |
| 736 |
else |
| 737 |
segoutarr[is] = drms_segment_lookupnum(outrec, 0); |
| 738 |
if (segoutarr[is] == NULL) |
| 739 |
{ |
| 740 |
fprintf(stderr, "ERROR: problem looking up outputput segment in series %s\n", serieslist[is]); |
| 741 |
return 0; |
| 742 |
} |
| 743 |
|
| 744 |
// the following is not needed if at least one of the first N-1 dimensions are 0 in the jsd, |
| 745 |
// or if the first N-1 dimensions in the jsd are always the dimensions desired. |
| 746 |
// it *is* needed any time one must override the jsd defaults |
| 747 |
/* |
| 748 |
switch(is) |
| 749 |
{ |
| 750 |
case TSOUT: |
| 751 |
case FFTOUT: |
| 752 |
segoutarr[is]->axis[0]=2*n_sampled_out; |
| 753 |
segoutarr[is]->axis[1]=lmode+1; |
| 754 |
break; |
| 755 |
case FFT1OUT: |
| 756 |
segoutarr[is]->axis[0]=2*npow; |
| 757 |
segoutarr[is]->axis[1]=2*lmode+1; |
| 758 |
break; |
| 759 |
case POWOUT: |
| 760 |
case ARPOWOUT: |
| 761 |
segoutarr[is]->axis[0]=npow; |
| 762 |
segoutarr[is]->axis[1]=2*lmode+1; |
| 763 |
break; |
| 764 |
case MAVGOUT: |
| 765 |
segoutarr[is]->axis[0]=npow; |
| 766 |
segoutarr[is]->axis[1]=1; |
| 767 |
break; |
| 768 |
default: |
| 769 |
; |
| 770 |
} |
| 771 |
*/ |
| 772 |
switch(is) |
| 773 |
{ |
| 774 |
case TSOUT: |
| 775 |
case FFTOUT: |
| 776 |
tstotal[1]=lmode+1; |
| 777 |
break; |
| 778 |
case FFT1OUT: |
| 779 |
fft1total[1]=2*lmode+1; |
| 780 |
break; |
| 781 |
case POWOUT: |
| 782 |
case ARPOWOUT: |
| 783 |
powtotal[1]=2*lmode+1; |
| 784 |
break; |
| 785 |
case MAVGOUT: |
| 786 |
default: |
| 787 |
; |
| 788 |
} |
| 789 |
|
| 790 |
} |
| 791 |
|
| 792 |
if (mavgout) |
| 793 |
for (i=0;i<n_sampled_out/2;i++) |
| 794 |
msum[i] = 0.0; |
| 795 |
|
| 796 |
/***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/ |
| 797 |
for (m=0; m<=lmode; m++) |
| 798 |
{ |
| 799 |
|
| 800 |
if (verbflag > 2) |
| 801 |
{ |
| 802 |
printf("starting m=%d\n", m); |
| 803 |
} |
| 804 |
|
| 805 |
instart[1]=m; |
| 806 |
inend[1]=m; |
| 807 |
inarr = drms_segment_readslice(segin, usetype, instart, inend, &status); |
| 808 |
if (status != DRMS_SUCCESS || inarr == NULL ) |
| 809 |
{ |
| 810 |
fprintf(stderr, "ERROR: problem reading input segment: status = %d, recnum = %lld\n", status, inrec->recnum); |
| 811 |
return 0; |
| 812 |
} |
| 813 |
in_data=(float *)(inarr->data); |
| 814 |
|
| 815 |
/* Extracts data copy */ |
| 816 |
extract_data(n_sampled_in, gaps, in_data, data); |
| 817 |
|
| 818 |
/* Detrend entire time series if required */ |
| 819 |
if (detrend_const != 0) |
| 820 |
{ |
| 821 |
detrend_order = floor(detrend_length/detrend_const)+2; |
| 822 |
cdetrend_discontig(n_sampled_in, (_Complex float *)data, gaps, |
| 823 |
detrend_order, detrend_length, detrend_skip, |
| 824 |
num_sections, sections, detrend_first); |
| 825 |
} |
| 826 |
|
| 827 |
/* Fill gaps in entire time series if required */ |
| 828 |
memcpy(gaps2, gaps, n_sampled_in*sizeof(int)); |
| 829 |
if (ifill != 0 && max_ar_order > 0) |
| 830 |
{ |
| 831 |
cfahlman_ulrych(n_sampled_in, (_Complex float *)data, gaps2, |
| 832 |
fu_min_percent, max_ar_order, fu_iter, fu_fit_edges, |
| 833 |
&ar_order, (_Complex float *)ar_coeff); |
| 834 |
} |
| 835 |
|
| 836 |
drms_free_array(inarr); |
| 837 |
|
| 838 |
memmove(data, &data[2*ifirst], 2*(n_sampled_in-ifirst)*sizeof(float)); |
| 839 |
if ((ifirst+n_sampled_out) >= n_sampled_in) |
| 840 |
memset(&data[2*(n_sampled_in-ifirst)], 0, 2*(n_sampled_out+ifirst-n_sampled_in)*sizeof(float)); |
| 841 |
|
| 842 |
if (tsout) |
| 843 |
{ |
| 844 |
// code to flip m on this output separately. in the present code this has already been accomplished by setting flipm. |
| 845 |
/* |
| 846 |
if (mfliparr[TSOUT] != mflipin) |
| 847 |
{ |
| 848 |
for (i = 0; i < n_sampled_out; i++) |
| 849 |
{ |
| 850 |
datatemp[2*i]=data[2*i]; |
| 851 |
datatemp[2*i+1]=-data[2*i+1]; |
| 852 |
} |
| 853 |
tmpptr=datatemp; |
| 854 |
} |
| 855 |
else |
| 856 |
tmpptr=data; |
| 857 |
*/ |
| 858 |
tsstart[1]=m; |
| 859 |
tsend[1]=m; |
| 860 |
outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status); |
| 861 |
outarr->bzero=segoutarr[TSOUT]->bzero; |
| 862 |
outarr->bscale=segoutarr[TSOUT]->bscale; |
| 863 |
status=drms_segment_writeslice_ext(segoutarr[TSOUT], outarr, tsstart, tsend, tstotal, 0); |
| 864 |
free(outarr); |
| 865 |
if (status != DRMS_SUCCESS) |
| 866 |
{ |
| 867 |
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]); |
| 868 |
return 0; |
| 869 |
} |
| 870 |
} |
| 871 |
|
| 872 |
if (fftout || fft1out || powout || mavgout) |
| 873 |
fftwf_execute(plan); |
| 874 |
|
| 875 |
if (fftout) |
| 876 |
{ |
| 877 |
// code to flip m on this output separately. in the present code this has already been accomplished by setting flipm. |
| 878 |
/* |
| 879 |
if (mfliparr[FFTOUT] != mflipin) |
| 880 |
{ |
| 881 |
datatemp[0]=data[0]; |
| 882 |
datatemp[1]=-data[1]; |
| 883 |
for (i = 1; i < n_sampled_out; i++) |
| 884 |
{ |
| 885 |
datatemp[2*i]=data[2*(n_sampled_out-i)]; |
| 886 |
datatemp[2*i+1]=-data[2*(n_sampled_out-i)+1]; |
| 887 |
} |
| 888 |
tmpptr=datatemp; |
| 889 |
} |
| 890 |
else |
| 891 |
tmpptr=data; |
| 892 |
*/ |
| 893 |
tsstart[1]=m; |
| 894 |
tsend[1]=m; |
| 895 |
outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status); |
| 896 |
outarr->bzero=segoutarr[FFTOUT]->bzero; |
| 897 |
outarr->bscale=segoutarr[FFTOUT]->bscale; |
| 898 |
status=drms_segment_writeslice_ext(segoutarr[FFTOUT], outarr, tsstart, tsend, tstotal, 0); |
| 899 |
free(outarr); |
| 900 |
if (status != DRMS_SUCCESS) |
| 901 |
{ |
| 902 |
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]); |
| 903 |
return 0; |
| 904 |
} |
| 905 |
} |
| 906 |
|
| 907 |
if (fft1out) |
| 908 |
{ |
| 909 |
|
| 910 |
fft1start[1]=lmode+m*msignarr[FFT1OUT]; |
| 911 |
fft1end[1]=lmode+m*msignarr[FFT1OUT]; |
| 912 |
outarr = drms_array_create(usetype, 2, fft1length, data+2*imin, &status); |
| 913 |
outarr->bzero=segoutarr[FFT1OUT]->bzero; |
| 914 |
outarr->bscale=segoutarr[FFT1OUT]->bscale; |
| 915 |
status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0); |
| 916 |
free(outarr); |
| 917 |
if (status != DRMS_SUCCESS) |
| 918 |
{ |
| 919 |
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]); |
| 920 |
return 0; |
| 921 |
} |
| 922 |
|
| 923 |
if (m > 0) |
| 924 |
{ |
| 925 |
/* Do negative m */ |
| 926 |
for (i=0; i<npow; i++) |
| 927 |
{ |
| 928 |
/* Conjugate for negative m */ |
| 929 |
datahold[2*i] = data[2*(n_sampled_out-imin-i)]; |
| 930 |
datahold[2*i+1] = -data[2*(n_sampled_out-imin-i)+1]; |
| 931 |
} |
| 932 |
fft1start[1]=lmode-m*msignarr[FFT1OUT]; |
| 933 |
fft1end[1]=lmode-m*msignarr[FFT1OUT]; |
| 934 |
outarr = drms_array_create(usetype, 2, fft1length, datahold, &status); |
| 935 |
outarr->bzero=segoutarr[FFT1OUT]->bzero; |
| 936 |
outarr->bscale=segoutarr[FFT1OUT]->bscale; |
| 937 |
status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0); |
| 938 |
free(outarr); |
| 939 |
if (status != DRMS_SUCCESS) |
| 940 |
{ |
| 941 |
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]); |
| 942 |
return 0; |
| 943 |
} |
| 944 |
} |
| 945 |
} |
| 946 |
|
| 947 |
if (powout || mavgout) |
| 948 |
for (i = 0; i < n_sampled_out; i++) |
| 949 |
pow[i] = data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]; |
| 950 |
|
| 951 |
if (powout) |
| 952 |
{ |
| 953 |
powstart[1]=lmode+m*msignarr[POWOUT]; |
| 954 |
powend[1]=lmode+m*msignarr[POWOUT]; |
| 955 |
outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status); |
| 956 |
outarr->bzero=segoutarr[POWOUT]->bzero; |
| 957 |
outarr->bscale=segoutarr[POWOUT]->bscale; |
| 958 |
status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0); |
| 959 |
free(outarr); |
| 960 |
if (status != DRMS_SUCCESS) |
| 961 |
{ |
| 962 |
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]); |
| 963 |
return 0; |
| 964 |
} |
| 965 |
|
| 966 |
if (m > 0) |
| 967 |
{ |
| 968 |
/* Do negative m */ |
| 969 |
if (imin) |
| 970 |
datahold[0] = pow[n_sampled_out-imin]; |
| 971 |
else |
| 972 |
datahold[0] = pow[0]; |
| 973 |
for (i=1; i<npow;i++) |
| 974 |
datahold[i] = pow[n_sampled_out-imin-i]; |
| 975 |
|
| 976 |
powstart[1]=lmode-m*msignarr[POWOUT]; |
| 977 |
powend[1]=lmode-m*msignarr[POWOUT]; |
| 978 |
outarr = drms_array_create(usetype, 2, powlength, datahold, &status); |
| 979 |
outarr->bzero=segoutarr[POWOUT]->bzero; |
| 980 |
outarr->bscale=segoutarr[POWOUT]->bscale; |
| 981 |
status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0); |
| 982 |
free(outarr); |
| 983 |
if (status != DRMS_SUCCESS) |
| 984 |
{ |
| 985 |
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]); |
| 986 |
return 0; |
| 987 |
} |
| 988 |
} |
| 989 |
} |
| 990 |
|
| 991 |
if (mavgout) |
| 992 |
{ |
| 993 |
mshift = floor(df1*splitting(lmode,m*msignarr[MAVGOUT])+0.5f); |
| 994 |
if (mshift >= 0) |
| 995 |
for (i=0; i<n_sampled_out/2-mshift; i++) |
| 996 |
msum[i] += pow[mshift+i]; |
| 997 |
else |
| 998 |
for (i=0; i<n_sampled_out/2+mshift; i++) |
| 999 |
msum[i-mshift] += pow[i]; |
| 1000 |
if (m > 0) |
| 1001 |
{ |
| 1002 |
/* Do negative m */ |
| 1003 |
mshift = floor(df1*splitting(lmode,-m*msignarr[MAVGOUT])+0.5f); |
| 1004 |
if (mshift >=0) |
| 1005 |
for (i=1; i<n_sampled_out/2-mshift;i++) |
| 1006 |
msum[i] += pow[n_sampled_out-mshift-i]; |
| 1007 |
else |
| 1008 |
for (i=1; i<n_sampled_out/2+mshift; i++) |
| 1009 |
msum[i-mshift] += pow[n_sampled_out-i]; |
| 1010 |
} |
| 1011 |
} |
| 1012 |
|
| 1013 |
if (arpowout) |
| 1014 |
{ |
| 1015 |
memmove(data, ar_coeff, 2*(ar_order+1)*sizeof(float)); |
| 1016 |
memset(&data[2*(ar_order+1)],0,2*(n_sampled_out-ar_order-1)*sizeof(float)); |
| 1017 |
fftwf_execute(plan); |
| 1018 |
for (i = 0; i < n_sampled_out; i++) |
| 1019 |
pow[i] = 1.0/(data[2*i]*data[2*i] + data[2*i+1]*data[2*i+1]); |
| 1020 |
|
| 1021 |
powstart[1]=lmode+m*msignarr[ARPOWOUT]; |
| 1022 |
powend[1]=lmode+m*msignarr[ARPOWOUT]; |
| 1023 |
outarr = drms_array_create(usetype, 2, powlength, pow+imin, &status); |
| 1024 |
outarr->bzero=segoutarr[ARPOWOUT]->bzero; |
| 1025 |
outarr->bscale=segoutarr[ARPOWOUT]->bscale; |
| 1026 |
status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0); |
| 1027 |
free(outarr); |
| 1028 |
if (status != DRMS_SUCCESS) |
| 1029 |
{ |
| 1030 |
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]); |
| 1031 |
return 0; |
| 1032 |
} |
| 1033 |
|
| 1034 |
if (m > 0) |
| 1035 |
{ |
| 1036 |
/* Do negative m */ |
| 1037 |
if (imin) |
| 1038 |
datahold[0] = pow[n_sampled_out-imin]; |
| 1039 |
else |
| 1040 |
datahold[0] = pow[0]; |
| 1041 |
for (i=1; i<npow;i++) |
| 1042 |
datahold[i] = pow[n_sampled_out-imin-i]; |
| 1043 |
|
| 1044 |
powstart[1]=lmode-m*msignarr[ARPOWOUT]; |
| 1045 |
powend[1]=lmode-m*msignarr[ARPOWOUT]; |
| 1046 |
outarr = drms_array_create(usetype, 2, powlength, datahold, &status); |
| 1047 |
outarr->bzero=segoutarr[ARPOWOUT]->bzero; |
| 1048 |
outarr->bscale=segoutarr[ARPOWOUT]->bscale; |
| 1049 |
status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0); |
| 1050 |
free(outarr); |
| 1051 |
if (status != DRMS_SUCCESS) |
| 1052 |
{ |
| 1053 |
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]); |
| 1054 |
return 0; |
| 1055 |
} |
| 1056 |
} |
| 1057 |
} |
| 1058 |
|
| 1059 |
} //end of loop over m |
| 1060 |
|
| 1061 |
if (mavgout) |
| 1062 |
{ |
| 1063 |
c=1.0f/(2*lmode+1); |
| 1064 |
for (i=0; i<npow; i++) |
| 1065 |
datahold[i] = msum[imin+i] * c; |
| 1066 |
outarr = drms_array_create(usetype, 2, powlength, datahold, &status); |
| 1067 |
outarr->bzero=segoutarr[MAVGOUT]->bzero; |
| 1068 |
outarr->bscale=segoutarr[MAVGOUT]->bscale; |
| 1069 |
status=drms_segment_write(segoutarr[MAVGOUT], outarr, 0); |
| 1070 |
free(outarr); |
| 1071 |
if (status != DRMS_SUCCESS) |
| 1072 |
{ |
| 1073 |
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); |
| 1074 |
return 0; |
| 1075 |
} |
| 1076 |
} |
| 1077 |
|
| 1078 |
for (is=0;is<6;is++) |
| 1079 |
{ |
| 1080 |
if (!flagarr[is]) |
| 1081 |
continue; |
| 1082 |
|
| 1083 |
outrec=outrecsetarr[is]->records[irec]; |
| 1084 |
drms_copykeys(outrec, inrec, 0, kDRMS_KeyClass_Explicit); |
| 1085 |
drms_setkey_int(outrec, "QUALITY", 0); |
| 1086 |
drms_setkey_time(outrec, "T_START", tstartout); |
| 1087 |
drms_setkey_time(outrec, "T_STOP", tstopout); |
| 1088 |
drms_setkey_time(outrec, "T_OBS", (tstartout + tstopout)/2); |
| 1089 |
drms_setkey_float(outrec, "T_STEP", tsamplex); |
| 1090 |
drms_setkey_int(outrec, "NDT", n_sampled_out); |
| 1091 |
if (tagflag) |
| 1092 |
drms_setkey_string(outrec, "TAG", tag); |
| 1093 |
if (verflag) |
| 1094 |
drms_setkey_string(outrec, "VERSION", version); |
| 1095 |
|
| 1096 |
tnow = (double)time(NULL); |
| 1097 |
tnow += UNIX_epoch; |
| 1098 |
drms_setkey_time(outrec, "DATE", tnow); |
| 1099 |
|
| 1100 |
// drms_close_record(outrec, DRMS_INSERT_RECORD); |
| 1101 |
|
| 1102 |
} |
| 1103 |
|
| 1104 |
} |
| 1105 |
|
| 1106 |
drms_close_records(inrecset, DRMS_FREE_RECORD); |
| 1107 |
for (is=0;is<6;is++) |
| 1108 |
{ |
| 1109 |
if (!flagarr[is]) |
| 1110 |
continue; |
| 1111 |
drms_close_records(outrecsetarr[is], DRMS_INSERT_RECORD); |
| 1112 |
} |
| 1113 |
|
| 1114 |
if(ar_coeff != NULL) |
| 1115 |
free(ar_coeff); |
| 1116 |
|
| 1117 |
if (fftout || fft1out || powout || arpowout || mavgout) |
| 1118 |
fftwf_destroy_plan(plan); |
| 1119 |
|
| 1120 |
wt=getwalltime(); |
| 1121 |
ct=getcputime(&ut, &st); |
| 1122 |
if (verbflag) |
| 1123 |
{ |
| 1124 |
printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n", |
| 1125 |
wt-wt0, ct-ct0); |
| 1126 |
} |
| 1127 |
|
| 1128 |
printf("module %s successful completion\n", cmdparams.argv[0]); |
| 1129 |
return 0; |
| 1130 |
} |
| 1131 |
|
| 1132 |
|
| 1133 |
static double splitting(int l, int m) |
| 1134 |
{ |
| 1135 |
double a1x[] = {406.0,407.0,408.0,410.5,412.0,411.5,409.0,406.0,406.0}; |
| 1136 |
double a3x[] = {10.0,10.0,15.0,19.5,21.0,21.3,21.3,21.3,21.8}; |
| 1137 |
double a5x[] = {0.0,0.0,0.0,-1.5,-2.5,-3.5,-4.0,-4.0,-4.5}; |
| 1138 |
/* f0 is the frequency used for generating the above a-coeff */ |
| 1139 |
double f0 = 1500.0; |
| 1140 |
/* fcol is the frequency for which to optimize the collaps */ |
| 1141 |
double fcol = 800.0; |
| 1142 |
|
| 1143 |
double ll,x,x3,x5,a1,a2,a3,a5,frac,lx; |
| 1144 |
int ix; |
| 1145 |
|
| 1146 |
if (l == 0) |
| 1147 |
ll = 1; |
| 1148 |
else |
| 1149 |
ll = sqrt(l*(l+1.)); |
| 1150 |
x = m/ll; |
| 1151 |
x3 = x*x*x; |
| 1152 |
x5 = x3*x*x; |
| 1153 |
lx = 5*log10(l*f0/fcol)-4; |
| 1154 |
ix = floor(lx); |
| 1155 |
frac = lx-ix; |
| 1156 |
if (lx <= 0) { |
| 1157 |
ix = 0; |
| 1158 |
frac = 0.0; |
| 1159 |
} |
| 1160 |
if (lx >=8) { |
| 1161 |
ix = 7; |
| 1162 |
frac = 1.0; |
| 1163 |
} |
| 1164 |
a1 = (1.0-frac)*a1x[ix]+frac*a1x[ix+1]; |
| 1165 |
a2 = 0.0; |
| 1166 |
a3 = (1.0-frac)*a3x[ix]+frac*a3x[ix+1]; |
| 1167 |
a5 = (1.0-frac)*a5x[ix]+frac*a5x[ix+1]; |
| 1168 |
|
| 1169 |
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; |
| 1170 |
} |
| 1171 |
|
| 1172 |
|
| 1173 |
|
| 1174 |
/* Extract the data samples actually used by skipping or averaging |
| 1175 |
data. Replace missing data that are not marked as gaps with zero. */ |
| 1176 |
void extract_data(int n_in, int *gaps, float *data_in, float *data_out) |
| 1177 |
{ |
| 1178 |
int i,j, nmiss = 0; |
| 1179 |
// the following check will be skipped in production executables since they will be built with NDEBUG defined. |
| 1180 |
// it doesn't matter because the check is already done in DoIt(). |
| 1181 |
assert(nskip==1 || navg==1); |
| 1182 |
if ((navg == 1) && (nskip == 1)) |
| 1183 |
{ |
| 1184 |
for (i=0; i<n_in; i++) |
| 1185 |
{ |
| 1186 |
if (gaps[i]==0 ) |
| 1187 |
{ |
| 1188 |
data_out[2*i] = 0.0f; |
| 1189 |
data_out[2*i+1] = 0.0f; |
| 1190 |
} |
| 1191 |
else |
| 1192 |
{ |
| 1193 |
if (isnan(data_in[2*i]) || isnan(data_in[2*i+1])) |
| 1194 |
{ |
| 1195 |
data_out[2*i] = 0.0f; |
| 1196 |
data_out[2*i+1] = 0.0f; |
| 1197 |
gaps[i] = 0; |
| 1198 |
nmiss++; |
| 1199 |
} |
| 1200 |
else |
| 1201 |
{ |
| 1202 |
data_out[2*i] = data_in[2*i]; |
| 1203 |
data_out[2*i+1] = flipm ? -data_in[2*i+1] : data_in[2*i+1]; |
| 1204 |
} |
| 1205 |
} |
| 1206 |
} |
| 1207 |
} |
| 1208 |
else if (nskip != 1) |
| 1209 |
{ |
| 1210 |
for (i = 0; i<n_in; i++) |
| 1211 |
{ |
| 1212 |
if (gaps[i]==0 ) |
| 1213 |
{ |
| 1214 |
data_out[2*i] = 0.0f; |
| 1215 |
data_out[2*i+1] = 0.0f; |
| 1216 |
} |
| 1217 |
else |
| 1218 |
{ |
| 1219 |
int ix = nskip*i+noff; |
| 1220 |
if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1])) |
| 1221 |
{ |
| 1222 |
data_out[2*i] = 0.0f; |
| 1223 |
data_out[2*i+1] = 0.0f; |
| 1224 |
gaps[i] = 0; |
| 1225 |
nmiss++; |
| 1226 |
} |
| 1227 |
else |
| 1228 |
{ |
| 1229 |
data_out[2*i] = data_in[2*ix]; |
| 1230 |
data_out[2*i+1] = flipm ? -data_in[2*ix+1] : data_in[2*ix+1]; |
| 1231 |
} |
| 1232 |
} |
| 1233 |
} |
| 1234 |
} |
| 1235 |
else if (navg != 1) |
| 1236 |
{ |
| 1237 |
for (i = 0; i<n_in; i++) |
| 1238 |
{ |
| 1239 |
if (gaps[i]==0 ) |
| 1240 |
{ |
| 1241 |
data_out[2*i] = 0.0f; |
| 1242 |
data_out[2*i+1] = 0.0f; |
| 1243 |
} |
| 1244 |
else |
| 1245 |
{ |
| 1246 |
float avgr = 0.0f; |
| 1247 |
float avgi = 0.0f; |
| 1248 |
for (j=0; j<navg; j++) |
| 1249 |
{ |
| 1250 |
int ix = navg*i+j; |
| 1251 |
if (isnan(data_in[2*ix]) || isnan(data_in[2*ix+1])) |
| 1252 |
{ |
| 1253 |
data_out[2*i] = 0.0f; |
| 1254 |
data_out[2*i+1] = 0.0f; |
| 1255 |
gaps[i] = 0; |
| 1256 |
nmiss++; |
| 1257 |
break; |
| 1258 |
} |
| 1259 |
else |
| 1260 |
{ |
| 1261 |
avgr += data_in[2*ix]; |
| 1262 |
avgi += data_in[2*ix+1]; |
| 1263 |
} |
| 1264 |
} |
| 1265 |
if (j == navg) |
| 1266 |
{ |
| 1267 |
data_out[2*i] = avgr/navg; |
| 1268 |
data_out[2*i+1] = avgi/navg; |
| 1269 |
} |
| 1270 |
} |
| 1271 |
} |
| 1272 |
} |
| 1273 |
} |
| 1274 |
|
| 1275 |
/* Extract the windows function for the samples actually used after |
| 1276 |
subsampling or averaging. Modify the list of continous sections |
| 1277 |
accordingly, and make sure we do not average across section boundaries. */ |
| 1278 |
static int extract_gaps(int n_in, int *gaps_in, int *n_out, int *gaps_out, |
| 1279 |
float tsample_in, float *tsample_out, |
| 1280 |
int *num_sections, int *sections) |
| 1281 |
{ |
| 1282 |
// the following check will be skipped in production executables since they will be built with NDEBUG defined. |
| 1283 |
// it doesn't matter because the check is already done in DoIt(). |
| 1284 |
assert(nskip==1 || navg==1); |
| 1285 |
int i,j,sect, start,stop; |
| 1286 |
|
| 1287 |
|
| 1288 |
if (*num_sections<1) |
| 1289 |
{ |
| 1290 |
fprintf(stderr,"Apparantly no sections of data available."); |
| 1291 |
return 1; |
| 1292 |
} |
| 1293 |
/* Mask out data that in between sections if this hasn't already been |
| 1294 |
done in gapfile. */ |
| 1295 |
for (i=0; i<sections[0] && i<n_in; i++) |
| 1296 |
gaps_in[i] = 0; |
| 1297 |
for (sect=1; sect<(*num_sections); sect++) |
| 1298 |
{ |
| 1299 |
for (i=sections[2*sect-1]+1; i<sections[2*sect] && i<n_in; i++) |
| 1300 |
gaps_in[i] = 0; |
| 1301 |
} |
| 1302 |
for (i=sections[2*(*num_sections-1)+1]+1; i<n_in; i++) |
| 1303 |
gaps_in[i] = 0; |
| 1304 |
|
| 1305 |
|
| 1306 |
if ((navg == 1) && (nskip == 1)) |
| 1307 |
{ |
| 1308 |
*n_out = n_in; |
| 1309 |
*tsample_out = tsample_in; |
| 1310 |
for (i=0; i<*n_out; i++) |
| 1311 |
gaps_out[i] = gaps_in[i]; |
| 1312 |
} |
| 1313 |
else if (nskip != 1) |
| 1314 |
{ |
| 1315 |
*n_out = n_in/nskip; |
| 1316 |
*tsample_out = nskip*tsample_in; |
| 1317 |
for (i=0; i<*n_out; i++) |
| 1318 |
{ |
| 1319 |
gaps_out[i] = gaps_in[nskip*i+noff]; |
| 1320 |
} |
| 1321 |
for (i=0; i<2*(*num_sections); i+=2) |
| 1322 |
{ |
| 1323 |
start = (int)ceilf(((float)(sections[i]-noff))/nskip); |
| 1324 |
stop = (int)floorf(((float)(sections[i+1]-noff))/nskip); |
| 1325 |
if (start <= stop) |
| 1326 |
{ |
| 1327 |
sections[i] = start; |
| 1328 |
sections[i+1] = stop; |
| 1329 |
} |
| 1330 |
else |
| 1331 |
{ |
| 1332 |
/* This section was skipped entirely. */ |
| 1333 |
for (j=i; j< 2*(*num_sections-1); j+=2) |
| 1334 |
{ |
| 1335 |
sections[j] = sections[j+2]; |
| 1336 |
sections[j+1] = sections[j+3]; |
| 1337 |
} |
| 1338 |
*num_sections -= 1; |
| 1339 |
} |
| 1340 |
} |
| 1341 |
} |
| 1342 |
else if (navg != 1) |
| 1343 |
{ |
| 1344 |
sect = 0; |
| 1345 |
*n_out = n_in/navg; |
| 1346 |
*tsample_out = tsample*navg; |
| 1347 |
for (i=0; i<*n_out; i++) |
| 1348 |
{ |
| 1349 |
int igx = 1; |
| 1350 |
while (sect < *num_sections && |
| 1351 |
!(navg*i >= sections[2*sect] && navg*i >= sections[2*sect+1])) |
| 1352 |
sect++; |
| 1353 |
|
| 1354 |
/* Don't allow averaging of data from different sections. */ |
| 1355 |
if (navg*i >= sections[2*sect] && (navg*i+navg-1)<=sections[2*sect+1]) |
| 1356 |
{ |
| 1357 |
for (j=0; j<navg; j++) |
| 1358 |
igx *= gaps_in[navg*i+j]; |
| 1359 |
gaps_out[i] = igx; |
| 1360 |
} |
| 1361 |
else |
| 1362 |
gaps_out[i] = 0; |
| 1363 |
} |
| 1364 |
for (i=0; i<2*(*num_sections); i+=2) |
| 1365 |
{ |
| 1366 |
start = (int)ceilf(((float)sections[i])/navg); |
| 1367 |
stop = (int)floorf(((float)sections[i+1])/navg); |
| 1368 |
if (start <= stop) |
| 1369 |
{ |
| 1370 |
sections[i] = start; |
| 1371 |
sections[i+1] = stop; |
| 1372 |
} |
| 1373 |
else |
| 1374 |
{ |
| 1375 |
/* This section was skipped entirely. */ |
| 1376 |
for (j=i; j< 2*(*num_sections-1); j+=2) |
| 1377 |
{ |
| 1378 |
sections[j] = sections[j+2]; |
| 1379 |
sections[j+1] = sections[j+3]; |
| 1380 |
} |
| 1381 |
*num_sections -= 1; |
| 1382 |
} |
| 1383 |
} |
| 1384 |
} |
| 1385 |
return 0; |
| 1386 |
} |
| 1387 |
|
| 1388 |
|
| 1389 |
/* |
| 1390 |
Read a list of sections [start_sample:end_sample] that should |
| 1391 |
contain continuous data. When detrending, jumps are allow to |
| 1392 |
occur between sections. The sections file should be a text file |
| 1393 |
of the form: |
| 1394 |
|
| 1395 |
num |
| 1396 |
start_1 end_1 |
| 1397 |
start_2 end_2 |
| 1398 |
... |
| 1399 |
start_num end_num |
| 1400 |
*/ |
| 1401 |
static int read_section_file(char *filename, int *num_sections, int **sections) |
| 1402 |
{ |
| 1403 |
FILE *fh; |
| 1404 |
int i; |
| 1405 |
|
| 1406 |
fh = fopen(filename,"r"); |
| 1407 |
if (fh!=NULL) |
| 1408 |
{ |
| 1409 |
fscanf(fh,"%d",num_sections); |
| 1410 |
*sections = (int *)malloc(2*(*num_sections)*sizeof(int)); |
| 1411 |
|
| 1412 |
for (i=0 ;i<2*(*num_sections); i+=2) |
| 1413 |
{ |
| 1414 |
fscanf(fh,"%d",&(*sections)[i]); |
| 1415 |
fscanf(fh,"%d",&(*sections)[i+1]); |
| 1416 |
|
| 1417 |
if ((*sections)[i]>(*sections)[i+1] || |
| 1418 |
(i>0 && (*sections)[i]<=(*sections)[i-1])) |
| 1419 |
{ |
| 1420 |
fprintf(stderr,"Invalid sections file, sections overlap.\n"); |
| 1421 |
return 2; |
| 1422 |
} |
| 1423 |
} |
| 1424 |
fclose(fh); |
| 1425 |
return 0; |
| 1426 |
} |
| 1427 |
else |
| 1428 |
return 1; |
| 1429 |
} |
| 1430 |
|
