JSOC/proj/globalhs/apps/jpkbgn.c - annotate

Return to jpkbgn.c CVS log

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

1 tplarson 1.1 /* this JSOC module is a port of an SOI module. 2 * ported by Tim Larson. 3 / 4 5 / 6 * peakbagging program - pkbgn24.c 7 */ 8 9 #include <stdio.h> 10 #include <stdlib.h> 11 #include <strings.h>
12 tplarson 1.5 #include <sys/types.h> 13 #include <sys/stat.h> 14 #include <unistd.h>
15 tplarson 1.1 16 #include "jsoc_main.h"
17 kehcheng 1.6 #include "fitsio.h"
18 tplarson 1.4 #include "drms_dsdsapi.h"
19 tplarson 1.1 20 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0])) 21 #define kNOTSPECIFIED "not specified" 22 23 char *module_name = "jpkbgn";
24 tplarson 1.13 char *cvsinfo_jpkbgn = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jpkbgn.c,v 1.12 2015/12/30 00:18:25 tplarson Exp $";
25 tplarson 1.1 26 ModuleArgs_t module_args[] = 27 { 28 {ARG_STRING, "in", "", "input data records"}, 29 {ARG_STRING, "out", "", "output data series"},
30 tplarson 1.9 {ARG_STRING, "ffttmp", "", ""},
31 tplarson 1.1 {ARG_STRING, "segin" , kNOTSPECIFIED, "name of input segment if not using segment 0"}, 32 {ARG_STRING, "segout", kNOTSPECIFIED, "name of output segment if not using segment 0"}, 33 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
34 tplarson 1.11 // {ARG_STRING, "srclink", "SRCDATA", "name of link to source data"}, 35 // {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
36 tplarson 1.8 {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},
37 tplarson 1.1 38 {ARG_STRING, "LOGFILE", "stdout", ""}, 39 {ARG_INT, "READFFT", "0", "", ""},
40 tplarson 1.2 {ARG_STRING, "GAPFILE", "", "record or file containing window function"},
41 tplarson 1.1 {ARG_INT, "IPAR", "1", "", ""}, 42 {ARG_STRING, "PARAMFILE", "nofile", ""}, 43 {ARG_STRING, "PAR1FILE", "", ""}, 44 {ARG_STRING, "RESFILE", "result", ""}, 45 {ARG_STRING, "CROSSFILE", "", ""}, 46 {ARG_INT, "NDELTAL", "6", "", ""}, 47 {ARG_INT, "NOISETYPE", "1", "", ""}, 48 {ARG_STRING, "CROSSNFILE", "/dev/null", ""}, 49 {ARG_INT, "NNOIB", "1", "", ""}, 50 51 {ARG_STRING, "NOIB", "40000 43000", ""}, 52 53 {ARG_INT, "NNOIM", "0", "", ""}, 54 {ARG_STRING, "NOIM", " ", ""}, 55 {ARG_FLOAT, "CSUBM", "0.5", "", ""}, 56 {ARG_INT, "IFOLLOW", "1", "", ""}, 57 {ARG_FLOAT, "ANOISE", "0.0", "", ""}, 58 {ARG_FLOAT, "AOTHER", "0.0", "", ""},
59 tplarson 1.2 // {ARG_INT, "LMODE", "", "", ""},
60 tplarson 1.1 {ARG_INT, "IODD", "1", "", ""}, 61 {ARG_INT, "NMIN", "0", "", ""}, 62 {ARG_INT, "NMAX", "50", "", ""}, 63 {ARG_FLOAT, "FMIN", "1000", "", ""}, 64 {ARG_FLOAT, "FMAX", "4000", "", ""}, 65 {ARG_INT, "IMFREQ", "0", "", ""}, 66 {ARG_INT, "ICASE", "3", "", ""}, 67 {ARG_INT, "NACOEFF", "5", "", ""}, 68 {ARG_INT, "NDT", "51840", "", ""},
69 tplarson 1.2 // {ARG_FLOAT, "TSAMPLE", "60", "", ""},
70 tplarson 1.1 {ARG_INT, "IDF", "1", "", ""}, 71 {ARG_INT, "NDF0", "20", "", ""}, 72 {ARG_INT, "NDF1", "50", "", ""}, 73 {ARG_FLOAT, "CDF", "5", "", ""}, 74 {ARG_INT, "NDM", "10", "", ""}, 75 {ARG_INT, "NDL", "6", "", ""}, 76 {ARG_FLOAT, "XSAFE", "10.0", "", ""}, 77 {ARG_FLOAT, "XSAFE1", "1.0", "", ""}, 78 {ARG_FLOAT, "DELMIN", "0.01", "", ""}, 79 {ARG_INT, "MAXCHI", "50", "", ""}, 80 {ARG_INT, "MAXGRAD", "30", "", ""}, 81 {ARG_INT, "MAXCOF", "4", "", ""}, 82 {ARG_INT, "IFIX", "0", "", ""}, 83 {ARG_INT, "IFIRST", "0", "", ""}, 84 {ARG_INT, "FLIPM", "0", "", ""}, 85 {ARG_INT, "NLOBES", "", "", ""}, 86 {ARG_INT, "NFOUR", "0", "", ""}, 87 {ARG_INT, "GONGFLAG", "0", "", ""}, 88 {ARG_INT, "PRINTFIT", "0", "", ""}, 89 {ARG_FLOAT, "CDETREND", "50.0", "", ""}, 90 {ARG_INT, "IFILL", "1", "", ""}, 91 tplarson 1.1 {ARG_INT, "FDIFF", "0", "", ""}, 92 {ARG_INT, "IASYM", "0", "", ""}, 93 {ARG_INT, "ICT", "0", "", ""}, 94 {ARG_INT, "IWOOD", "0", "", ""},
95 tplarson 1.11 {ARG_FLOAT, "WOODB1", "0.0", "", ""}, 96 {ARG_FLOAT, "WOODB2", "0.0", "", ""},
97 tplarson 1.1 {ARG_INT, "CTX", "1.0", "", ""}, 98 99 {ARG_END} 100 }; 101 102 #include "saveparm.c"
103 tplarson 1.9 #include "set_history.c"
104 tplarson 1.1
105 tplarson 1.9 void detrend(float data, float gaps, int length, float cdetrend); 106 void fill_gaps(float data, float gaps, float *newgaps, int length);
107 tplarson 1.10 //char *getpkbgapsversion(void);
108 tplarson 1.1
109 tplarson 1.11 void sub24_(char logfile,char fftfile, float rgaps, int ipar, char paramfile, char par1file, char resfile, 110 char crossfile, int ndeltal, int noisetype, char crossnfile, int nnoib, int noib, int nnoim, int noim, 111 float csubm, int ifol, float anoise, float aother, int lin, int ioddin, int iasymin, int ictin, float ctxin, 112 int iwoodin, float woodb1, float woodb2, int nmin, int nmax, float fmin, float fmax, 113 int imfin, int icasein, int nain, int idtin, int ndtin, float tsin, int idf, int ndf0, int ndf1, 114 float cdf, int ndmin, int ndl, 115 float xsin, float xs1in, float delx, int maxchi, int maxgrad, int maxcof, int recin, 116 int nlobesin, int nfourin, int gongin, int printflag,
117 tplarson 1.1 int, int, int, int, int, int, int); 118
119 tplarson 1.3 void cffti_(int , float ); 120 void cfftb_(int , float , float *); 121
122 tplarson 1.10 //void getpkbgnversion_(char *, int);
123 tplarson 1.9
124 tplarson 1.1 125 int DoIt(void) 126 { 127 int status = DRMS_SUCCESS; 128 int newstat = 0; 129 int noib[10], noim[10]; 130 131 int fstatus = 0; 132 fitsfile fitsptr; 133 long fdims[1]; 134 long fpixel[]={1}; 135 int anynul=0; 136 137 int zero=0; 138 int i, j; 139 float c,sum; 140 char ptr; 141 142 FILE fileptr; 143 144 int readfft, savefft; 145 tplarson 1.1 int in_nsets; 146 147 int fdiff,flipm,m,navail,nread,nx1; 148 float isign,*data; 149 TIME epoch, step, start, stop; 150 int gapsize, istart, istop;
151 tplarson 1.9 char ffttmp, fftfull, *infft;
152 tplarson 1.1 float cdetrend; 153 char *syscall; 154
155 tplarson 1.2 char gapfile, logfile, paramfile, par1file, resfile, crossfile, *crossnfile;
156 tplarson 1.1 int loglen, fftlen, paramlen, par1len, reslen, crosslen, crossnlen; 157 int ipar, ndeltal, noisetype, nnoib, nnoim; 158 char noibstr, noimstr; 159 float csubm, anoise, aother, tsample, cdf, xsafe, xsafe1, delmin, fmin, fmax; 160 int ifollow, lmode, iodd, nmin, nmax, imfreq, icase, nacoeff; 161 int ndt, idf, ndf0, ndf1, ndm, ndl, maxchi, maxgrad, maxcof; 162 163 int ifix, ifirst, ifill; 164 int nlobes,nfour,gongflag,printfit; 165 int iasym, ict, iwood; 166 float ctx;
167 tplarson 1.11 float woodb1, woodb2;
168 tplarson 1.1 169 DRMS_Segment_t segin = NULL; 170 DRMS_Segment_t segout = NULL; 171 DRMS_Array_t inarr = NULL; 172 DRMS_Array_t outarr = NULL; 173 DRMS_Record_t inrec = NULL; 174 DRMS_Record_t outrec = NULL; 175 DRMS_Type_t usetype = DRMS_TYPE_FLOAT; 176 DRMS_RecLifetime_t lifetime; 177 long long histrecnum = -1;
178 tplarson 1.2 DRMS_Segment_t gapseg = NULL; 179 DRMS_Array_t gaparr = NULL;
180 tplarson 1.9 DRMS_RecordSet_t ffttmprecset = NULL; 181 DRMS_Record_t ffttmprec = NULL; 182 DRMS_Segment_t ffttmpseg = NULL; 183 FILE ffttmpfile = NULL;
184 tplarson 1.2 char tstartstr[100]; 185 char fftfile[DRMS_MAXPATHLEN+1]; 186 char dirstr[DRMS_MAXPATHLEN+1];
187 tplarson 1.1 188 TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC / 189 double tstart, tstartsave, tstop, tstopsave, tstep; 190 191 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); 192 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); 193 194 /#ifdef __linux__/ 195 / Used to have to use this: 196 int regsz = 4; 197 After upgrade have to use this: / 198 #ifdef __ia64 199 int regsz = 4; 200 #else 201 int regsz = 1; 202 #endif 203 204 float rgaps, igaps; 205 float gaps, gaps1, gaps2, gaps3; 206 float help; 207 float datr, dati, data3; 208 tplarson 1.1 float datarr; 209
210 tplarson 1.9 ffttmp = (char *)cmdparams_save_str (&cmdparams, "ffttmp", &newstat);
211 tplarson 1.3 logfile = (char *)cmdparams_save_str (&cmdparams, "LOGFILE", &newstat); 212 readfft = cmdparams_save_int (&cmdparams, "READFFT", &newstat);
213 tplarson 1.9 // savefft = cmdparams_save_int (&cmdparams, "SAVEFFT", &newstat);
214 tplarson 1.3 gapfile = (char )cmdparams_save_str (&cmdparams, "GAPFILE", &newstat); 215 ipar = cmdparams_save_int (&cmdparams, "IPAR", &newstat); 216 paramfile = (char )cmdparams_save_str (&cmdparams, "PARAMFILE", &newstat); 217 par1file = (char )cmdparams_save_str (&cmdparams, "PAR1FILE", &newstat); 218 resfile = (char )cmdparams_save_str (&cmdparams, "RESFILE", &newstat); 219 crossfile = (char )cmdparams_save_str (&cmdparams, "CROSSFILE", &newstat); 220 ndeltal = cmdparams_save_int (&cmdparams, "NDELTAL", &newstat); 221 noisetype = cmdparams_save_int (&cmdparams, "NOISETYPE", &newstat); 222 crossnfile = (char )cmdparams_save_str (&cmdparams, "CROSSNFILE", &newstat); 223 nnoib = cmdparams_save_int (&cmdparams, "NNOIB", &newstat); 224 noibstr = (char )cmdparams_save_str (&cmdparams, "NOIB", &newstat); 225 nnoim = cmdparams_save_int (&cmdparams, "NNOIM", &newstat); 226 noimstr = (char )cmdparams_save_str (&cmdparams, "NOIM", &newstat); 227 csubm = cmdparams_save_float (&cmdparams, "CSUBM", &newstat); 228 ifollow = cmdparams_save_int (&cmdparams, "IFOLLOW", &newstat); 229 anoise = cmdparams_save_float (&cmdparams, "ANOISE", &newstat); 230 aother = cmdparams_save_float (&cmdparams, "AOTHER", &newstat);
231 tplarson 1.2 // lmode = cmdparams_save_int (&cmdparams, "LMODE", &newstat);
232 tplarson 1.3 iodd = cmdparams_save_int (&cmdparams, "IODD", &newstat); 233 nmin = cmdparams_save_int (&cmdparams, "NMIN", &newstat); 234 nmax = cmdparams_save_int (&cmdparams, "NMAX", &newstat); 235 fmin = cmdparams_save_float (&cmdparams, "FMIN", &newstat); 236 fmax = cmdparams_save_float (&cmdparams, "FMAX", &newstat); 237 imfreq = cmdparams_save_int (&cmdparams, "IMFREQ", &newstat); 238 icase = cmdparams_save_int (&cmdparams, "ICASE", &newstat); 239 nacoeff = cmdparams_save_int (&cmdparams, "NACOEFF", &newstat); 240 ndt = cmdparams_save_int (&cmdparams, "NDT", &newstat);
241 tplarson 1.2 // tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat);
242 tplarson 1.3 idf = cmdparams_save_int (&cmdparams, "IDF", &newstat); 243 ndf0 = cmdparams_save_int (&cmdparams, "NDF0", &newstat); 244 ndf1 = cmdparams_save_int (&cmdparams, "NDF1", &newstat); 245 cdf = cmdparams_save_float (&cmdparams, "CDF", &newstat); 246 ndm = cmdparams_save_int (&cmdparams, "NDM", &newstat); 247 ndl = cmdparams_save_int (&cmdparams, "NDL", &newstat); 248 xsafe = cmdparams_save_float (&cmdparams, "XSAFE", &newstat); 249 xsafe1 = cmdparams_save_float (&cmdparams, "XSAFE1", &newstat); 250 delmin = cmdparams_save_float (&cmdparams, "DELMIN", &newstat); 251 maxchi = cmdparams_save_int (&cmdparams, "MAXCHI", &newstat); 252 maxgrad = cmdparams_save_int (&cmdparams, "MAXGRAD", &newstat); 253 maxcof = cmdparams_save_int (&cmdparams, "MAXCOF", &newstat); 254 ifix = cmdparams_save_int (&cmdparams, "IFIX", &newstat); 255 ifirst = cmdparams_save_int (&cmdparams, "IFIRST", &newstat); 256 flipm = cmdparams_save_int (&cmdparams, "FLIPM", &newstat); 257 nlobes = cmdparams_save_int (&cmdparams, "NLOBES", &newstat); 258 nfour = cmdparams_save_int (&cmdparams, "NFOUR", &newstat); 259 gongflag = cmdparams_save_int (&cmdparams, "GONGFLAG", &newstat); 260 printfit = cmdparams_save_int (&cmdparams, "PRINTFIT", &newstat); 261 cdetrend = cmdparams_save_float (&cmdparams, "CDETREND", &newstat); 262 ifill = cmdparams_save_int (&cmdparams, "IFILL", &newstat); 263 tplarson 1.3 fdiff = cmdparams_save_int (&cmdparams, "FDIFF", &newstat); 264 iasym = cmdparams_save_int (&cmdparams, "IASYM", &newstat); 265 ict = cmdparams_save_int (&cmdparams, "ICT", &newstat); 266 iwood = cmdparams_save_int (&cmdparams, "IWOOD", &newstat);
267 tplarson 1.11 woodb1 = cmdparams_save_float (&cmdparams, "WOODB1", &newstat); 268 woodb2 = cmdparams_save_float (&cmdparams, "WOODB2", &newstat);
269 tplarson 1.3 ctx = cmdparams_save_int (&cmdparams, "CTX", &newstat);
270 tplarson 1.13 //printf("ndl = %d\n",ndl);
271 tplarson 1.3 char inrecquery = (char )cmdparams_save_str(&cmdparams, "in", &newstat); 272 char outseries = (char )cmdparams_save_str(&cmdparams, "out", &newstat); 273 char segnamein = (char )cmdparams_save_str(&cmdparams, "segin", &newstat); 274 char segnameout = (char )cmdparams_save_str(&cmdparams, "segout", &newstat);
275 tplarson 1.1 int seginflag = strcmp(kNOTSPECIFIED, segnamein); 276 int segoutflag = strcmp(kNOTSPECIFIED, segnameout);
277 tplarson 1.3 char histlinkname = (char )cmdparams_save_str(&cmdparams, "histlink", &newstat);
278 tplarson 1.9 // char srclinkname = (char )cmdparams_save_str(&cmdparams, "srclink", &newstat); 279 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat); 280 /* 281 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
282 tplarson 1.1 if (permflag) 283 lifetime = DRMS_PERMANENT; 284 else 285 lifetime = DRMS_TRANSIENT;
286 tplarson 1.9 */
287 tplarson 1.1
288 tplarson 1.9 if (newstat)
289 tplarson 1.1 {
290 tplarson 1.9 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat); 291 cpsave_decode_error(newstat); 292 return 1; 293 } 294 else if (savestrlen != strlen(savestr))
295 tplarson 1.1 {
296 tplarson 1.9 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr)); 297 return 1;
298 tplarson 1.1 }
299 tplarson 1.9 300 301 int histflag = strncasecmp("none", outseries, 4); 302 if (histflag)
303 tplarson 1.1 {
304 tplarson 1.9 long long histrecnum; 305 DRMS_Record_t tempoutrec = drms_create_record(drms_env, 306 outseries, 307 DRMS_TRANSIENT, 308 &status); 309 310 if (status != DRMS_SUCCESS) 311 { 312 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status); 313 return 1; 314 } 315 316 DRMS_Link_t histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname); 317 if (histlink != NULL) 318 {
319 tplarson 1.10 histrecnum=set_history(histlink);
320 tplarson 1.9 if (histrecnum < 0) 321 { 322 fprintf(stderr, "ERROR: problem creating record in history dataseries for output dataseries %s\n", outseries); 323 drms_close_record(tempoutrec, DRMS_FREE_RECORD); 324 return 1; 325 } 326 printf("histrecnum=%ld\n",histrecnum); 327 } 328 else 329 { 330 fprintf(stderr,"ERROR: could not find history link in output dataseries %s\n", outseries); 331 return 1; 332 } 333 334 drms_close_record(tempoutrec, DRMS_FREE_RECORD); 335 printf("module %s successful completion\n", cmdparams.argv[0]); 336 return 0; 337
338 tplarson 1.1 }
339 tplarson 1.2
340 tplarson 1.1 341 /* Set PARAMFILE to PAR1FILE if not set */ 342 343 if (strcmp(paramfile,"nofile") == 0) 344 paramfile=par1file; 345
346 tplarson 1.5 char crosspath; 347 char dir, sub, hold; 348 struct stat crosstat; 349 if (fstatus=stat(crossfile, &crosstat)) 350 { 351 DRMS_RecordSet_t crossrecset = drms_open_records(drms_env, crossfile, &status); 352 if (status != DRMS_SUCCESS \|\| crossrecset == NULL) 353 { 354 fprintf(stderr,"ERROR: problem finding leakage matrix: file status = %d, DRMS status = %d\n",fstatus,status); 355 return 1; 356 } 357 if (crossrecset->n == 0) 358 { 359 fprintf(stderr,"ERROR: crossfile recordset contains no records\n"); 360 return 1; 361 } 362 if (crossrecset->n > 1) 363 { 364 fprintf(stderr,"ERROR: crossfile recordset with more than 1 record not supported.\n"); 365 return 1; 366 } 367 tplarson 1.5 DRMS_Segment_t crosseg = drms_segment_lookupnum(crossrecset->records[0], 0); 368 crosspath = (char *)malloc(DRMS_MAXPATHLEN+1); 369 if (crosseg != NULL && crosspath != NULL) 370 status=drms_record_directory(crossrecset->records[0],crosspath,0); 371 if (crosseg == NULL \|\| status != DRMS_SUCCESS \|\| crosspath[0] == '\0') 372 { 373 fprintf(stderr, "ERROR: problem finding leakage matrix segment: status = %d\n", status); 374 return 1; 375 } 376 dir=drms_getkey_string(crossrecset->records[0], "dirname", &status); 377 sub=strtok(dir,"/"); 378 while ((hold=strtok(NULL,"/")) != NULL) 379 sub=hold; 380 strcat(crosspath,"/"); 381 strcat(crosspath,sub); 382 strcat(crosspath,"/");
383 tplarson 1.8 if (verbflag) 384 printf("crosspath = %s\n", crosspath);
385 tplarson 1.5 } 386 else 387 { 388 crosspath=strdup(crossfile); 389 }
390 tplarson 1.2 391 /* Read window function / 392 if (fits_open_file(&fitsptr, gapfile, READONLY, &fstatus)) 393 { 394 DRMS_RecordSet_t gaprecset = drms_open_records(drms_env, gapfile, &status); 395 if (status != DRMS_SUCCESS \|\| gaprecset == NULL) 396 { 397 fprintf(stderr,"ERROR: problem reading gaps: file status = %d, DRMS status = %d\n",fstatus,status);
398 tplarson 1.3 return 1;
399 tplarson 1.2 } 400 if (gaprecset->n == 0) 401 { 402 fprintf(stderr,"ERROR: gapfile recordset contains no records\n");
403 tplarson 1.3 return 1;
404 tplarson 1.2 } 405 if (gaprecset->n > 1) 406 { 407 fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n");
408 tplarson 1.3 return 1;
409 tplarson 1.2 } 410 gapseg = drms_segment_lookupnum(gaprecset->records[0], 0); 411 if (gapseg != NULL) 412 gaparr = drms_segment_read(gapseg, DRMS_TYPE_FLOAT, &status); 413 if (status != DRMS_SUCCESS \|\| gaparr == NULL \|\| gapseg == NULL) 414 { 415 fprintf(stderr, "ERROR: problem reading gaps segment: status = %d\n", status);
416 tplarson 1.3 return 1;
417 tplarson 1.2 } 418 rgaps=(float *)(gaparr->data); 419 gapsize=gaparr->axis[0]; 420 }
421 tplarson 1.1 else 422 {
423 tplarson 1.2 fits_get_img_size(fitsptr, 1, fdims, &fstatus); 424 gapsize=fdims[0]; 425 rgaps=(float )(malloc(gapsizesizeof(float))); 426 fits_read_pix(fitsptr, TFLOAT, fpixel, gapsize, NULL, rgaps, anynul, &fstatus); 427 fits_close_file(fitsptr, &fstatus); 428 if (fstatus) 429 { 430 fits_report_error(stderr, fstatus);
431 tplarson 1.3 return 1;
432 tplarson 1.2 }
433 tplarson 1.1 } 434
435 tplarson 1.8 if (verbflag) 436 printf("gapfile read, gapsize = %d\n",gapsize);
437 tplarson 1.1 438 igaps=(float )(malloc(gapsizesizeof(float))); 439 gaps=(float )(malloc(gapsizesizeof(float))); 440 gaps1=(float )(malloc(gapsizesizeof(float))); 441 gaps2=(float )(malloc(gapsizesizeof(float))); 442 datr=(float )(malloc(gapsizesizeof(float))); 443 dati=(float )(malloc(gapsizesizeof(float))); 444 /* 445 help=(float )(malloc((4gapsize+30)sizeof(float))); 446 Changed 20061022 / 447 help=(float )(malloc((4ndt+30)sizeof(float))); 448 gaps3=(float )(malloc(ndtsizeof(float))); 449 data3=(float )(malloc(2ndtsizeof(float))); 450 451 isign = 1.0;
452 tplarson 1.3 if (flipm != 0)
453 tplarson 1.1 isign=-1.0; 454 sum=0.0; 455 for (i=0; i<gapsize; i++) 456 { 457 sum=sum+rgaps[i]; 458 igaps[i] = isign*rgaps[i]; 459 if (rgaps[i] == 0.0) 460 gaps[i]=0.0; 461 else 462 gaps[i]=1.0; 463 }
464 tplarson 1.8 if (verbflag) 465 printf("%f\n",sum);
466 tplarson 1.1 467 cffti_(&ndt, help); 468 469
470 tplarson 1.3 DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status); 471 if (status != DRMS_SUCCESS \|\| inrecset == NULL)
472 tplarson 1.2 { 473 fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status);
474 tplarson 1.3 return 1;
475 tplarson 1.2 } 476
477 tplarson 1.3 if (inrecset->n == 0)
478 tplarson 1.2 {
479 tplarson 1.8 fprintf(stderr,"ERROR: input recordset contains no records\n");
480 tplarson 1.3 return 1;
481 tplarson 1.2 } 482
483 tplarson 1.3 if (inrecset->n > 1)
484 tplarson 1.2 {
485 tplarson 1.8 fprintf(stderr,"ERROR: nrecs > 1 not yet supported.\n");
486 tplarson 1.3 return 1;
487 tplarson 1.2 } 488
489 tplarson 1.1
490 tplarson 1.3 inrec=inrecset->records[0];
491 tplarson 1.2 int itest; 492 char *inchecklist[] = {"T_START", "T_STOP", "LMIN", "LMAX", "T_STEP"};
493 tplarson 1.1
494 tplarson 1.2 for (itest=0; itest < ARRLENGTH(inchecklist); itest++) 495 { 496 DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]);
497 tplarson 1.3 if (inkeytest == NULL)
498 tplarson 1.2 { 499 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
500 tplarson 1.3 drms_close_records(inrecset, DRMS_FREE_RECORD); 501 return 1;
502 tplarson 1.2 } 503 }
504 tplarson 1.1 505 tstart=tstartsave=drms_getkey_time(inrec, "T_START", NULL); 506 tstop =tstopsave =drms_getkey_time(inrec, "T_STOP", NULL); 507 tstep=drms_getkey_float(inrec, "T_STEP", NULL); 508 509 navail=(tstop-tstart)/tstep; 510 if (navail != gapsize) 511 {
512 tplarson 1.3 fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, navail=%d\n", gapsize, navail);
513 tplarson 1.1 return 0; 514 } 515 516 nread=navail;
517 tplarson 1.2 tsample=tstep;
518 tplarson 1.1 519 if ((ifirst+ndt) < navail) 520 nread=ifirst+ndt; 521
522 tplarson 1.8 if (verbflag) 523 printf("%d %d\n",navail,nread);
524 tplarson 1.1
525 tplarson 1.2 int lmin=drms_getkey_int(inrec, "LMIN", NULL); 526 int lmax=drms_getkey_int(inrec, "LMAX", NULL); 527 if (lmin != lmax) 528 {
529 tplarson 1.8 fprintf(stderr,"ERROR: lmin != lmax not yet supported.\n");
530 tplarson 1.2 return 0; 531 } 532 lmode=lmin; 533
534 tplarson 1.7 char *tag = drms_getkey_string(inrec, "TAG", &status); 535 if (status != DRMS_SUCCESS) 536 tag=strdup("none");
537 tplarson 1.2 if (readfft == 0) 538 {
539 tplarson 1.9 ffttmprec = drms_create_record(drms_env, ffttmp, DRMS_PERMANENT, &status);
540 tplarson 1.2 541 if (status != DRMS_SUCCESS) 542 {
543 tplarson 1.9 fprintf(stderr,"ERROR: couldn't create a record in ffttmp dataseries %s, status = %d\n", ffttmp, status);
544 tplarson 1.2 return 0; 545 } 546
547 tplarson 1.9 drms_copykeys(ffttmprec, inrec, 0, kDRMS_KeyClass_Explicit); 548 drms_setkey_int(ffttmprec, "LMODE", lmode); 549 drms_setkey_int(ffttmprec, "NDT", ndt);
550 tplarson 1.8 tnow = (double)time(NULL); 551 tnow += UNIX_epoch;
552 tplarson 1.9 drms_setkey_time(ffttmprec, "DATE", tnow);
553 tplarson 1.8
554 tplarson 1.9 ffttmpseg = drms_segment_lookupnum(ffttmprec, 0);
555 tplarson 1.2 556 sprintf(fftfile, "fft%d", lmode);
557 tplarson 1.9 ffttmpfile = drms_segment_fopen(ffttmpseg, fftfile, 0, &status);
558 tplarson 1.2 if (status != DRMS_SUCCESS) 559 {
560 tplarson 1.9 fprintf(stderr,"ERROR: couldn't open a file to write ffttmp record, status = %d\n", status);
561 tplarson 1.2 return 0; 562 } 563
564 tplarson 1.9 drms_segment_filename(ffttmpseg,fftfile);
565 tplarson 1.2 566 } 567 else 568 {
569 tplarson 1.9 char *ffttmpquery = malloc(100);
570 tplarson 1.2 sprint_time(tstartstr, tstart, "TAI", 0);
571 tplarson 1.9 sprintf(ffttmpquery, "%s[%s][%d][%d][%s]", ffttmp,tstartstr,lmode,ndt,tag); 572 //printf("ffttmpquery = %s\n", ffttmpquery); 573 ffttmprecset = drms_open_records(drms_env, ffttmpquery, &status);
574 tplarson 1.12 if (status != DRMS_SUCCESS \|\| ffttmprecset == NULL \|\| ffttmprecset->n == 0)
575 tplarson 1.2 {
576 tplarson 1.9 fprintf(stderr,"ERROR: couldn't open ffttmp record %s, status = %d\n", ffttmpquery, status);
577 tplarson 1.2 return 0; 578 } 579
580 tplarson 1.9 ffttmprec = ffttmprecset->records[0]; 581 ffttmpseg = drms_segment_lookupnum(ffttmprec, 0);
582 tplarson 1.2
583 tplarson 1.9 drms_segment_filename(ffttmpseg,fftfile);
584 tplarson 1.2 585 } 586
587 tplarson 1.8 if (verbflag) 588 printf("READFFT=%d,fftfull=%s\n",readfft,fftfile);
589 tplarson 1.2
590 tplarson 1.8 //printf("cdetrend=%f,ifill=%d\n",cdetrend,ifill);
591 tplarson 1.2
592 tplarson 1.1 /* Loop over m and do FFT's */ 593 if (readfft == 0) 594 {
595 tplarson 1.3 if (seginflag) 596 segin = drms_segment_lookup(inrec, segnamein); 597 else 598 segin = drms_segment_lookupnum(inrec, 0); 599 if (segin) 600 inarr = drms_segment_read(segin, usetype, &status); 601 if (status != DRMS_SUCCESS \|\| inarr == NULL \|\| segin == NULL) 602 {
603 tplarson 1.8 fprintf(stderr, "problem reading input segment: status = %d\n", status);
604 tplarson 1.3 return 0; 605 } 606 607 datarr=(float *)(inarr->data); 608
609 tplarson 1.1 for (m=0; m<= lmode; m++) 610 { 611
612 tplarson 1.2 // printf("%d %d\n",ndt,m);
613 tplarson 1.1 614 data=datarr + m2gapsize; 615 616 for (j=0; j<navail; j++) 617 { 618 if (isnan(data[2j])) 619 datr[j]=0.0; 620 else 621 datr[j]=rgaps[j]data[2j]; 622 if (isnan(data[2j+1])) 623 dati[j]=0.0; 624 else 625 dati[j]=igaps[j]data[2j+1]; 626 } 627 628 if (cdetrend != 0.0) 629 { 630 detrend(datr,gaps,navail,cdetrend); 631 detrend(dati,gaps,navail,cdetrend); 632 } 633 634 tplarson 1.1 if (ifill == 0) 635 { 636 for (i=0; i<navail; i++) 637 { 638 gaps1[i]=gaps[i]; 639 } 640 } 641 else 642 { 643 fill_gaps(datr,gaps,gaps1,navail); 644 if (m !=0 ) 645 { 646 fill_gaps(dati,gaps,gaps1,navail); 647 } 648 } 649 650 if (fdiff!=0) 651 { 652 for (j=0; j<(navail-1); j++) 653 { 654 if ((gaps1[j]==1) & (gaps1[j+1]==1)) 655 tplarson 1.1 { 656 gaps2[j]=1.0; 657 datr[j]=datr[j+1]-datr[j]; 658 dati[j]=dati[j+1]-dati[j]; 659 } 660 else 661 { 662 gaps2[j]=0.0; 663 datr[j]=0.0; 664 dati[j]=0.0; 665 } 666 } 667 gaps2[navail-1]=0.0; 668 datr[navail-1]=0.0; 669 dati[navail-1]=0.0; 670 } 671 else 672 { 673 for (j=0; j<navail; j++) 674 { 675 gaps2[j]=gaps1[j]; 676 tplarson 1.1 } 677 } 678 679 /* Make complex array. 680 Now use data3 instead of reusing data in order to allow ndt>navail. / 681 nx1=ndt; 682 if (nread < (ifirst+ndt)) 683 nx1=nread-ifirst; 684 685 for (j=0; j<nx1; j++) 686 { 687 gaps3[j]=gaps2[j+ifirst]; 688 data3[2j]=datr[j+ifirst]; 689 data3[2j+1]=dati[j+ifirst]; 690 } 691 for (j = nx1; j < ndt; j++) 692 { 693 gaps3[j] = 0; 694 data3[2j] = 0; 695 data3[2j+1] = 0; 696 } 697 tplarson 1.1 698 if ((ifix == 1) && (m == 0)) 699 { 700 c = sqrt(2.0); 701 for (j = 0; j < ndt; j++) 702 { 703 data3[2j] = c; 704 data3[2j+1] *= c; 705 } 706 } 707 708 cfftb_ (&ndt, data3, help); 709
710 tplarson 1.9 fwrite ( (char)data3, 8ndt, 1, ffttmpfile);
711 tplarson 1.1 712 713 } /* End of m loop */ 714
715 tplarson 1.3 /* Pad the file for buffered reading in Fortran when NDT is small */
716 tplarson 1.9 fwrite ( (char)data3, 8ndt, 1, ffttmpfile); 717 drms_segment_fclose(ffttmpfile);
718 tplarson 1.8 if (verbflag) 719 printf("fft file written\n");
720 tplarson 1.1 } 721 else 722 { 723 /* Fill dummy time-series to make gaps */ 724 for (j=0; j<navail; j++) 725 { 726 datr[j]=j; 727 } 728 if (cdetrend != 0.0) 729 { 730 detrend(datr,gaps,navail,cdetrend); 731 } 732 if (ifill == 0) 733 { 734 for (i=0; i<navail; i++) 735 { 736 gaps1[i]=gaps[i]; 737 } 738 } 739 else 740 { 741 tplarson 1.1 fill_gaps(datr,gaps,gaps1,navail); 742 } 743 if (fdiff!=0) 744 { 745 for (j=0; j<(navail-1); j++) 746 { 747 if ((gaps1[j]==1) & (gaps1[j+1]==1)) 748 { 749 gaps2[j]=1.0; 750 } 751 else 752 { 753 gaps2[j]=0.0; 754 } 755 } 756 gaps2[navail-1]=0.0; 757 } 758 else 759 { 760 for (j=0; j<navail; j++) 761 { 762 tplarson 1.1 gaps2[j]=gaps1[j]; 763 } 764 } 765 nx1=ndt; 766 if (nread < (ifirst+ndt)) 767 nx1=nread-ifirst; 768 for (j=0; j<nx1; j++) 769 { 770 gaps3[j]=gaps2[j+ifirst]; 771 } 772 for (j = nx1; j < ndt; j++) 773 { 774 gaps3[j] = 0; 775 } 776 } 777 778 drms_free_array(inarr);
779 tplarson 1.3 if (inrecset)
780 tplarson 1.1 {
781 tplarson 1.3 drms_close_records(inrecset, DRMS_FREE_RECORD);
782 tplarson 1.1 }
783 tplarson 1.8 if (verbflag) 784 printf("input recordset closed\n");
785 tplarson 1.1 786 i = -1; 787 for ( ptr = noibstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) 788 noib[++i] = atoi(ptr); 789 790 i = -1; 791 for ( ptr = noimstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) 792 noim[++i] = atoi(ptr); 793 794 loglen = strlen(logfile);
795 tplarson 1.2 fftlen = strlen(fftfile);
796 tplarson 1.1 paramlen = strlen(paramfile); 797 par1len = strlen(par1file); 798 reslen = strlen(resfile);
799 tplarson 1.5 crosslen = strlen(crosspath);
800 tplarson 1.1 crossnlen = strlen(crossnfile);
801 tplarson 1.9 printf("crosspath=%s, crosslen=%d\n",crosspath,crosslen);
802 tplarson 1.13 printf("ndl = %d\n",ndl);
803 tplarson 1.1 /* Call Fortran program */
804 tplarson 1.2 sub24_(logfile,fftfile, gaps3, &ipar,
805 tplarson 1.1 paramfile, par1file, resfile,
806 tplarson 1.5 crosspath, &ndeltal, &noisetype, crossnfile, &nnoib, noib, &nnoim, noim,
807 tplarson 1.1 &csubm, &ifollow, &anoise, &aother, &lmode, &iodd, &iasym, &ict, &ctx,
808 tplarson 1.11 &iwood, &woodb1, &woodb2, &nmin, &nmax, &fmin, &fmax,
809 tplarson 1.1 &imfreq, &icase, &nacoeff, &zero, &ndt, &tsample, &idf, &ndf0, &ndf1, 810 &cdf, &ndm, &ndl, 811 &xsafe, &xsafe1, &delmin, &maxchi, &maxgrad, &maxcof, &regsz, 812 &nlobes, &nfour, 813 &gongflag, &printfit, 814 loglen, fftlen, paramlen, par1len, reslen, crosslen, crossnlen); 815
816 tplarson 1.3 if (readfft == 0)
817 tplarson 1.9 drms_close_record(ffttmprec,DRMS_INSERT_RECORD);
818 tplarson 1.3 else
819 tplarson 1.9 drms_close_records(ffttmprecset,DRMS_FREE_RECORD);
820 tplarson 1.1
821 tplarson 1.9 printf("module %s successful completion\n", cmdparams.argv[0]);
822 tplarson 1.1 return 0; 823 } 824 825 826 void detrend( 827 float data, 828 float gaps, 829 int length, 830 float cdetrend 831 ) 832 { 833 #define chunksz 1440 834 #define maxpols 50 835 836 extern float sdot_(int , float , int , float , int ); 837 extern float saxpy_(int , float , float , int , float , int ); 838 extern float sscal_(int , float , float , int ); 839 840 int i,j,k,l,n_good,ndegree; 841 int goodlist[chunksz]; 842 float pols[maxpols][chunksz]; 843 tplarson 1.1 float a,cg,rms2,sum,x[chunksz]; 844 float help; 845 int one=1; 846 847 for (i=0; i<length ; i=i+chunksz) { 848 help=data+i; 849 n_good=0; 850 cg=0.0; 851 for (j=i; j<i+chunksz; j++) { 852 if (gaps[j] != 0) { 853 goodlist[n_good]=j-i; 854 cg=cg+j-i; 855 n_good++; 856 } 857 } 858 if (n_good != 0) { 859 cg=cg/n_good; 860 rms2=0.0; 861 for (j=0; j<n_good; j++) {rms2=rms2+(goodlist[j]-cg)(goodlist[j]-cg);} 862 ndegree = floor((goodlist[n_good-1]-goodlist[0])/cdetrend)+2; 863 for (k=0; k<n_good; k++) { 864 tplarson 1.1 x[k]=(goodlist[k]-cg)/sqrt(rms2); 865 } 866 for (k=0; k<n_good; k++) {pols[0][k]=1.0/sqrt(n_good);} 867 for (k=0; k<n_good; k++) {pols[1][k]=x[k];} 868 for (j=2; j<ndegree; j++) { 869 for (k=0; k<n_good; k++) { 870 pols[j][k]=(2.0j)/jx[k]pols[j-1][k]-(j-1.0)/jpols[j-2][k]; 871 } 872 for (l=0; l<j; l++) { 873 a=-sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one); 874 saxpy_(&n_good,&a,&pols[l][0],&one,&pols[j][0],&one); 875 } 876 a=1.0/sqrt(sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one)); 877 sscal_(&n_good,&a,&pols[j][0],&one); 878 } 879 for (j=0; j<ndegree; j++) { 880 sum=0.0; 881 for (k=0; k<n_good; k++) { 882 sum=sum+pols[j][k]help[goodlist[k]]; 883 } 884 for (k=0; k<n_good; k++) { 885 tplarson 1.1 help[goodlist[k]] -= sumpols[j][k]; 886 } 887 } 888 } / if n_good / 889 890 } / for i */ 891 892 } 893
894 tplarson 1.9 void memcof(float , int, int, float , float ); 895 void fixrts(float , int); 896 void predic(float , int, float , int, float *, int);
897 tplarson 1.1 898 void fill_gaps( 899 float data, 900 float gaps, 901 float *newgaps, 902 int length 903 ) 904 {
905 tplarson 1.10 #include "nr.h" 906 #include "nrutil.h"
907 tplarson 1.1 #define chunksz 1440 908 #define maxpoles 10 909 #define maxgap 5 910 911 int i,j,k; 912 int n_good, gap_start, gap_end, gap_size, type; 913 float help,good_points[chunksz]; 914 int goodlist[chunksz]; 915 int npoles; 916 float pm,cof[maxpoles]; 917 float input[maxpoles]; 918 float f_predic[maxgap],b_predic[maxgap]; 919 920 npoles=6; 921 for (i=0; i<length ; i=i+chunksz) { 922 help=data+i; 923 n_good=0; 924 for (j=i; j<i+chunksz; j++) { 925 if (gaps[j] != 0) { 926 goodlist[n_good]=j-i; 927 good_points[n_good]=data[j]; 928 tplarson 1.1 n_good++; 929 newgaps[j]=gaps[j]; 930 } 931 } 932 if (n_good != 0) { 933 memcof(good_points-1,n_good,npoles,&pm,cof-1); 934 fixrts(cof-1,npoles); 935 for (j=1; j<n_good; j++) { 936 if ((goodlist[j]-goodlist[j-1]) > 1) { 937 gap_start=goodlist[j-1]+1; 938 gap_end=goodlist[j]-1; 939 gap_size=gap_end-gap_start+1; 940 if (gap_size <= maxgap) { 941 type=0; 942 / Forward prediction / 943 / Check that there are npoles good points before gap / 944 if ((j >= npoles) && (goodlist[j-npoles] == goodlist[j-1]-npoles+1)) { 945 for (k=0; k<npoles; k++) input[k]=help[gap_start-npoles+k]; 946 predic(input-1,npoles,cof-1,npoles,f_predic-1,gap_size); 947 type +=1; 948 } 949 tplarson 1.1 / Backwards prediction / 950 / Check that there are npoles good points after gap / 951 if ((j+npoles <= n_good) && (goodlist[j+npoles-1] == goodlist[j]+npoles-1)) { 952 for (k=0; k<npoles; k++) input[k]=help[gap_end+npoles-k]; 953 / for (k=0; k<npoles; k++) printf("%d %d %f\n",gap_end,k,help[gap_end+npoles-k]); / 954 predic(input-1,npoles,cof-1,npoles,b_predic-1,gap_size); 955 type +=2; 956 } 957 if (type == 3) { 958 / Average predictions / 959 for (k=0; k<gap_size; k++) 960 help[gap_start+k] = (f_predic[k]+b_predic[gap_size-k-1])/2.; 961 } 962 if (type == 1) { 963 / Use forward prediction / 964 for (k=0; k<gap_size; k++) 965 help[gap_start+k] = f_predic[k]; 966 } 967 if (type == 2) { 968 / Use backwards prediction / 969 for (k=0; k<gap_size; k++) 970 tplarson 1.1 help[gap_start+k] = b_predic[gap_size-k-1]; 971 } 972 if (type !=0 ) { 973 for (k=0; k<gap_size; k++) 974 newgaps[gap_start+k+i]=1; 975 } 976 } 977 } 978 } / for j / 979 } / if n_good / 980 } / for i */ 981 982 }

Karen Tian