| 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 |
#include <sys/types.h> |
| 13 |
#include <sys/stat.h> |
| 14 |
#include <unistd.h> |
| 15 |
|
| 16 |
#include "jsoc_main.h" |
| 17 |
#include "fitsio.h" |
| 18 |
#include "drms_dsdsapi.h" |
| 19 |
|
| 20 |
#define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0])) |
| 21 |
#define kNOTSPECIFIED "not specified" |
| 22 |
|
| 23 |
char *module_name = "jpkbgn"; |
| 24 |
char *cvsinfo_jpkbgn = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/jpkbgn.c,v 1.14 2019/10/15 17:29:21 schou Exp $"; |
| 25 |
|
| 26 |
ModuleArgs_t module_args[] = |
| 27 |
{ |
| 28 |
{ARG_STRING, "in", "", "input data records"}, |
| 29 |
{ARG_STRING, "out", "", "output data series"}, |
| 30 |
{ARG_STRING, "ffttmp", "", ""}, |
| 31 |
{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 |
// {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 |
{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 |
|
| 38 |
{ARG_STRING, "LOGFILE", "stdout", ""}, |
| 39 |
{ARG_INT, "READFFT", "0", "", ""}, |
| 40 |
{ARG_STRING, "GAPFILE", "", "record or file containing window function"}, |
| 41 |
{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 |
// {ARG_INT, "LMODE", "", "", ""}, |
| 60 |
{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 |
// {ARG_FLOAT, "TSAMPLE", "60", "", ""}, |
| 70 |
{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 |
{ARG_INT, "FDIFF", "0", "", ""}, |
| 92 |
{ARG_INT, "IASYM", "0", "", ""}, |
| 93 |
{ARG_INT, "ICT", "0", "", ""}, |
| 94 |
{ARG_INT, "IWOOD", "0", "", ""}, |
| 95 |
{ARG_FLOAT, "WOODB1", "0.0", "", ""}, |
| 96 |
{ARG_FLOAT, "WOODB2", "0.0", "", ""}, |
| 97 |
{ARG_INT, "CTX", "1.0", "", ""}, |
| 98 |
|
| 99 |
{ARG_END} |
| 100 |
}; |
| 101 |
|
| 102 |
#include "saveparm.c" |
| 103 |
#include "set_history.c" |
| 104 |
|
| 105 |
void detrend(float *data, float *gaps, int length, float cdetrend); |
| 106 |
void fill_gaps(float *data, float *gaps, float *newgaps, int length); |
| 107 |
//char *getpkbgapsversion(void); |
| 108 |
|
| 109 |
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 |
int, int, int, int, int, int, int); |
| 118 |
|
| 119 |
void cffti_(int *, float *); |
| 120 |
void cfftb_(int *, float *, float *); |
| 121 |
|
| 122 |
//void getpkbgnversion_(char *, int); |
| 123 |
|
| 124 |
|
| 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 |
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 |
char *ffttmp, *fftfull, *infft; |
| 152 |
float cdetrend; |
| 153 |
char *syscall; |
| 154 |
|
| 155 |
char *gapfile, *logfile, *paramfile, *par1file, *resfile, *crossfile, *crossnfile; |
| 156 |
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 |
float woodb1, woodb2; |
| 168 |
|
| 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 |
DRMS_Segment_t *gapseg = NULL; |
| 179 |
DRMS_Array_t *gaparr = NULL; |
| 180 |
DRMS_RecordSet_t *ffttmprecset = NULL; |
| 181 |
DRMS_Record_t *ffttmprec = NULL; |
| 182 |
DRMS_Segment_t *ffttmpseg = NULL; |
| 183 |
FILE *ffttmpfile = NULL; |
| 184 |
char tstartstr[100]; |
| 185 |
char fftfile[DRMS_MAXPATHLEN+1]; |
| 186 |
char dirstr[DRMS_MAXPATHLEN+1]; |
| 187 |
|
| 188 |
TIME tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */ |
| 189 |
double tstart, tstartout, 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 |
float *datarr; |
| 209 |
|
| 210 |
ffttmp = (char *)cmdparams_save_str (&cmdparams, "ffttmp", &newstat); |
| 211 |
logfile = (char *)cmdparams_save_str (&cmdparams, "LOGFILE", &newstat); |
| 212 |
readfft = cmdparams_save_int (&cmdparams, "READFFT", &newstat); |
| 213 |
// savefft = cmdparams_save_int (&cmdparams, "SAVEFFT", &newstat); |
| 214 |
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 |
// lmode = cmdparams_save_int (&cmdparams, "LMODE", &newstat); |
| 232 |
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 |
// tsample = cmdparams_save_float (&cmdparams, "TSAMPLE", &newstat); |
| 242 |
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 |
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 |
woodb1 = cmdparams_save_float (&cmdparams, "WOODB1", &newstat); |
| 268 |
woodb2 = cmdparams_save_float (&cmdparams, "WOODB2", &newstat); |
| 269 |
ctx = cmdparams_save_int (&cmdparams, "CTX", &newstat); |
| 270 |
//printf("ndl = %d\n",ndl); |
| 271 |
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 |
int seginflag = strcmp(kNOTSPECIFIED, segnamein); |
| 276 |
int segoutflag = strcmp(kNOTSPECIFIED, segnameout); |
| 277 |
char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat); |
| 278 |
// 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 |
if (permflag) |
| 283 |
lifetime = DRMS_PERMANENT; |
| 284 |
else |
| 285 |
lifetime = DRMS_TRANSIENT; |
| 286 |
*/ |
| 287 |
|
| 288 |
if (newstat) |
| 289 |
{ |
| 290 |
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 |
{ |
| 296 |
fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr)); |
| 297 |
return 1; |
| 298 |
} |
| 299 |
|
| 300 |
|
| 301 |
int histflag = strncasecmp("none", outseries, 4); |
| 302 |
if (histflag) |
| 303 |
{ |
| 304 |
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 |
histrecnum=set_history(histlink); |
| 320 |
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 |
} |
| 339 |
|
| 340 |
|
| 341 |
/* Set PARAMFILE to PAR1FILE if not set */ |
| 342 |
|
| 343 |
if (strcmp(paramfile,"nofile") == 0) |
| 344 |
paramfile=par1file; |
| 345 |
|
| 346 |
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 |
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 |
if (verbflag) |
| 384 |
printf("crosspath = %s\n", crosspath); |
| 385 |
} |
| 386 |
else |
| 387 |
{ |
| 388 |
crosspath=strdup(crossfile); |
| 389 |
} |
| 390 |
|
| 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 |
return 1; |
| 399 |
} |
| 400 |
if (gaprecset->n == 0) |
| 401 |
{ |
| 402 |
fprintf(stderr,"ERROR: gapfile recordset contains no records\n"); |
| 403 |
return 1; |
| 404 |
} |
| 405 |
if (gaprecset->n > 1) |
| 406 |
{ |
| 407 |
fprintf(stderr,"ERROR: gapfile recordset with more than 1 record not yet supported.\n"); |
| 408 |
return 1; |
| 409 |
} |
| 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 |
return 1; |
| 417 |
} |
| 418 |
rgaps=(float *)(gaparr->data); |
| 419 |
gapsize=gaparr->axis[0]; |
| 420 |
} |
| 421 |
else |
| 422 |
{ |
| 423 |
fits_get_img_size(fitsptr, 1, fdims, &fstatus); |
| 424 |
gapsize=fdims[0]; |
| 425 |
rgaps=(float *)(malloc(gapsize*sizeof(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 |
return 1; |
| 432 |
} |
| 433 |
} |
| 434 |
|
| 435 |
if (verbflag) |
| 436 |
printf("gapfile read, gapsize = %d\n",gapsize); |
| 437 |
|
| 438 |
igaps=(float *)(malloc(gapsize*sizeof(float))); |
| 439 |
gaps=(float *)(malloc(gapsize*sizeof(float))); |
| 440 |
gaps1=(float *)(malloc(gapsize*sizeof(float))); |
| 441 |
gaps2=(float *)(malloc(gapsize*sizeof(float))); |
| 442 |
datr=(float *)(malloc(gapsize*sizeof(float))); |
| 443 |
dati=(float *)(malloc(gapsize*sizeof(float))); |
| 444 |
/* |
| 445 |
help=(float *)(malloc((4*gapsize+30)*sizeof(float))); |
| 446 |
Changed 20061022 */ |
| 447 |
help=(float *)(malloc((4*ndt+30)*sizeof(float))); |
| 448 |
gaps3=(float *)(malloc(ndt*sizeof(float))); |
| 449 |
data3=(float *)(malloc(2*ndt*sizeof(float))); |
| 450 |
|
| 451 |
isign = 1.0; |
| 452 |
if (flipm != 0) |
| 453 |
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 |
if (verbflag) |
| 465 |
printf("%f\n",sum); |
| 466 |
|
| 467 |
cffti_(&ndt, help); |
| 468 |
|
| 469 |
|
| 470 |
DRMS_RecordSet_t *inrecset = drms_open_records(drms_env, inrecquery, &status); |
| 471 |
if (status != DRMS_SUCCESS || inrecset == NULL) |
| 472 |
{ |
| 473 |
fprintf(stderr,"ERROR: problem opening input recordset: status = %d\n",status); |
| 474 |
return 1; |
| 475 |
} |
| 476 |
|
| 477 |
if (inrecset->n == 0) |
| 478 |
{ |
| 479 |
fprintf(stderr,"ERROR: input recordset contains no records\n"); |
| 480 |
return 1; |
| 481 |
} |
| 482 |
|
| 483 |
if (inrecset->n > 1) |
| 484 |
{ |
| 485 |
fprintf(stderr,"ERROR: nrecs > 1 not yet supported.\n"); |
| 486 |
return 1; |
| 487 |
} |
| 488 |
|
| 489 |
|
| 490 |
inrec=inrecset->records[0]; |
| 491 |
int itest; |
| 492 |
char *inchecklist[] = {"T_START", "T_STOP", "LMIN", "LMAX", "T_STEP"}; |
| 493 |
|
| 494 |
for (itest=0; itest < ARRLENGTH(inchecklist); itest++) |
| 495 |
{ |
| 496 |
DRMS_Keyword_t *inkeytest = hcon_lookup_lower(&inrec->keywords, inchecklist[itest]); |
| 497 |
if (inkeytest == NULL) |
| 498 |
{ |
| 499 |
fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]); |
| 500 |
drms_close_records(inrecset, DRMS_FREE_RECORD); |
| 501 |
return 1; |
| 502 |
} |
| 503 |
} |
| 504 |
|
| 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 |
tstartout=tstart+ifirst*tstep; |
| 509 |
|
| 510 |
navail=(tstop-tstart)/tstep; |
| 511 |
if (navail != gapsize) |
| 512 |
{ |
| 513 |
fprintf(stderr, "ERROR: gaps seem incompatible with timeseries, gapsize=%d, navail=%d\n", gapsize, navail); |
| 514 |
return 0; |
| 515 |
} |
| 516 |
|
| 517 |
nread=navail; |
| 518 |
tsample=tstep; |
| 519 |
|
| 520 |
if ((ifirst+ndt) < navail) |
| 521 |
nread=ifirst+ndt; |
| 522 |
|
| 523 |
if (verbflag) |
| 524 |
printf("%d %d\n",navail,nread); |
| 525 |
|
| 526 |
int lmin=drms_getkey_int(inrec, "LMIN", NULL); |
| 527 |
int lmax=drms_getkey_int(inrec, "LMAX", NULL); |
| 528 |
if (lmin != lmax) |
| 529 |
{ |
| 530 |
fprintf(stderr,"ERROR: lmin != lmax not yet supported.\n"); |
| 531 |
return 0; |
| 532 |
} |
| 533 |
lmode=lmin; |
| 534 |
|
| 535 |
char *tag = drms_getkey_string(inrec, "TAG", &status); |
| 536 |
if (status != DRMS_SUCCESS) |
| 537 |
tag=strdup("none"); |
| 538 |
if (readfft == 0) |
| 539 |
{ |
| 540 |
ffttmprec = drms_create_record(drms_env, ffttmp, DRMS_PERMANENT, &status); |
| 541 |
|
| 542 |
if (status != DRMS_SUCCESS) |
| 543 |
{ |
| 544 |
fprintf(stderr,"ERROR: couldn't create a record in ffttmp dataseries %s, status = %d\n", ffttmp, status); |
| 545 |
return 0; |
| 546 |
} |
| 547 |
|
| 548 |
drms_copykeys(ffttmprec, inrec, 0, kDRMS_KeyClass_Explicit); |
| 549 |
drms_setkey_time(ffttmprec, "T_START", tstartout); |
| 550 |
drms_setkey_int(ffttmprec, "LMODE", lmode); |
| 551 |
drms_setkey_int(ffttmprec, "NDT", ndt); |
| 552 |
tnow = (double)time(NULL); |
| 553 |
tnow += UNIX_epoch; |
| 554 |
drms_setkey_time(ffttmprec, "DATE", tnow); |
| 555 |
|
| 556 |
ffttmpseg = drms_segment_lookupnum(ffttmprec, 0); |
| 557 |
|
| 558 |
sprintf(fftfile, "fft%d", lmode); |
| 559 |
ffttmpfile = drms_segment_fopen(ffttmpseg, fftfile, 0, &status); |
| 560 |
if (status != DRMS_SUCCESS) |
| 561 |
{ |
| 562 |
fprintf(stderr,"ERROR: couldn't open a file to write ffttmp record, status = %d\n", status); |
| 563 |
return 0; |
| 564 |
} |
| 565 |
|
| 566 |
drms_segment_filename(ffttmpseg,fftfile); |
| 567 |
|
| 568 |
} |
| 569 |
else |
| 570 |
{ |
| 571 |
char *ffttmpquery = malloc(100); |
| 572 |
sprint_time(tstartstr, tstartout, "TAI", 0); |
| 573 |
sprintf(ffttmpquery, "%s[%s][%d][%d][%s]", ffttmp,tstartstr,lmode,ndt,tag); |
| 574 |
//printf("ffttmpquery = %s\n", ffttmpquery); |
| 575 |
ffttmprecset = drms_open_records(drms_env, ffttmpquery, &status); |
| 576 |
if (status != DRMS_SUCCESS || ffttmprecset == NULL || ffttmprecset->n == 0) |
| 577 |
{ |
| 578 |
fprintf(stderr,"ERROR: couldn't open ffttmp record %s, status = %d\n", ffttmpquery, status); |
| 579 |
return 0; |
| 580 |
} |
| 581 |
|
| 582 |
ffttmprec = ffttmprecset->records[0]; |
| 583 |
ffttmpseg = drms_segment_lookupnum(ffttmprec, 0); |
| 584 |
|
| 585 |
drms_segment_filename(ffttmpseg,fftfile); |
| 586 |
|
| 587 |
} |
| 588 |
|
| 589 |
if (verbflag) |
| 590 |
printf("READFFT=%d,fftfull=%s\n",readfft,fftfile); |
| 591 |
|
| 592 |
//printf("cdetrend=%f,ifill=%d\n",cdetrend,ifill); |
| 593 |
|
| 594 |
/* Loop over m and do FFT's */ |
| 595 |
if (readfft == 0) |
| 596 |
{ |
| 597 |
if (seginflag) |
| 598 |
segin = drms_segment_lookup(inrec, segnamein); |
| 599 |
else |
| 600 |
segin = drms_segment_lookupnum(inrec, 0); |
| 601 |
if (segin) |
| 602 |
inarr = drms_segment_read(segin, usetype, &status); |
| 603 |
if (status != DRMS_SUCCESS || inarr == NULL || segin == NULL) |
| 604 |
{ |
| 605 |
fprintf(stderr, "problem reading input segment: status = %d\n", status); |
| 606 |
return 0; |
| 607 |
} |
| 608 |
|
| 609 |
datarr=(float *)(inarr->data); |
| 610 |
|
| 611 |
for (m=0; m<= lmode; m++) |
| 612 |
{ |
| 613 |
|
| 614 |
// printf("%d %d\n",ndt,m); |
| 615 |
|
| 616 |
data=datarr + m*2l*gapsize; // Force long computation |
| 617 |
|
| 618 |
for (j=0; j<navail; j++) |
| 619 |
{ |
| 620 |
if (isnan(data[2*j])) |
| 621 |
datr[j]=0.0; |
| 622 |
else |
| 623 |
datr[j]=rgaps[j]*data[2*j]; |
| 624 |
if (isnan(data[2*j+1])) |
| 625 |
dati[j]=0.0; |
| 626 |
else |
| 627 |
dati[j]=igaps[j]*data[2*j+1]; |
| 628 |
} |
| 629 |
|
| 630 |
if (cdetrend != 0.0) |
| 631 |
{ |
| 632 |
detrend(datr,gaps,navail,cdetrend); |
| 633 |
detrend(dati,gaps,navail,cdetrend); |
| 634 |
} |
| 635 |
|
| 636 |
if (ifill == 0) |
| 637 |
{ |
| 638 |
for (i=0; i<navail; i++) |
| 639 |
{ |
| 640 |
gaps1[i]=gaps[i]; |
| 641 |
} |
| 642 |
} |
| 643 |
else |
| 644 |
{ |
| 645 |
fill_gaps(datr,gaps,gaps1,navail); |
| 646 |
if (m !=0 ) |
| 647 |
{ |
| 648 |
fill_gaps(dati,gaps,gaps1,navail); |
| 649 |
} |
| 650 |
} |
| 651 |
|
| 652 |
if (fdiff!=0) |
| 653 |
{ |
| 654 |
for (j=0; j<(navail-1); j++) |
| 655 |
{ |
| 656 |
if ((gaps1[j]==1) & (gaps1[j+1]==1)) |
| 657 |
{ |
| 658 |
gaps2[j]=1.0; |
| 659 |
datr[j]=datr[j+1]-datr[j]; |
| 660 |
dati[j]=dati[j+1]-dati[j]; |
| 661 |
} |
| 662 |
else |
| 663 |
{ |
| 664 |
gaps2[j]=0.0; |
| 665 |
datr[j]=0.0; |
| 666 |
dati[j]=0.0; |
| 667 |
} |
| 668 |
} |
| 669 |
gaps2[navail-1]=0.0; |
| 670 |
datr[navail-1]=0.0; |
| 671 |
dati[navail-1]=0.0; |
| 672 |
} |
| 673 |
else |
| 674 |
{ |
| 675 |
for (j=0; j<navail; j++) |
| 676 |
{ |
| 677 |
gaps2[j]=gaps1[j]; |
| 678 |
} |
| 679 |
} |
| 680 |
|
| 681 |
/* Make complex array. |
| 682 |
Now use data3 instead of reusing data in order to allow ndt>navail. */ |
| 683 |
nx1=ndt; |
| 684 |
if (nread < (ifirst+ndt)) |
| 685 |
nx1=nread-ifirst; |
| 686 |
|
| 687 |
for (j=0; j<nx1; j++) |
| 688 |
{ |
| 689 |
gaps3[j]=gaps2[j+ifirst]; |
| 690 |
data3[2*j]=datr[j+ifirst]; |
| 691 |
data3[2*j+1]=dati[j+ifirst]; |
| 692 |
} |
| 693 |
for (j = nx1; j < ndt; j++) |
| 694 |
{ |
| 695 |
gaps3[j] = 0; |
| 696 |
data3[2*j] = 0; |
| 697 |
data3[2*j+1] = 0; |
| 698 |
} |
| 699 |
|
| 700 |
if ((ifix == 1) && (m == 0)) |
| 701 |
{ |
| 702 |
c = sqrt(2.0); |
| 703 |
for (j = 0; j < ndt; j++) |
| 704 |
{ |
| 705 |
data3[2*j] *= c; |
| 706 |
data3[2*j+1] *= c; |
| 707 |
} |
| 708 |
} |
| 709 |
|
| 710 |
cfftb_ (&ndt, data3, help); |
| 711 |
|
| 712 |
fwrite ( (char*)data3, 8*ndt, 1, ffttmpfile); |
| 713 |
|
| 714 |
|
| 715 |
} /* End of m loop */ |
| 716 |
|
| 717 |
/* Pad the file for buffered reading in Fortran when NDT is small */ |
| 718 |
fwrite ( (char*)data3, 8*ndt, 1, ffttmpfile); |
| 719 |
drms_segment_fclose(ffttmpfile); |
| 720 |
if (verbflag) |
| 721 |
printf("fft file written\n"); |
| 722 |
} |
| 723 |
else |
| 724 |
{ |
| 725 |
/* Fill dummy time-series to make gaps */ |
| 726 |
for (j=0; j<navail; j++) |
| 727 |
{ |
| 728 |
datr[j]=j; |
| 729 |
} |
| 730 |
if (cdetrend != 0.0) |
| 731 |
{ |
| 732 |
detrend(datr,gaps,navail,cdetrend); |
| 733 |
} |
| 734 |
if (ifill == 0) |
| 735 |
{ |
| 736 |
for (i=0; i<navail; i++) |
| 737 |
{ |
| 738 |
gaps1[i]=gaps[i]; |
| 739 |
} |
| 740 |
} |
| 741 |
else |
| 742 |
{ |
| 743 |
fill_gaps(datr,gaps,gaps1,navail); |
| 744 |
} |
| 745 |
if (fdiff!=0) |
| 746 |
{ |
| 747 |
for (j=0; j<(navail-1); j++) |
| 748 |
{ |
| 749 |
if ((gaps1[j]==1) & (gaps1[j+1]==1)) |
| 750 |
{ |
| 751 |
gaps2[j]=1.0; |
| 752 |
} |
| 753 |
else |
| 754 |
{ |
| 755 |
gaps2[j]=0.0; |
| 756 |
} |
| 757 |
} |
| 758 |
gaps2[navail-1]=0.0; |
| 759 |
} |
| 760 |
else |
| 761 |
{ |
| 762 |
for (j=0; j<navail; j++) |
| 763 |
{ |
| 764 |
gaps2[j]=gaps1[j]; |
| 765 |
} |
| 766 |
} |
| 767 |
nx1=ndt; |
| 768 |
if (nread < (ifirst+ndt)) |
| 769 |
nx1=nread-ifirst; |
| 770 |
for (j=0; j<nx1; j++) |
| 771 |
{ |
| 772 |
gaps3[j]=gaps2[j+ifirst]; |
| 773 |
} |
| 774 |
for (j = nx1; j < ndt; j++) |
| 775 |
{ |
| 776 |
gaps3[j] = 0; |
| 777 |
} |
| 778 |
} |
| 779 |
|
| 780 |
drms_free_array(inarr); |
| 781 |
if (inrecset) |
| 782 |
{ |
| 783 |
drms_close_records(inrecset, DRMS_FREE_RECORD); |
| 784 |
} |
| 785 |
if (verbflag) |
| 786 |
printf("input recordset closed\n"); |
| 787 |
|
| 788 |
i = -1; |
| 789 |
for ( ptr = noibstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) |
| 790 |
noib[++i] = atoi(ptr); |
| 791 |
|
| 792 |
i = -1; |
| 793 |
for ( ptr = noimstr; (ptr = strtok(ptr, " ")) != NULL; ptr = NULL ) |
| 794 |
noim[++i] = atoi(ptr); |
| 795 |
|
| 796 |
loglen = strlen(logfile); |
| 797 |
fftlen = strlen(fftfile); |
| 798 |
paramlen = strlen(paramfile); |
| 799 |
par1len = strlen(par1file); |
| 800 |
reslen = strlen(resfile); |
| 801 |
crosslen = strlen(crosspath); |
| 802 |
crossnlen = strlen(crossnfile); |
| 803 |
printf("crosspath=%s, crosslen=%d\n",crosspath,crosslen); |
| 804 |
printf("ndl = %d\n",ndl); |
| 805 |
/* Call Fortran program */ |
| 806 |
sub24_(logfile,fftfile, gaps3, &ipar, |
| 807 |
paramfile, par1file, resfile, |
| 808 |
crosspath, &ndeltal, &noisetype, crossnfile, &nnoib, noib, &nnoim, noim, |
| 809 |
&csubm, &ifollow, &anoise, &aother, &lmode, &iodd, &iasym, &ict, &ctx, |
| 810 |
&iwood, &woodb1, &woodb2, &nmin, &nmax, &fmin, &fmax, |
| 811 |
&imfreq, &icase, &nacoeff, &zero, &ndt, &tsample, &idf, &ndf0, &ndf1, |
| 812 |
&cdf, &ndm, &ndl, |
| 813 |
&xsafe, &xsafe1, &delmin, &maxchi, &maxgrad, &maxcof, ®sz, |
| 814 |
&nlobes, &nfour, |
| 815 |
&gongflag, &printfit, |
| 816 |
loglen, fftlen, paramlen, par1len, reslen, crosslen, crossnlen); |
| 817 |
|
| 818 |
if (readfft == 0) |
| 819 |
drms_close_record(ffttmprec,DRMS_INSERT_RECORD); |
| 820 |
else |
| 821 |
drms_close_records(ffttmprecset,DRMS_FREE_RECORD); |
| 822 |
|
| 823 |
printf("module %s successful completion\n", cmdparams.argv[0]); |
| 824 |
return 0; |
| 825 |
} |
| 826 |
|
| 827 |
|
| 828 |
void detrend( |
| 829 |
float *data, |
| 830 |
float *gaps, |
| 831 |
int length, |
| 832 |
float cdetrend |
| 833 |
) |
| 834 |
{ |
| 835 |
#define chunksz 1440 |
| 836 |
#define maxpols 50 |
| 837 |
|
| 838 |
extern float sdot_(int *, float *, int *, float *, int *); |
| 839 |
extern float saxpy_(int *, float *, float *, int *, float *, int *); |
| 840 |
extern float sscal_(int *, float *, float *, int *); |
| 841 |
|
| 842 |
int i,j,k,l,n_good,ndegree; |
| 843 |
int goodlist[chunksz]; |
| 844 |
float pols[maxpols][chunksz]; |
| 845 |
float a,cg,rms2,sum,x[chunksz]; |
| 846 |
float *help; |
| 847 |
int one=1; |
| 848 |
|
| 849 |
for (i=0; i<length ; i=i+chunksz) { |
| 850 |
help=data+i; |
| 851 |
n_good=0; |
| 852 |
cg=0.0; |
| 853 |
for (j=i; j<i+chunksz; j++) { |
| 854 |
if (gaps[j] != 0) { |
| 855 |
goodlist[n_good]=j-i; |
| 856 |
cg=cg+j-i; |
| 857 |
n_good++; |
| 858 |
} |
| 859 |
} |
| 860 |
if (n_good != 0) { |
| 861 |
cg=cg/n_good; |
| 862 |
rms2=0.0; |
| 863 |
for (j=0; j<n_good; j++) {rms2=rms2+(goodlist[j]-cg)*(goodlist[j]-cg);} |
| 864 |
ndegree = floor((goodlist[n_good-1]-goodlist[0])/cdetrend)+2; |
| 865 |
for (k=0; k<n_good; k++) { |
| 866 |
x[k]=(goodlist[k]-cg)/sqrt(rms2); |
| 867 |
} |
| 868 |
for (k=0; k<n_good; k++) {pols[0][k]=1.0/sqrt(n_good);} |
| 869 |
for (k=0; k<n_good; k++) {pols[1][k]=x[k];} |
| 870 |
for (j=2; j<ndegree; j++) { |
| 871 |
for (k=0; k<n_good; k++) { |
| 872 |
pols[j][k]=(2.0*j)/j*x[k]*pols[j-1][k]-(j-1.0)/j*pols[j-2][k]; |
| 873 |
} |
| 874 |
for (l=0; l<j; l++) { |
| 875 |
a=-sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one); |
| 876 |
saxpy_(&n_good,&a,&pols[l][0],&one,&pols[j][0],&one); |
| 877 |
} |
| 878 |
a=1.0/sqrt(sdot_(&n_good,&pols[j][0],&one,&pols[l][0],&one)); |
| 879 |
sscal_(&n_good,&a,&pols[j][0],&one); |
| 880 |
} |
| 881 |
for (j=0; j<ndegree; j++) { |
| 882 |
sum=0.0; |
| 883 |
for (k=0; k<n_good; k++) { |
| 884 |
sum=sum+pols[j][k]*help[goodlist[k]]; |
| 885 |
} |
| 886 |
for (k=0; k<n_good; k++) { |
| 887 |
help[goodlist[k]] -= sum*pols[j][k]; |
| 888 |
} |
| 889 |
} |
| 890 |
} /* if n_good */ |
| 891 |
|
| 892 |
} /* for i */ |
| 893 |
|
| 894 |
} |
| 895 |
|
| 896 |
void memcof(float *, int, int, float *, float *); |
| 897 |
void fixrts(float *, int); |
| 898 |
void predic(float *, int, float *, int, float *, int); |
| 899 |
|
| 900 |
void fill_gaps( |
| 901 |
float *data, |
| 902 |
float *gaps, |
| 903 |
float *newgaps, |
| 904 |
int length |
| 905 |
) |
| 906 |
{ |
| 907 |
#include "nr.h" |
| 908 |
#include "nrutil.h" |
| 909 |
#define chunksz 1440 |
| 910 |
#define maxpoles 10 |
| 911 |
#define maxgap 5 |
| 912 |
|
| 913 |
int i,j,k; |
| 914 |
int n_good, gap_start, gap_end, gap_size, type; |
| 915 |
float *help,good_points[chunksz]; |
| 916 |
int goodlist[chunksz]; |
| 917 |
int npoles; |
| 918 |
float pm,cof[maxpoles]; |
| 919 |
float input[maxpoles]; |
| 920 |
float f_predic[maxgap],b_predic[maxgap]; |
| 921 |
|
| 922 |
npoles=6; |
| 923 |
for (i=0; i<length ; i=i+chunksz) { |
| 924 |
help=data+i; |
| 925 |
n_good=0; |
| 926 |
for (j=i; j<i+chunksz; j++) { |
| 927 |
if (gaps[j] != 0) { |
| 928 |
goodlist[n_good]=j-i; |
| 929 |
good_points[n_good]=data[j]; |
| 930 |
n_good++; |
| 931 |
newgaps[j]=gaps[j]; |
| 932 |
} |
| 933 |
} |
| 934 |
if (n_good != 0) { |
| 935 |
memcof(good_points-1,n_good,npoles,&pm,cof-1); |
| 936 |
fixrts(cof-1,npoles); |
| 937 |
for (j=1; j<n_good; j++) { |
| 938 |
if ((goodlist[j]-goodlist[j-1]) > 1) { |
| 939 |
gap_start=goodlist[j-1]+1; |
| 940 |
gap_end=goodlist[j]-1; |
| 941 |
gap_size=gap_end-gap_start+1; |
| 942 |
if (gap_size <= maxgap) { |
| 943 |
type=0; |
| 944 |
/* Forward prediction */ |
| 945 |
/* Check that there are npoles good points before gap */ |
| 946 |
if ((j >= npoles) && (goodlist[j-npoles] == goodlist[j-1]-npoles+1)) { |
| 947 |
for (k=0; k<npoles; k++) input[k]=help[gap_start-npoles+k]; |
| 948 |
predic(input-1,npoles,cof-1,npoles,f_predic-1,gap_size); |
| 949 |
type +=1; |
| 950 |
} |
| 951 |
/* Backwards prediction */ |
| 952 |
/* Check that there are npoles good points after gap */ |
| 953 |
if ((j+npoles <= n_good) && (goodlist[j+npoles-1] == goodlist[j]+npoles-1)) { |
| 954 |
for (k=0; k<npoles; k++) input[k]=help[gap_end+npoles-k]; |
| 955 |
/* for (k=0; k<npoles; k++) printf("%d %d %f\n",gap_end,k,help[gap_end+npoles-k]); */ |
| 956 |
predic(input-1,npoles,cof-1,npoles,b_predic-1,gap_size); |
| 957 |
type +=2; |
| 958 |
} |
| 959 |
if (type == 3) { |
| 960 |
/* Average predictions */ |
| 961 |
for (k=0; k<gap_size; k++) |
| 962 |
help[gap_start+k] = (f_predic[k]+b_predic[gap_size-k-1])/2.; |
| 963 |
} |
| 964 |
if (type == 1) { |
| 965 |
/* Use forward prediction */ |
| 966 |
for (k=0; k<gap_size; k++) |
| 967 |
help[gap_start+k] = f_predic[k]; |
| 968 |
} |
| 969 |
if (type == 2) { |
| 970 |
/* Use backwards prediction */ |
| 971 |
for (k=0; k<gap_size; k++) |
| 972 |
help[gap_start+k] = b_predic[gap_size-k-1]; |
| 973 |
} |
| 974 |
if (type !=0 ) { |
| 975 |
for (k=0; k<gap_size; k++) |
| 976 |
newgaps[gap_start+k+i]=1; |
| 977 |
} |
| 978 |
} |
| 979 |
} |
| 980 |
} /* for j */ |
| 981 |
} /* if n_good */ |
| 982 |
} /* for i */ |
| 983 |
|
| 984 |
} |
