1 tplarson 1.1 /* this JSOC module is a port of an SOI module written by Rasmus Larsen.
2 * ported by Tim Larson.
|
3 tplarson 1.3 * tim doubts this works correctly with nskip or navg
|
4 tplarson 1.1 */
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 kehcheng 1.7 #include "fitsio.h"
|
23 tplarson 1.1
24 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
25 #define kNOTSPECIFIED "not specified"
26
27 char *module_name = "jtsfiddle";
|
28 tplarson 1.15 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 tplarson 1.1
30 /* Command line arguments: */
31 ModuleArgs_t module_args[] =
32 {
33 {ARG_STRING, "in", "", "input data records"},
|
34 tplarson 1.3 {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 tplarson 1.1 {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 tplarson 1.13 {ARG_STRING, "TAG", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, usually used as an additional primekey"},
|
48 tplarson 1.15 {ARG_STRING, "VERSION", kNOTSPECIFIED, "this parameter sets a keyword of the same name to the same value, useful for selecting obsolete records"},
|
49 tplarson 1.13
|
50 tplarson 1.1 {ARG_STRING, "LOGFILE", "stdout", ""},
|
51 tplarson 1.2 {ARG_STRING, "GAPFILE", "", "record or file containing window function"},
52 {ARG_STRING, "SECTIONFILE", kNOTSPECIFIED, "record or file specifying continuous sections of data"},
53
|
54 tplarson 1.1 {ARG_INT, "NDT", "-1", "", ""},
|
55 tplarson 1.2 // {ARG_FLOAT, "TSAMPLE", "60", "", ""},
|
56 tplarson 1.1 {ARG_INT, "IFILL", "1", "", ""},
57 {ARG_INT, "IFIRST", "0", "", ""},
|
58 tplarson 1.3 // {ARG_INT, "FLIPM", "0", "", ""},
|
59 tplarson 1.1 {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 tplarson 1.3 {ARG_INT, "DETREND_FIRST", "0", "", ""},
67 /*
|
68 tplarson 1.1 {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 tplarson 1.3 */
|
75 tplarson 1.1 {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 tplarson 1.3 #include "set_history.c"
|
86 tplarson 1.1
|
87 tplarson 1.12 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 tplarson 1.14 //char *getdetrendversion(void);
96 //char *getgapfillversion(void);
|
97 tplarson 1.1
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 tplarson 1.3 static int detrend_length, detrend_skip, detrend_first;
|
104 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.3 float *datarr, *datptr, *outptr;
138 float *datahold;
139 char tstartstr[100];
140 int lmin, lmax;
|
141 tplarson 1.1
142 int fstatus = 0;
143 fitsfile *fitsptr;
144 long fdims[1];
145 long fpixel[]={1};
146 int *anynul=0;
147
|
148 tplarson 1.3 int instart[2], inend[2];
|
149 tplarson 1.10 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 tplarson 1.1 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 tplarson 1.2 DRMS_Segment_t *gapseg = NULL;
165 DRMS_Array_t *gaparr = NULL;
|
166 tplarson 1.1
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 tplarson 1.6 double tstart, tstartsave, tstop, tstopsave, cadence, tstartout, tstopout;
|
175 tplarson 1.1 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 tplarson 1.3 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 tplarson 1.1 n_sampled_out = cmdparams_save_int (&cmdparams, "NDT", &newstat);
|
188 tplarson 1.2 // tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat);
|
189 tplarson 1.1 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 tplarson 1.3 // flipm = cmdparams_save_int (&cmdparams, "FLIPM", &newstat);
|
196 tplarson 1.1 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 tplarson 1.3 detrend_first = cmdparams_save_int (&cmdparams, "DETREND_FIRST", &newstat);
200 /*
|
201 tplarson 1.1 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 tplarson 1.3 */
|
208 tplarson 1.1 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 tplarson 1.3 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 tplarson 1.1 int seginflag = strcmp(kNOTSPECIFIED, segnamein);
219 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
|
220 tplarson 1.3 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
221 char *srclinkname = (char *)cmdparams_save_str(&cmdparams, "srclink", &newstat);
|
222 tplarson 1.1 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 tplarson 1.13 char *tag = (char *)cmdparams_save_str(&cmdparams, "TAG", &newstat);
230 int tagflag = strcmp(kNOTSPECIFIED, tag);
|
231 tplarson 1.15 char *version = (char *)cmdparams_save_str(&cmdparams, "VERSION", &newstat);
232 int verflag = strcmp(kNOTSPECIFIED, version);
|
233 tplarson 1.13
|
234 tplarson 1.3 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 tplarson 1.3 serieslist[POWOUT]=powseries;
256 serieslist[ARPOWOUT]=arpowseries;
257 serieslist[MAVGOUT]=mavgseries;
258 long long histrecnumarr[6]={-1, -1, -1, -1, -1, -1};
|
259 tplarson 1.5 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 tplarson 1.1
263 /**** Sanity check of input parameters. ****/
264
265 if ((tsout == 0) && (fftout == 0) && (fft1out == 0) && (powout == 0) && (mavgout == 0) && (arpowout == 0))
266 {
|
267 tplarson 1.3 fprintf(stderr, "ERROR: must specify an output dataseries\n");
268 return 1;
|
269 tplarson 1.1 }
|
270 tplarson 1.3
|
271 tplarson 1.1 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 tplarson 1.3 return 1;
|
279 tplarson 1.1 }
280 if (noff < 0 || noff >= nskip)
281 {
282 fprintf(stderr, "ERROR: noff must be >= 0 and < nskip\n");
|
283 tplarson 1.3 return 1;
|
284 tplarson 1.1 }
|
285 tplarson 1.3 if (arpowout && !ifill)
|
286 tplarson 1.1 {
287 fprintf(stderr, "ERROR: must have nonzero ifill with arpowout\n");
|
288 tplarson 1.3 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 tplarson 1.1 }
307
308
|
309 tplarson 1.3 // 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 tplarson 1.5 char *outchecklist[] = {"T_START", "QUALITY", "LMIN", "LMAX", "NDT"};
|
312 tplarson 1.11 int is, ishold, itest, mflip;
|
313 tplarson 1.3 char *holdseries="";
|
314 tplarson 1.14 /*
|
315 tplarson 1.11 char *cvsinfo;
316 cvsinfo = (char *)malloc(1024);
|
317 tplarson 1.15 strcpy(cvsinfo,"$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jtsfiddle.c,v 1.14 2013/04/28 08:05:21 tplarson Exp $");
|
318 tplarson 1.11 strcat(cvsinfo,"\n");
319 strcat(cvsinfo,getdetrendversion());
320 strcat(cvsinfo,"\n");
321 strcat(cvsinfo,getgapfillversion());
|
322 tplarson 1.14 */
|
323 tplarson 1.3 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 tplarson 1.1
|
342 tplarson 1.3
343 // set up ancillary dataseries for processing metadata
344 if (histlink != NULL)
345 {
346 if (strcmp(holdseries, histlink->info->target_series))
347 {
|
348 tplarson 1.11 ishold=is;
|
349 tplarson 1.3 holdseries=strdup(histlink->info->target_series);
|
350 tplarson 1.14 histrecnumarr[is]=set_history(histlink);
|
351 tplarson 1.3 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 tplarson 1.11 histrecnumarr[is]=histrecnumarr[ishold];
|
361 tplarson 1.3 }
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 tplarson 1.2 {
|
370 tplarson 1.3 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 tplarson 1.2 }
|
378 tplarson 1.1
|
379 tplarson 1.3 mflip = drms_getkey_int(tempoutrec, "MFLIPPED", &status);
380 if (status == DRMS_SUCCESS)
381 mfliparr[is]=mflip;
|
382 tplarson 1.1
|
383 tplarson 1.3 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
384 }
|
385 tplarson 1.1
|
386 tplarson 1.3 // 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 tplarson 1.1
|
393 tplarson 1.2 if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus))
|
394 tplarson 1.1 {
|
395 tplarson 1.2 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 tplarson 1.3 return 1;
|
400 tplarson 1.2 }
401 if (gaprecset->n == 0)
402 {
403 fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
|
404 tplarson 1.3 return 1;
|
405 tplarson 1.2 }
406 if (gaprecset->n > 1)
407 {
408 fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
|
409 tplarson 1.3 return 1;
|
410 tplarson 1.2 }
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 tplarson 1.3 return 1;
|
418 tplarson 1.2 }
419 agaps=(int *)(gaparr->data);
420 gapsize=gaparr->axis[0];
|
421 tplarson 1.1 }
|
422 tplarson 1.2 else
|
423 tplarson 1.1 {
|
424 tplarson 1.2 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 tplarson 1.3 return 1;
|
433 tplarson 1.2 }
|
434 tplarson 1.1 }
435
|
436 tplarson 1.9 if (verbflag)
437 printf("gapfile read, gapsize = %d\n",gapsize);
|
438 tplarson 1.1
439 data_dim=gapsize;
440 if (n_sampled_out>gapsize)
441 data_dim=n_sampled_out;
442
|
443 tplarson 1.3 /* 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 tplarson 1.1 /* Read location of jumps. */
|
453 tplarson 1.2 if (!strcmp(sectionfile,kNOTSPECIFIED))
|
454 tplarson 1.1 {
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 tplarson 1.2 int secstat;
465 if (secstat=read_section_file(sectionfile, &num_sections, §ions))
|
466 tplarson 1.1 {
|
467 tplarson 1.2 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 tplarson 1.3 return 1;
|
472 tplarson 1.2 }
473 if (secrecset->n == 0)
474 {
475 fprintf(stderr,"ERROR: sectionfile recordset contains no records\n");
|
476 tplarson 1.3 return 1;
|
477 tplarson 1.2 }
478 if (secrecset->n > 1)
479 {
480 fprintf(stderr,"ERROR: sectionfile recordset with more than 1 record not yet supported.\n");
|
481 tplarson 1.3 return 1;
|
482 tplarson 1.2 }
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 tplarson 1.3 return 1;
|
488 tplarson 1.2 }
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 tplarson 1.3 return 1;
|
503 tplarson 1.2 }
504 }
505 free(sectkey);
|
506 tplarson 1.1 }
507 }
508
|
509 tplarson 1.9 if (verbflag)
510 printf("num_sections = %d\n",num_sections);
|
511 tplarson 1.1
|
512 tplarson 1.3 DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status);
|
513 tplarson 1.2
|
514 tplarson 1.3 if (status != DRMS_SUCCESS || inrecset == NULL)
|
515 tplarson 1.2 {
516 fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
|
517 tplarson 1.3 return 1;
|
518 tplarson 1.2 }
|
519 tplarson 1.3 if (inrecset->n == 0)
|
520 tplarson 1.2 {
|
521 tplarson 1.4 fprintf(stderr, "ERROR: input recordset contains no records\n");
|
522 tplarson 1.3 return 1;
|
523 tplarson 1.2 }
524
|
525 tplarson 1.3 inrec=inrecset->records[0];
526 char *inchecklist[] = {"T_START", "T_STOP", "QUALITY", "LMIN", "LMAX", "T_STEP"};
|
527 tplarson 1.2 for (itest=0; itest < ARRLENGTH(inchecklist); itest++)
528 {
529 DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
|
530 tplarson 1.3 if (inkeytest == NULL)
|
531 tplarson 1.2 {
532 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
|
533 tplarson 1.3 drms_close_records(inrecset, DRMS_FREE_RECORD);
534 return 1;
|
535 tplarson 1.2 }
536 }
537
|
538 tplarson 1.6 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 tplarson 1.2
|
542 tplarson 1.3 int mflipin=drms_getkey_int(inrec, "MFLIPPED", &status);
543 if (status != DRMS_SUCCESS) mflipin=0;
|
544 tplarson 1.2
|
545 tplarson 1.3 // 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 tplarson 1.2
|
553 tplarson 1.3 for (is=2;is<5;is++)
|
554 tplarson 1.1 {
|
555 tplarson 1.3 if (!flagarr[is])
556 continue;
557 if (mfliparr[is] != (mflipin || flipm))
558 msignarr[is]=-1;
559 else
560 msignarr[is]=1;
|
561 tplarson 1.1 }
562
|
563 tplarson 1.3 if (mflipin)
|
564 tplarson 1.8 msignarr[MAVGOUT]=1;
565 else
|
566 tplarson 1.3 msignarr[MAVGOUT]=-1;
|
567 tplarson 1.2
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 tplarson 1.1 {
571 fprintf(stderr, "ERROR: problem in extract_gaps\n");
|
572 tplarson 1.3 return 1;
|
573 tplarson 1.1 }
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 tplarson 1.3 }
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 tplarson 1.1
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 tplarson 1.2 /* Set sum to zero before starting */
600 if (mavgout)
|
601 tplarson 1.15 msum = (float *)malloc((n_sampled_out/2)*sizeof(float));
|
602 tplarson 1.2
|
603 tplarson 1.3 if (fft1out || powout || mavgout || arpowout)
|
604 tplarson 1.1 {
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 tplarson 1.3 datahold=(float *)malloc(2*npow*sizeof(float));
614 }
|
615 tplarson 1.1
|
616 tplarson 1.3 // 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 tplarson 1.10 tslength[0]=tstotal[0]=2*n_sampled_out;
|
623 tplarson 1.3 tslength[1]=1;
624 tsstart[0]=0;
625 tsend[0]=2*n_sampled_out-1;
626 }
627 if (fft1out)
628 {
|
629 tplarson 1.10 fft1length[0]=fft1total[0]=2*npow;
|
630 tplarson 1.3 fft1length[1]=1;
631 fft1start[0]=0;
632 fft1end[0]=2*npow-1;
633 }
634 if (powout || arpowout || mavgout)
635 {
|
636 tplarson 1.10 powlength[0]=powtotal[0]=npow;
|
637 tplarson 1.3 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 tplarson 1.5 for (is=0; is<6; is++)
652 {
653 if (!flagarr[is])
654 continue;
|
655 tplarson 1.3
|
656 tplarson 1.5 outrecsetarr[is] = drms_create_records(drms_env, inrecset->n, serieslist[is], DRMS_PERMANENT, &status);
|
657 tplarson 1.3
|
658 tplarson 1.5 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 tplarson 1.3 }
664
|
665 tplarson 1.5 int irec;
666 for (irec=0;irec<inrecset->n;irec++)
|
667 tplarson 1.3 {
668
|
669 tplarson 1.5 inrec = inrecset->records[irec];
|
670 tplarson 1.1
|
671 tplarson 1.5 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 tplarson 1.1
|
680 tplarson 1.10 if (verbflag > 1)
681 {
682 printf("starting l=%d\n", lmode);
683 }
684
|
685 tplarson 1.6 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 tplarson 1.1
|
696 tplarson 1.6 n_in=(tstop-tstart)/cadence;
|
697 tplarson 1.5 if (n_in != gapsize)
|
698 tplarson 1.3 {
|
699 tplarson 1.5 fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, n_in=%d\n", gapsize, n_in);
700 return 0;
|
701 tplarson 1.3 }
702
|
703 tplarson 1.5 tstartout=tstart+ifirst*tsample;
704 tstopout=tstartout+df1;
705 sprint_time(tstartstr, tstartout, "TAI", 0);
|
706 tplarson 1.1
|
707 tplarson 1.5 if (seginflag)
708 segin = drms_segment_lookup(inrec, segnamein);
|
709 tplarson 1.3 else
|
710 tplarson 1.5 segin = drms_segment_lookupnum(inrec, 0);
711 if (segin == NULL)
|
712 tplarson 1.3 {
|
713 tplarson 1.5 fprintf(stderr, "ERROR: problem looking up input segment: recnum = %lld\n", inrec->recnum);
|
714 tplarson 1.3 return 0;
715 }
|
716 tplarson 1.1
|
717 tplarson 1.5 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 tplarson 1.5 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 tplarson 1.3 // 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 tplarson 1.5 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 tplarson 1.5 ;
770 }
771 */
|
772 tplarson 1.10 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 tplarson 1.3 }
|
791 tplarson 1.1
|
792 tplarson 1.15 if (mavgout)
793 for (i=0;i<n_sampled_out/2;i++)
794 msum[i] = 0.0;
795
|
796 tplarson 1.1 /***** Main Loop over m: Detrend, gapfill and/or compute FFT/power spectrum *******/
|
797 tplarson 1.5 for (m=0; m<=lmode; m++)
798 {
|
799 tplarson 1.3
|
800 tplarson 1.10 if (verbflag > 2)
|
801 tplarson 1.5 {
|
802 tplarson 1.10 printf("starting m=%d\n", m);
|
803 tplarson 1.5 }
|
804 tplarson 1.3
|
805 tplarson 1.5 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 tplarson 1.1
815 /* Extracts data copy */
|
816 tplarson 1.5 extract_data(n_sampled_in, gaps, in_data, data);
|
817 tplarson 1.1
818 /* Detrend entire time series if required */
|
819 tplarson 1.5 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 tplarson 1.1 detrend_order, detrend_length, detrend_skip,
|
824 tplarson 1.3 num_sections, sections, detrend_first);
|
825 tplarson 1.5 }
|
826 tplarson 1.1
827 /* Fill gaps in entire time series if required */
|
828 tplarson 1.5 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 tplarson 1.1 fu_min_percent, max_ar_order, fu_iter, fu_fit_edges,
833 &ar_order, (_Complex float *)ar_coeff);
|
834 tplarson 1.5 }
|
835 tplarson 1.3
|
836 tplarson 1.5 drms_free_array(inarr);
|
837 tplarson 1.3
|
838 tplarson 1.5 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 tplarson 1.3
|
842 tplarson 1.5 if (tsout)
843 {
|
844 tplarson 1.3 // code to flip m on this output separately. in the present code this has already been accomplished by setting flipm.
|
845 tplarson 1.2 /*
|
846 tplarson 1.5 if (mfliparr[TSOUT] != mflipin)
|
847 tplarson 1.3 {
|
848 tplarson 1.5 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 tplarson 1.3 }
|
855 tplarson 1.5 else
856 tmpptr=data;
|
857 tplarson 1.2 */
|
858 tplarson 1.5 tsstart[1]=m;
859 tsend[1]=m;
860 outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
|
861 tplarson 1.9 outarr->bzero=segoutarr[TSOUT]->bzero;
862 outarr->bscale=segoutarr[TSOUT]->bscale;
|
863 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[TSOUT], outarr, tsstart, tsend, tstotal, 0);
|
864 tplarson 1.5 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 tplarson 1.3 }
|
871 tplarson 1.1
|
872 tplarson 1.5 if (fftout || fft1out || powout || mavgout)
873 fftwf_execute(plan);
|
874 tplarson 1.1
|
875 tplarson 1.5 if (fftout)
876 {
|
877 tplarson 1.3 // code to flip m on this output separately. in the present code this has already been accomplished by setting flipm.
878 /*
|
879 tplarson 1.5 if (mfliparr[FFTOUT] != mflipin)
|
880 tplarson 1.3 {
|
881 tplarson 1.5 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 tplarson 1.3 }
890 else
|
891 tplarson 1.5 tmpptr=data;
|
892 tplarson 1.3 */
|
893 tplarson 1.5 tsstart[1]=m;
894 tsend[1]=m;
895 outarr = drms_array_create(usetype, 2, tslength, tmpptr, &status);
|
896 tplarson 1.9 outarr->bzero=segoutarr[FFTOUT]->bzero;
897 outarr->bscale=segoutarr[FFTOUT]->bscale;
|
898 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[FFTOUT], outarr, tsstart, tsend, tstotal, 0);
|
899 tplarson 1.5 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 tplarson 1.3 }
906
|
907 tplarson 1.5 if (fft1out)
|
908 tplarson 1.3 {
909
|
910 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[FFT1OUT]->bzero;
914 outarr->bscale=segoutarr[FFT1OUT]->bscale;
|
915 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
|
916 tplarson 1.3 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 tplarson 1.1
|
923 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[FFT1OUT]->bzero;
936 outarr->bscale=segoutarr[FFT1OUT]->bscale;
|
937 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[FFT1OUT], outarr, fft1start, fft1end, fft1total, 0);
|
938 tplarson 1.5 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 tplarson 1.3 }
|
946 tplarson 1.1
|
947 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[POWOUT]->bzero;
957 outarr->bscale=segoutarr[POWOUT]->bscale;
|
958 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
|
959 tplarson 1.3 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 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[POWOUT]->bzero;
980 outarr->bscale=segoutarr[POWOUT]->bscale;
|
981 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[POWOUT], outarr, powstart, powend, powtotal, 0);
|
982 tplarson 1.5 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 tplarson 1.1 }
|
990 tplarson 1.3
|
991 tplarson 1.5 if (mavgout)
|
992 tplarson 1.1 {
|
993 tplarson 1.5 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 tplarson 1.1 }
|
1012 tplarson 1.3
|
1013 tplarson 1.5 if (arpowout)
|
1014 tplarson 1.1 {
|
1015 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[ARPOWOUT]->bzero;
1025 outarr->bscale=segoutarr[ARPOWOUT]->bscale;
|
1026 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
|
1027 tplarson 1.3 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 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[ARPOWOUT]->bzero;
1048 outarr->bscale=segoutarr[ARPOWOUT]->bscale;
|
1049 tplarson 1.10 status=drms_segment_writeslice_ext(segoutarr[ARPOWOUT], outarr, powstart, powend, powtotal, 0);
|
1050 tplarson 1.5 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 tplarson 1.3 }
|
1057 tplarson 1.1 }
1058
|
1059 tplarson 1.5 } //end of loop over m
|
1060 tplarson 1.1
|
1061 tplarson 1.5 if (mavgout)
|
1062 tplarson 1.3 {
|
1063 tplarson 1.5 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 tplarson 1.9 outarr->bzero=segoutarr[MAVGOUT]->bzero;
1068 outarr->bscale=segoutarr[MAVGOUT]->bscale;
|
1069 tplarson 1.5 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 tplarson 1.3 }
|
1077 tplarson 1.2
|
1078 tplarson 1.5 for (is=0;is<6;is++)
1079 {
1080 if (!flagarr[is])
1081 continue;
|
1082 tplarson 1.1
|
1083 tplarson 1.5 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 tplarson 1.13 if (tagflag)
1092 drms_setkey_string(outrec, "TAG", tag);
|
1093 tplarson 1.15 if (verflag)
1094 drms_setkey_string(outrec, "VERSION", version);
1095
|
1096 tplarson 1.5 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 tplarson 1.3
|
1102 tplarson 1.5 }
|
1103 tplarson 1.1
1104 }
1105
|
1106 tplarson 1.5 drms_close_records(inrecset, DRMS_FREE_RECORD);
1107 for (is=0;is<6;is++)
|
1108 tplarson 1.1 {
|
1109 tplarson 1.5 if (!flagarr[is])
1110 continue;
1111 drms_close_records(outrecsetarr[is], DRMS_INSERT_RECORD);
|
1112 tplarson 1.1 }
1113
|
1114 tplarson 1.5 if(ar_coeff != NULL)
|
1115 tplarson 1.1 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 tplarson 1.10 printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n",
|
1125 tplarson 1.1 wt-wt0, ct-ct0);
1126 }
1127
|
1128 tplarson 1.3 printf("module %s successful completion\n", cmdparams.argv[0]);
|
1129 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1
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 tplarson 1.3 // 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.3 // 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 tplarson 1.1 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 tplarson 1.1
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 tplarson 1.1 {
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 tplarson 1.1 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 tplarson 1.1 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 tplarson 1.1 /*
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 tplarson 1.1 *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 tplarson 1.2 return 2;
|
1422 tplarson 1.1 }
1423 }
1424 fclose(fh);
1425 return 0;
1426 }
1427 else
1428 return 1;
1429 }
1430
|
Return to
jtsfiddle.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps