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