ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/JSOC/proj/lev0/apps/build_lev1_hmi.c
Revision: 1.18
Committed: Mon Aug 18 20:59:30 2014 UTC (9 years, 1 month ago) by prodtest
Content type: text/plain
Branch: MAIN
CVS Tags: Ver_8-8, Ver_8-6, Ver_8-7, Ver_8-10
Changes since 1.17: +7 -0 lines
Log Message:
Added quality check for camera minimum value error

File Contents

# Content
1 /*-----------------------------------------------------------------------------
2 * cvs/JSOC/proj/lev1/apps/build_lev1.c
3 *-----------------------------------------------------------------------------
4 *
5 * This is a module that runs with DRMS and processes lev0
6 * filtergrams to lev1.
7 * It is scheduled by build_lev1_mgr either by qsub to the cluster
8 * or by fork to run on the local machine. It is called with
9 * 12 (NUMRECLEV1 defined in lev0lev1.h) lev0 images at a time.
10 * See the -h nice_intro() for details of the calling options.
11 */
12
13 #include <jsoc_main.h>
14 #include <cmdparams.h>
15 #include <drms.h>
16 #include <drms_names.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <ctype.h>
20 #include <strings.h>
21 #include <sys/types.h>
22 #include <sys/time.h>
23 #include <sys/stat.h> //for umask(2)
24 #include <unistd.h> //for alarm(2) among other things...
25 #include <printk.h>
26 #include <astro.h>
27 #include <fresize.h>
28 #include <gapfill.h>
29 //#include </home/jsoc/include/fftw3.h>
30 #include "fftw3.h"
31 #include "imgdecode.h"
32 #include "lev0lev1.h"
33 #include "quallev1.h"
34 #include "limb_fit.h"
35 #include "cosmic_ray.h"
36 #include "get_pointing_info.c"
37
38 //default in and out data series
39 //#define LEV0SERIESNAMEHMI "hmi.lev0e"
40 //#define LEV0SERIESNAMEAIA "aia.lev0e"
41 #define LEV0SERIESNAMEHMI "su_production.lev0f_hmi"
42 #define LEV0SERIESNAMEAIA "su_production.lev0f_aia"
43 //#define LEV1SERIESNAMEHMI "su_production.hmi_lev1e" //temp test case
44 #define LEV1SERIESNAMEHMI "hmi.lev1"
45 #define LEV1SERIESNAMEAIA "su_production.aia_lev1e" //temp test case
46 //#define DSFFNAME "su_richard.flatfield" //temp test case
47 //#define DSFFNAMEHMI "su_production.hmi_flatfield" //temp test case
48 #define DSFFNAMEHMI "hmi.flatfield"
49 //#define DSFFNAMEAIA "su_production.aia_flatfield" //temp test case
50 //#define DSFFNAMEAIA "aia_test.flatfield" //temp test case
51 #define DSFFNAMEAIA "aia.flatfield"
52
53 #define kDpath "dpath"
54 #define kDpathDef "/home/jsoc/cvs/Development/JSOC"
55
56 #define LEV1LOG_BASEDIR "/usr/local/logs/lev1"
57 #define H1LOGFILE "/usr/local/logs/lev1/build_lev1_hmi.%s.log"
58 #define NUMTIMERS 8 //number of seperate timers avail
59 #define NOTSPECIFIED "***NOTSPECIFIED***"
60 #define LOGTEST 0
61 #define CAL_HCFTID 17 //image is cal mode
62 #define STOP_FILE "/usr/local/logs/lev1/build_mgr_stop_hmi"
63
64 int compare_rptr(const void *a, const void *b);
65 static TIME SDO_to_DRMS_time(int sdo_s, int sdo_ss);
66
67 // List of default parameter values.
68 ModuleArgs_t module_args[] = {
69 {ARG_STRING, "instru", NOTSPECIFIED, "instrument. either hmi or aia"},
70 {ARG_STRING, "mode", NOTSPECIFIED, "either recnum or fsn"},
71 {ARG_STRING, "dsin", NOTSPECIFIED, "dataset of lev0 filtergrams"},
72 {ARG_STRING, "dsout", NOTSPECIFIED, "dataset of lev1 output"},
73 {ARG_STRING, "dsff", NOTSPECIFIED, "dataset of darks flats bads"},
74 {ARG_STRING, "dsaiabad", NOTSPECIFIED, "dataset of AIA bad pixels"},
75 {ARG_STRING, "logfile", NOTSPECIFIED, "optional log file name. Will create one if not given"},
76 {ARG_INTS, "brec", "0", "first lev0 rec# to process. 0=error must be given by build_lev1_mgr"},
77 {ARG_INTS, "erec", "0", "last lev0 rec# to process. 0=error must be given by build_lev1_mgr"},
78 {ARG_INTS, "bfsn", "0", "first lev0 fsn# to process. 0=error must be given by build_lev1_mgr"},
79 {ARG_INTS, "efsn", "0", "last lev0 fsn# to process. 0=error must be given by build_lev1_mgr"},
80 {ARG_INTS, "quicklook", "1", "1=quick look, 0 = definitive mode"},
81 {ARG_FLAG, "v", "0", "verbose flag"},
82 {ARG_FLAG, "h", "0", "help flag"},
83 {ARG_FLAG, "r", "0", "restart flag"},
84 {ARG_STRING, kDpath, kDpathDef, "path to the source code tree (the JSOC root directory)"},
85 {ARG_END}
86 };
87
88 CmdParams_t cmdparams;
89 // Module name presented to DRMS.
90 char *module_name = "build_lev1_hmi";
91
92 FILE *h1logfp; // fp for h1 ouput log for this run
93 static struct stat stbuf;
94 //static IMG Image0, Image1;
95 //static IMG *Img0 = &Image0;
96 //static IMG *Img1 = &Image1;
97 static LEV0LEV1 lev0lev1;
98 static LEV0LEV1 *l0l1 = &lev0lev1;
99 //static CCSDS_Packet_t *Hk;
100 static DRMS_Record_t *rs;
101 static DRMS_Record_t *rs0, *rs1, *rsff, *rsbad_pix, *rec_bad_aia;
102 static DRMS_Record_t *rptr;
103 static DRMS_Segment_t *segment;
104 static DRMS_Segment_t *segmentff;
105 static DRMS_Segment_t *darkseg;
106 static DRMS_Segment_t *badseg;
107 static DRMS_Segment_t *badoutpixseg;
108 static DRMS_Segment_t *spikeseg;
109 static DRMS_Array_t *segArray;
110 static DRMS_RecordSet_t *rset0, *rset1, *rsetff, *rsbad_aia;
111 static DRMS_Array_t *Array0;
112 static DRMS_Array_t *Arrayff;
113 static DRMS_Array_t *ArrayDark;
114 static DRMS_Array_t *ArrayBad;
115 static DRMS_Array_t *ArraySpike;
116 static TIME sdo_epoch;
117 static PTINFO *ptinfo = NULL;
118 static PTINFO ptdata;
119 static char bld_vers[16];
120 static char datestr[32];
121 static char open_dsname[256];
122 static char dsffname[128];
123 static char open_aiabad_dsname[256];
124 static char path[DRMS_MAXPATHLEN], bad_pix_path[DRMS_MAXPATHLEN];
125 static char bad_aia_path[DRMS_MAXPATHLEN];
126 static char rs1_path[DRMS_MAXPATHLEN];
127 static struct timeval first[NUMTIMERS], second[NUMTIMERS];
128 static char *orbseries = "sdo.fds_orbit_vectors";
129 //static char *orbseries = "sdo_ground.fds_orbit_vectors";
130
131 static int nspikes, respike, fid, aiftsid, *oldvalues, *spikelocs, *newvalues;
132 static int hcftid, aiagp6;
133 static short aifcps;
134 double aiascale = 1.0;
135
136 IORBIT_Info_t *IOinfo = NULL;
137 IORBIT_Info_t IOdata;
138 LIBASTRO_Error_t IOstatus;
139 unsigned int fsnarray[NUMRECLEV1];
140 unsigned int fsnx = 0;
141 //short data1[MAXPIXELS];
142 //int data1[MAXPIXELS];
143 float data1[MAXPIXELS]; //floats for HMI
144 float ftmp;
145 int data1A[MAXPIXELS]; //ints for AIA
146 int array_cosmic[16777216]; //4096x4096
147 double tgttimes[NUMRECLEV1];
148
149 long long brec, erec, bfsn, efsn; //begin and end lev0 rec/fsn. must be same data type
150 long long bnumx, enumx; //has either the brec/erec of bfsn/efsn pair accord to mode
151 int verbose;
152 int hmiaiaflg = 0; //0=hmi, 1=aia
153 int modeflg = 0; //0=fsn, 1=recnum
154 int imagecnt = 0; // num of images since last commit
155 int restartflg = 0; // set when build_lev1 is called for a restart
156 int abortflg = 0;
157 int sigalrmflg = 0; // set on signal so prog will know
158 int ignoresigalrmflg = 0; // set after a close_image()
159 int quicklook;
160 //global quality flags
161 int flatmiss[NUMRECLEV1];
162 int orbmiss[NUMRECLEV1];
163 int asdmiss[NUMRECLEV1];
164 int mpdmiss[NUMRECLEV1];
165 int limbmiss[NUMRECLEV1];
166 int noimage[NUMRECLEV1];
167 int missflg[NUMRECLEV1];
168
169 char logname[128];
170 char argdsin[128], argdsout[128], arglogfile[128], arginstru[80];
171 char argbx[80], argex[80], argquick[80], argmode[80], argdsaiabad[80];
172 char timetag[32];
173 char tlmseriesname[128]; // e.g. hmi.tlm
174 char lev0seriesname[128]; // e.g. hmi.lev0
175 char *username; // from getenv("USER")
176 char *logfile; // optional log name passed in
177 char *instru; // instument. hmi or aia
178 char *mode; // given mode. recnum or fsn
179 char *dsin; // lev0 input dataset
180 char *dsout; // lev1 output dataset
181 char *dsff; // flat field dataset
182 char *dsaiabad; // AIA bad pixel dataset (kludge around bad ff)
183
184 //!!TEMP
185 typedef struct {
186 float rsun_lf;
187 float x0_lf;
188 float y0_lf;
189 } LIMB_SOMETHING;
190
191
192
193 int get_nspikes() { return nspikes; }
194 int get_respike(void) { return respike; }
195 int *get_spikelocs() { return spikelocs; }
196 int *get_oldvalues() { return oldvalues; }
197 int *get_newvalues() { return newvalues; }
198 void set_nspikes(int new_nspikes) { nspikes = new_nspikes; }
199 void set_spikelocs(int *new_spikelocs) { spikelocs = new_spikelocs; }
200 void set_oldvalues(int *new_oldvalues) { oldvalues = new_oldvalues; }
201 void set_newvalues(int *new_newvalues) { newvalues = new_newvalues; }
202
203 //Set the QUALITY keyword for lev1
204 void do_quallev1(DRMS_Record_t *rs0, DRMS_Record_t *rs1, int inx, unsigned int fsn)
205 {
206 int quallev1 = 0;
207 int rstatus, quallev0, qual_select;
208 char *pchar;
209
210 quallev1 = missflg[inx];
211 if(flatmiss[inx]) quallev1 = quallev1 | Q_NOFLAT;
212 if(orbmiss[inx]) quallev1 = quallev1 | Q_NOORB;
213 if(limbmiss[inx]) quallev1 = quallev1 | Q_NOLIMB;
214 if(asdmiss[inx]) {
215 quallev1 = quallev1 | Q_NOASD;
216 drms_setkey_string(rs1, "ASD_REC", DRMS_MISSING_STRING);
217 }
218 if(mpdmiss[inx]) quallev1 = quallev1 | Q_NOMPD;
219 if(noimage[inx]) quallev1 = quallev1 | Q_MISSALL;
220 if(ptinfo) {
221 ptdata = ptinfo[inx];
222 if(strcmp(ptdata.acs_mode, "SCIENCE")) { //not in science pointing mode
223 quallev1 = quallev1 | Q_NOACS_SCI;
224 }
225 if(strcmp(ptdata.acs_eclp, "NO")) { //in eclipse
226 quallev1 = quallev1 | Q_ACS_ECLP;
227 }
228 if(strcmp(ptdata.acs_sunp, "YES")) { //not in sun presence flag
229 quallev1 = quallev1 | Q_ACS_SUNP;
230 }
231 if(strcmp(ptdata.acs_safe, "NO")) { //safemode set
232 quallev1 = quallev1 | Q_ACS_SAFE;
233 }
234 }
235 pchar = drms_getkey_string(rs0, "IMG_TYPE", &rstatus);
236 if(rstatus) {
237 printk("ERROR: in drms_getkey_string(IMG_TYPE) fsn=%u\n", fsn);
238 }
239 else {
240 if(strcmp(pchar, "LIGHT")) { //dark image
241 quallev1 = quallev1 | Q_IMG_TYPE;
242 }
243 }
244 if(hmiaiaflg) { //aia
245 pchar = drms_getkey_string(rs0, "AISTATE", &rstatus);
246 if(aiftsid >= 0xc000) quallev1 = quallev1 | Q_CAL_IMG;
247 if((aifcps <= -20) ||(aifcps >=100)) quallev1 = quallev1 | Q_AIA_FOOR;
248 if(aiagp6 != 0) quallev1 = quallev1 | Q_AIA_REGF;
249 }
250 else { //hmi
251 pchar = drms_getkey_string(rs0, "HWLTNSET", &rstatus);
252 if((fid >= 1) && (fid <=9999)) {
253 quallev1 = quallev1 | Q_CAL_IMG; //cal image
254 }
255 if(hcftid == CAL_HCFTID) quallev1 = quallev1 | Q_CALM_IMG; //hmi cal mode
256 }
257 if(rstatus) {
258 printk("ERROR: in drms_getkey_string(HWLTNSET or AISTATE) fsn=%u\n", fsn);
259 }
260 else {
261 if(!strcmp(pchar, "OPEN")) { //ISS loop open
262 quallev1 = quallev1 | Q_LOOP_OPEN;
263 }
264 if(!hmiaiaflg) { //for HMI thermal recovery and lunar transit
265 quallev0 = drms_getkey_int(rs0, "QUALITY", &rstatus);
266 qual_select = (quallev0 >> 28) & 0x3;
267 switch(qual_select) {
268 case 1:
269 quallev1 = quallev1 | Q_THERM_RECOV;
270 break;
271 case 2:
272 if(!strcmp(pchar, "CLOSED"))
273 quallev1 = quallev1 | Q_THERM_RECOV;
274 break;
275 case 3:
276 quallev1 = quallev1 | Q_LUNAR_TRAN;
277 break;
278 }
279 }
280 }
281
282 if(quicklook) {
283 quallev1 = quallev1 | Q_NRT;
284 }
285
286 //for camera minimum value error
287 quallev0 = drms_getkey_int(rs0, "QUALITY", &rstatus);
288 if ( (quallev0 & 0x8000) != 0) {
289 quallev1 = quallev1 | Q_CAM_ANOM1;
290 }
291
292 drms_setkey_int(rs1, "QUALITY", quallev1);
293 drms_setkey_string(rs1, "BLD_VERS", bld_vers); //build vers to every record
294 }
295
296 int nice_intro ()
297 {
298 int usage = cmdparams_get_int (&cmdparams, "h", NULL);
299 if (usage)
300 {
301 printf ("Usage:\nbuild_lev1_hmi [-vhr] "
302 "mode=<recnum|fsn> instru=<hmi|aia> dsin=<lev0> dsout=<lev1>\n"
303 "brec=<rec#>|bfsn=<fsn#> erec=<rec#>|efsn=<fsn#>\n"
304 "quicklook=<0|1> [logfile=<file>]\n"
305 " -h: help - show this message then exit\n"
306 " -v: verbose\n"
307 " -r: restart. only used when we restart our selves periodically\n"
308 "mode= recnum: brec and erec have the record # range to process \n"
309 " fsn: bfsn and efsn have the fsn # range to process\n"
310 " For safety, the mode and arg name used must be consistent\n"
311 "instru= instrument. must be 'hmi' or 'aia'\n"
312 "dsin= data set name of lev0 input\n"
313 " default hmi=hmi.lev0e aia=aia.lev0e\n"
314 "dsout= data set name of lev1 output\n"
315 " default hmi=su_production.hmi_lev1e aia=su_production.aia_lev1e\n"
316 "dsff= data set name of AIA flat field series\n"
317 " default aia_test.flatfield\n"
318 "dsaiabad= AIA bad pixel series, default is to use BPL from dsff,\n"
319 "but bugs need to be fixed\n"
320 "brec= first lev0 rec# to process. 0=error must be given by build_lev1_mgr\n"
321 "erec= last lev0 rec# to process. 0=error must be given by build_lev1_mgr\n"
322 "bfsn= first fsn# to process. 0=error must be given by build_lev1_mgr\n"
323 "efsn= last fsn# to process. 0=error must be given by build_lev1_mgr\n"
324 "quicklook= 1 = quicklook mode, 0 = definitive mode\n"
325 "logfile= optional log file name. If not given uses:\n"
326 " /usr/local/logs/lev1/build_lev1_hmi.<time_stamp>.log\n");
327 return(1);
328 }
329 verbose = cmdparams_get_int (&cmdparams, "v", NULL);
330 restartflg = cmdparams_get_int (&cmdparams, "r", NULL);
331 return (0);
332 }
333
334 TIME SDO_to_DRMS_time(int sdo_s, int sdo_ss)
335 {
336 static int firstcall = 1;
337 if (firstcall)
338 {
339 firstcall = 0;
340 }
341 /* XXX fix build 3/18/2008, arta */
342 return(sdo_epoch + (TIME)sdo_s + (TIME)(sdo_ss & 0xFFFF)/65536.0);
343 }
344
345 // Setup global datestr[] like: 2008.07.14_08:29:31
346 char *do_datestr() {
347 time_t tval;
348 struct tm *t_ptr;
349
350 tval = time(NULL);
351 t_ptr = localtime(&tval);
352 sprintf(datestr, "%d.%02d.%02d_%02d:%02d:%02d",
353 (t_ptr->tm_year+1900), (t_ptr->tm_mon+1),
354 t_ptr->tm_mday, t_ptr->tm_hour, t_ptr->tm_min, t_ptr->tm_sec);
355 return(datestr);
356 }
357
358 // Returns a time tag like yyyy.mm.dd.hhmmss
359 char *gettimetag()
360 {
361 struct timeval tvalr;
362 struct tm *t_ptr;
363
364 gettimeofday(&tvalr, NULL);
365 t_ptr = localtime((const time_t *)&tvalr);
366 sprintf(timetag, "%04d.%02d.%02d.%02d%02d%02d",
367 (t_ptr->tm_year+1900), (t_ptr->tm_mon+1), t_ptr->tm_mday, t_ptr->tm_hour, t_ptr->tm_min, t_ptr->tm_sec);
368 return(timetag);
369 }
370
371
372 void BeginTimer(int n)
373 {
374 gettimeofday (&first[n], NULL);
375 }
376
377 float EndTimer(int n)
378 {
379 gettimeofday (&second[n], NULL);
380 if (first[n].tv_usec > second[n].tv_usec) {
381 second[n].tv_usec += 1000000;
382 second[n].tv_sec--;
383 }
384 return (float) (second[n].tv_sec-first[n].tv_sec) +
385 (float) (second[n].tv_usec-first[n].tv_usec)/1000000.0;
386 }
387
388 // Outputs the variable format message (re: printf) to the log file.
389 int h1log(const char *fmt, ...)
390 {
391 va_list args;
392 char string[32768];
393
394 va_start(args, fmt);
395 vsprintf(string, fmt, args);
396 if(h1logfp) {
397 fprintf(h1logfp, string);
398 fflush(h1logfp);
399 }
400 else { // couldn't open log
401 printf(string); // also print to stdout
402 fflush(stdout);
403 }
404 va_end(args);
405 return(0);
406 }
407
408 int send_mail(char *fmt, ...)
409 {
410 va_list args;
411 char string[1024], cmd[1024];
412
413 va_start(args, fmt);
414 vsprintf(string, fmt, args);
415 sprintf(cmd, "echo \"%s\" | Mail -s \"build_lev1_hmi mail\" lev0_user", string);
416 system(cmd);
417 va_end(args);
418 return(0);
419 }
420
421 // Got a fatal error.
422 void abortit(int stat)
423 {
424 printk("***Abort in progress ...\n");
425 printk("**Exit build_lev1_hmi w/ status = %d\n", stat);
426 if (h1logfp) fclose(h1logfp);
427 exit(stat);
428 }
429
430 //!!TBD Keh-Cheng
431 int rdout_mode_correct()
432 {
433 return(0);
434 }
435
436 //!!TBD Art
437 int orbit_calc()
438 {
439 return(0);
440 }
441
442 //#include "aia_despike.c"
443 #include "do_flat.c"
444 #include "do_patch.c"
445 #include "get_image_location.c"
446 #include "limb_fit_function.c"
447 #include "cosmic_ray.c"
448 #include "heightformation.c"
449
450 //Called with the range to do. The args are either rec#s or fsn#s accord to modeflg
451 int do_ingest(long long bbrec, long long eerec, const char *dpath)
452 {
453 //FILE *fwt;
454 Image_Location *p_imageloc;
455 Image_Location imageloc[NUMRECLEV1];
456 TIME t_obs0;
457 TIME tobs[NUMRECLEV1];
458 float percentd;
459 float cdelt1, rsun, crpix1, crpix2, crota2;
460 double rsun_lf, x0_lf, y0_lf;
461 int rstatus, dstatus, lstatus, ncnt, fcnt, i, j, k, qualint, nobs;
462 int hshiexp, hcamid, nbad, n_cosmic;
463 int *spikedata, status, axes[2], nbytes;
464 uint32_t missvals, totvals;
465 long long recnum0, recnum1, recnumff;
466 char recrange[128], lev0name[128], flatrec[128];
467 //char tmpname[80];
468
469 if(modeflg) sprintf(recrange, ":#%lld-#%lld", bbrec, eerec);
470 else sprintf(recrange, "%lld-%lld", bbrec, eerec);
471 sprintf(open_dsname, "%s[%s]", dsin, recrange);
472 printk("open_dsname = %s\n", open_dsname);
473 printk("#levnum recnum fsn\n");
474
475 rset0 = drms_open_records(drms_env, open_dsname, &rstatus); //open lev0
476 if(!rset0 || (rset0->n == 0) || rstatus) {
477 printk("Can't do drms_open_records(%s)\n", open_dsname);
478 return(1);
479 }
480 drms_stage_records(rset0, 1, 0);
481 ncnt = rset0->n;
482 rptr = (DRMS_Record_t *)malloc(ncnt * sizeof(DRMS_Record_t));
483 if(rptr == NULL) {
484 printk("Can't malloc() for DRMS_Record_t sort\n");
485 return(1);
486 }
487 //make opened records sequential in mem for sort
488 for(i=0; i < ncnt; i++) {
489 memcpy(&rptr[i], rset0->records[i], sizeof(DRMS_Record_t));
490 }
491 //Must sort to get in ascending t_obs order
492 qsort(rptr, ncnt, sizeof(DRMS_Record_t), &compare_rptr);
493
494 //this loop is for the benefit of the lev1view display to show
495 //info on all lev0 records opened
496 for(i=0; i < ncnt; i++) {
497 flatmiss[i] = 0; //init quality flags
498 orbmiss[i] = 0;
499 limbmiss[i] = 0;
500 asdmiss[i] = 0;
501 mpdmiss[i] = 0;
502 noimage[i] = 0;
503 missflg[i] = 0;
504 rs0 = &rptr[i];
505 recnum0 = rs0->recnum;
506 //must get fsn in case got record from last lev1 record
507 fsnx = drms_getkey_int(rs0, "FSN", &rstatus);
508 fsnarray[i] = fsnx;
509 printk("*0 %u %u\n", recnum0, fsnx);
510 //also set up call for get_pointing_info() and iorbit_getinfo()
511 tobs[i] = drms_getkey_time(rs0, "t_obs", &rstatus);
512 if(rstatus) {
513 printk("Error on drms_getkey_time() fsn=%u. Use DRMS_MISSING_TIME\n",
514 fsnx);
515 tobs[i] = DRMS_MISSING_TIME;
516 }
517 }
518 if(rstatus = get_pointing_info(drms_env, tobs, ncnt, &ptinfo)) {
519 printk("**ERROR: get_pointing_info() status = %d fsn tobs ASD:\n",
520 rstatus);
521 for(j=0; j < ncnt; j++) { //!!TEMP debuf stuff
522 printk("%u %10.5f ", fsnarray[j], tobs[j]);
523 if(ptinfo) {
524 ptdata = ptinfo[j];
525 printk("%s\n", ptdata.asd_rec);
526 }
527 asdmiss[j] = 1; //set for QUALITY for ea image
528 }
529 //return(1); //!!No,press on. New 2/22/2011
530 }
531 if ((IOstatus = iorbit_getinfo(drms_env,
532 orbseries,
533 NULL,
534 IORBIT_Alg_Quadratic,
535 tobs,
536 ncnt,
537 kIORBIT_CacheAction_DontCache,
538 &IOinfo)) != kLIBASTRO_Success)
539 {
540 if(IOstatus == kLIBASTRO_InsufficientData) {
541 printk("***ERROR in iorbit_getinfo: kLIBASTRO_InsufficientData\n");
542 }
543 else {
544 printk("***ERROR in iorbit_getinfo() status=%d\n", IOstatus);
545 }
546 for(j=0; j < ncnt; j++) { //set qual bits
547 orbmiss[j] = 1;
548 }
549 return(1); //new 2/22/2011
550 }
551 rset1 = drms_create_records(drms_env, ncnt, dsout, DRMS_PERMANENT,&dstatus);
552 if(dstatus) {
553 printk("**ERROR: Can't create records for %s\n", dsout);
554 for(j=0; j < ncnt; j++) { //set qual bits
555 noimage[j] = 1;
556 }
557 return(1); //new 2/22/2011
558 }
559 //Now fill in info for call to Carl's get_image_location()
560 for(i=0; i < ncnt; i++) {
561 rs0 = &rptr[i];
562 imageloc[i].tobs = tobs[i];
563 imageloc[i].camera = drms_getkey_int(rs0, "CAMERA", &rstatus);
564 if(rstatus) {
565 printk("ERROR: in drms_getkey_int(CAMERA) fsn=%u\n", fsnarray[i]);
566 }
567 imageloc[i].wavelength = drms_getkey_int(rs0, "WAVELNTH", &rstatus);
568 if(rstatus) {
569 printk("ERROR: in drms_getkey_int(WAVELNTH) fsn=%u\n", fsnarray[i]);
570 }
571 snprintf(imageloc[i].telescope, 10, "%s",
572 drms_getkey_string(rs0, "TELESCOP", &rstatus));
573 if(rstatus) {
574 printk("ERROR: in drms_getkey_string(TELESCOP) fsn=%u\n", fsnarray[i]);
575 }
576 }
577 p_imageloc = imageloc;
578 rstatus = get_image_location(drms_env, ncnt, &p_imageloc);
579 if(rstatus) { //error
580 printk("ERROR: get_image_location() returns status=%d\n", rstatus);
581 for(j=0; j < ncnt; j++) { //set qual bits
582 mpdmiss[i] = 0;
583 }
584 return(1); //new 2/22/2011
585 }
586
587 for(i=0; i < ncnt; i++) { //do for all the sorted lev0 records
588 //StartTimer(2); //!!TEMP
589 rs0 = &rptr[i];
590 recnum0 = rs0->recnum;
591 fsnx = fsnarray[i];
592 sprintf(lev0name, "%s[%u]", dsin, fsnx);
593 if(drms_getkey_int(rs0, "QUALITY", 0) < 0) {
594 printk("Bad QUALITY for %s, no lev1 made\n", lev0name);
595 noimage[i] = 1;
596 //continue; //make an image anyway so the lev1 fsn will exist
597 }
598 segment = drms_segment_lookupnum(rs0, 0);
599 Array0 = drms_segment_read(segment, DRMS_TYPE_SHORT, &rstatus);
600 if(!Array0) {
601 printk("Can't do drms_segment_read() %s status=%d\n",
602 lev0name, rstatus);
603 noimage[i] = 1;
604 return(1); //return until we learn
605 continue;
606 }
607 l0l1->adata0 = (short *)Array0->data; //free at end
608 l0l1->rs0 = rs0;
609 l0l1->recnum0 = recnum0;
610 l0l1->fsn = fsnx;
611 if(hmiaiaflg) { //aia
612 l0l1->dat1.adata1A = &data1A;
613 //l0l1->himgcfid = drms_getkey_int(rs0, "AIFDBID", &rstatus);
614 l0l1->himgcfid = 90; //!!TEMP force uncropped, no overscan
615 aiftsid = drms_getkey_int(rs0, "AIFTSID", &rstatus);
616 aifcps = drms_getkey_short(rs0, "AIFCPS", &rstatus);
617 aiagp6 = drms_getkey_int(rs0, "AIAGP6", &rstatus);
618 }
619 else {
620 // Apply image corruption patch(es)
621 // Patch 1 - crop table corruption Dec 2011 - Jan 2012
622 if (NEED_PATCH1(fsnx))
623 do_patch1(l0l1->adata0);
624 // Patch 2 - camera 1 lookup table corruption 30 March 2014
625 if (NEED_PATCH2(fsnx))
626 do_patch2(l0l1->adata0);
627
628 l0l1->dat1.adata1 = &data1;
629 l0l1->himgcfid = drms_getkey_int(rs0, "HIMGCFID", &rstatus);
630 }
631 if(rstatus) {
632 printk("Can't do drms_getkey_int(HIMGCFID) for fsn %u\n", fsnx);
633 //!!TEMP continue on for testing of AIA lev1
634 l0l1->himgcfid = 104; //!!TEMP force a value for now
635 //return(1); //return until we learn
636 //continue; //maybe cleanup and continue here
637 }
638
639 sprintf(open_dsname, "%s[%u]", dsout, fsnx);
640 rs = rset1->records[i];
641 drms_record_directory(rs, rs1_path, 0);
642 if(!*rs1_path) {
643 printk("***ERROR: No path to segment for %s\n", open_dsname);
644 noimage[i] = 1;
645 continue;
646 }
647 //printf("\npath to lev1 = %s\n", rs1_path); //!!TEMP
648 dstatus = drms_setkey_int(rs, "FSN", fsnx);
649 dstatus = drms_setkey_string(rs, "LEV0SERIES", lev0name);
650 if(!(segment = drms_segment_lookup(rs, "image_lev1"))) {
651 printk("No drms_segment_lookup(rs, image_lev1) for %s\n", open_dsname);
652 noimage[i] = 1;
653 continue;
654 }
655 if(hmiaiaflg) {
656 segArray = drms_array_create(DRMS_TYPE_INT,
657 segment->info->naxis,
658 segment->axis,
659 &data1A,
660 &dstatus);
661 }
662 else {
663 segArray = drms_array_create(DRMS_TYPE_FLOAT,
664 segment->info->naxis,
665 segment->axis,
666 &data1,
667 &dstatus);
668 }
669 //transer the lev0 keywords
670 rstatus = drms_copykeys(rs, rs0, 0, kDRMS_KeyClass_Explicit);
671 if(rstatus != DRMS_SUCCESS) {
672 printk("Error %d in drms_copykeys() for fsn %u\n", fsnx);
673 return(1); //new 2/22/2011
674 continue;
675 }
676 qualint = drms_getkey_int(rs0, "QUALITY", &rstatus);
677 drms_setkey_int(rs, "QUALLEV0", qualint);
678 fid = drms_getkey_int(rs0, "FID", &rstatus);
679
680 drms_setkey_time(rs, "T_OBS", tobs[i]);
681 printk("t_obs for lev0 = %10.5f\n", tobs[i]); //!!TEMP
682 drms_setkey_double(rs, "DATE", CURRENT_SYSTEM_TIME);
683 if(ptinfo) {
684 ptdata = ptinfo[i];
685 drms_setkey_float(rs, "SAT_Y0", ptdata.sat_y0);
686 drms_setkey_float(rs, "SAT_Z0", ptdata.sat_z0);
687 drms_setkey_float(rs, "SAT_ROT", ptdata.sat_rot);
688 drms_setkey_string(rs, "ACS_MODE", ptdata.acs_mode);
689 drms_setkey_string(rs, "ACS_ECLP", ptdata.acs_eclp);
690 drms_setkey_string(rs, "ACS_SUNP", ptdata.acs_sunp);
691 drms_setkey_string(rs, "ACS_SAFE", ptdata.acs_safe);
692 drms_setkey_string(rs, "ASD_REC", ptdata.asd_rec);
693 drms_setkey_string(rs, "ACS_CGT", ptdata.acs_cgt);
694 }
695 if(IOinfo) {
696 IOdata = IOinfo[i];
697 drms_setkey_double(rs, "HAEX_OBS", IOdata.hciX);
698 drms_setkey_double(rs, "HAEY_OBS", IOdata.hciY);
699 drms_setkey_double(rs, "HAEZ_OBS", IOdata.hciZ);
700 drms_setkey_double(rs, "GAEX_OBS", IOdata.gciX);
701 drms_setkey_double(rs, "GAEY_OBS", IOdata.gciY);
702 drms_setkey_double(rs, "GAEZ_OBS", IOdata.gciZ);
703 //drms_setkey_float(rs, "DSUN_OBS", (float)IOdata.dsun_obs);
704 drms_setkey_double(rs, "DSUN_OBS", IOdata.dsun_obs);
705 drms_setkey_double(rs, "OBS_VR", IOdata.obs_vr);
706 drms_setkey_double(rs, "OBS_VW", IOdata.obs_vw);
707 drms_setkey_double(rs, "OBS_VN", IOdata.obs_vn);
708 drms_setkey_double(rs, "RSUN_OBS", IOdata.rsun_obs);
709 drms_setkey_float(rs, "CRLN_OBS", (float)IOdata.crln_obs);
710 drms_setkey_float(rs, "CRLT_OBS", (float)IOdata.crlt_obs);
711 drms_setkey_float(rs, "HGLT_OBS", (float)IOdata.crlt_obs);//!!TEMP
712 drms_setkey_float(rs, "HGLN_OBS", 0.0);//!!TEMP
713 drms_setkey_int(rs, "CAR_ROT", (int)IOdata.car_rot);
714 drms_setkey_string(rs, "ORB_REC", IOdata.orb_rec);
715 }
716 drms_setkey_float(rs, "X0_MP", imageloc[i].x);
717 drms_setkey_float(rs, "Y0_MP", imageloc[i].y);
718 drms_setkey_float(rs, "INST_ROT", imageloc[i].instrot);
719 //drms_setkey_float(rs, "INST_ROT", 180.0); //force this for now
720 drms_setkey_float(rs, "IMSCL_MP", imageloc[i].imscale);
721 drms_setkey_string(rs, "MPO_REC", imageloc[i].mpo_rec);
722
723 int camera = drms_getkey_int(rs0, "CAMERA", &rstatus);
724 if(rstatus) {
725 printk("Can't do drms_getkey_int() for fsn %u\n", fsnx);
726 noimage[i] = 1;
727 goto TEMPSKIP;
728 //return(1);
729 }
730 //Now figure out what flat field to use
731 if(drms_ismissing_time(tobs[i])) {
732 printk("DRMS_MISSING_TIME for fsn=%u. Continue...\n", fsnx);
733 noimage[i] = 1;
734 goto TEMPSKIP;
735 }
736 if(!hmiaiaflg) { //HMI
737 /***OLD - Now just get the pzt flat. the flatfield_version is typically = 0
738 // if(quicklook) {
739 // sprintf(open_dsname, "%s[? t_start=(select max(t_start) from %s where t_start <= %10.5f and t_stop > %10.5f and CAMERA=%d) and CAMERA=%d ?]",
740 // dsffname, dsffname, tobs[i], tobs[i], camera, camera);
741 // }
742 // else {
743 // sprintf(open_dsname, "%s[? t_start <= %10.5f and t_stop > %10.5f and CAMERA=%d and flatfield_version >= 1 ?]",
744 // dsffname, tobs[i], tobs[i], camera);
745 // }
746 *****END OLD*******************************************************/
747 sprintf(open_dsname, "%s[? t_start <= %10.5f and t_stop > %10.5f and CAMERA=%d ?]",
748 dsffname, tobs[i], tobs[i], camera);
749
750 }
751 else { //AIA
752 char *wavstr = drms_getkey_string(rs0, "WAVE_STR", &rstatus);
753 if(rstatus) {
754 printk("Can't do drms_getkey_string() for WAVE_STR\n");
755 return(1);
756 }
757 /***OLD - Now just get the pzt flat. The flatfield_version is typically = 0
758 // if(quicklook) {
759 // sprintf(open_dsname, "%s[? t_start=(select max(t_start) from %s where t_start <= %10.5f and t_stop > %10.5f and WAVE_STR='%s') and WAVE_STR='%s' ?]",
760 // dsffname, dsffname, tobs[i], tobs[i], wavstr, wavstr);
761 // }
762 // else {
763 // sprintf(open_dsname, "%s[? t_start <= %10.5f and t_stop > %10.5f and WAVE_STR='%s' and flatfield_version >= 1 ?]",
764 // dsffname, tobs[i], tobs[i], wavstr);
765 // }
766 *****END OLD**************************************************************/
767 sprintf(open_dsname, "%s[? t_start <= %10.5f and t_stop > %10.5f and WAVE_STR=%s ?]",
768 dsffname, tobs[i], tobs[i], wavstr);
769
770 }
771 //printf("!!TEMP Flat field query: %s\n", open_dsname); //!!TEMP
772 rsetff = drms_open_records(drms_env, open_dsname, &rstatus); //open FF
773 if(!rsetff || (rsetff->n == 0) || rstatus) {
774 printk("Can't do drms_open_records(%s)\n", open_dsname);
775 flatmiss[i] = 1; noimage[i] = 1;
776 goto TEMPSKIP;
777 return(1); //new 2/22/2011
778 }
779 fcnt = rsetff->n;
780 if(fcnt > 1) {
781 printk("More than one FF found for %s?\n", open_dsname);
782 printk("Use last one of %d found\n", fcnt); //!!TEMP
783 //return(1); //!!TBD
784 }
785 //rsff = rsetff->records[0];
786 rsff = rsetff->records[fcnt-1];
787 recnumff = rsff->recnum;
788 sprintf(flatrec, "%s[:#%lld]", dsffname, recnumff);
789 if(dstatus = drms_setkey_string(rs, "FLAT_REC", flatrec )) {
790 printk("**ERROR on drms_setkey_string() for %s\n", flatrec);
791 }
792 drms_record_directory(rsff, path, 1);
793 if(!*path) {
794 printk("***ERROR: No path to segment for %s\n", open_dsname);
795 //goto TEMPSKIP; //!!!TEMP until have good flatfield
796 return(1);
797 }
798 //printf("\npath to FF = %s\n", path); //!!TEMP
799 segmentff = drms_segment_lookup(rsff, "flatfield");
800 Arrayff = drms_segment_read(segmentff, DRMS_TYPE_FLOAT, &rstatus);
801 if(!Arrayff) {
802 printk("Can't do drms_segment_read() for Flat Field status=%d\n",
803 rstatus);
804 return(1);
805 }
806 l0l1->adataff = (float *)Arrayff->data; //!!TBD free at end
807
808 darkseg = drms_segment_lookup(rsff, "DARK");
809 ArrayDark = drms_segment_read(darkseg, DRMS_TYPE_FLOAT, &rstatus);
810 if(!ArrayDark) {
811 printk("Can't do drms_segment_read() for DARK. status=%d\n", rstatus);
812 return(1);
813 }
814 l0l1->adatadark = (float *)ArrayDark->data; //free at end
815
816 if (strcmp(dsaiabad, NOTSPECIFIED)) {
817 char *wavstr = drms_getkey_string(rs0, "WAVE_STR", &rstatus);
818 sprintf(open_aiabad_dsname, "%s[%s][]", dsaiabad, wavstr);
819 rsbad_aia = drms_open_records(drms_env, open_aiabad_dsname, &rstatus);
820 if(!rsbad_aia || rsbad_aia->n == 0 || rstatus)
821 printk("drms_open_records for rsbad_aia failed\n");
822 fcnt = rsbad_aia->n;
823 if(fcnt > 1) {
824 printk("More than one bpl found for %s?\n", open_aiabad_dsname);
825 }
826 rec_bad_aia = rsbad_aia->records[fcnt-1];
827 badseg = drms_segment_lookup(rec_bad_aia, "bad_pixel_list");
828 } else {
829 badseg = drms_segment_lookup(rsff, "BAD_PIXEL");
830 }
831 ArrayBad = drms_segment_read(badseg, DRMS_TYPE_INT, &rstatus);
832 nbad = drms_array_size(ArrayBad)/sizeof(int);
833 if(!ArrayBad) {
834 printk("Can't do drms_segment_read() for BAD_PIXEL. status=%d\n",
835 rstatus);
836 return(1);
837 }
838 l0l1->adatabad = (int *)ArrayBad->data; //free at end
839 //!!NEW 3/9/10
840 /************************************
841 if(!(badoutpixseg = drms_segment_lookup(rs, "BAD_PIXEL"))) BRACKET
842 if(!(badoutpixseg = drms_segment_lookup(rs, "bad_pixel"))) BRACKET
843 ***************************************/
844 if(!hmiaiaflg) { //HMI. !!TBD make the aia & hmi jsd agree
845 if(!(badoutpixseg = drms_segment_lookup(rs, "bad_pixel_list"))) {
846 printk("No drms_segment_lookup(rs, bad_pixel_list) for lev1\n");
847 return(1);
848 }
849 }
850 else {
851 if(!(badoutpixseg = drms_segment_lookup(rs, "bad_pixel"))) {
852 printk("No drms_segment_lookup(rs, bad_pixel) for lev1\n");
853 return(1);
854 }
855 }
856 //move the segment write to after cosmic_rays() updates it
857 /**************************************************************
858 dstatus = drms_segment_write(badoutpixseg, ArrayBad, 0);
859 if (dstatus) {
860 printk("ERROR: drms_segment_write error=%d for lev1 bad_pixel_list\n");
861 //noimage[i] = 1;
862 return(1);
863 }
864 **************************************************************/
865 l0l1->rs1 = rs;
866 l0l1->rsff = rsff;
867 l0l1->recnum1 = rs->recnum;
868 l0l1->darkflag = 0;
869 hshiexp = drms_getkey_int(rs, "HSHIEXP", &rstatus);
870 hcamid = drms_getkey_int(rs, "HCAMID", &rstatus);
871 if(hmiaiaflg) {
872 int aimgshce = drms_getkey_int(rs, "AIMGSHCE", &rstatus);
873 if(aimgshce == 0) l0l1->darkflag = 1;
874 if(rstatus = do_flat_aia(l0l1)) {
875 printk("***ERROR in do_flat_aia() status=%d\n", rstatus);
876 printf("***ERROR in do_flat_aia() status=%d\n", rstatus);
877 flatmiss[i] = 1; noimage[i] = 1;
878 //return(1); //!!TBD what to do?
879 goto FLATERR;
880 }
881 }
882 else {
883 if(hshiexp == 0) l0l1->darkflag = 1;
884 //StartTimer(1); //!!TEMP
885 if(rstatus = do_flat(l0l1)) {
886 printk("***ERROR in do_flat() status=%d\n", rstatus);
887 printf("***ERROR in do_flat() status=%d\n", rstatus);
888 flatmiss[i] = 1; noimage[i] = 1;
889 //return(1); //!!TBD what to do?
890 goto FLATERR;
891 }
892 //ftmp = StopTimer(1);
893 //printf( "\nTime sec for hmi do_flat() = %f\n\n", ftmp );
894 }
895 //sprintf(tmpname, "/tmp/data_lev1.%u", fsnx);
896 //fwt = fopen(tmpname, "w");
897 //int wsize = 4096 * 4096;
898 //fwrite(l0l1->adata1, sizeof(float), wsize, fwt);
899 //fclose(fwt);
900 drms_setkey_float(rs, "OSCNMEAN", l0l1->oscnmean);
901 drms_setkey_float(rs, "OSCNRMS", l0l1->oscnrms);
902 drms_setkey_int(rs, "DATAMIN", l0l1->datamin);
903 drms_setkey_int(rs, "DATAMAX", l0l1->datamax);
904 drms_setkey_int(rs, "DATAMEDN", l0l1->datamedn);
905 drms_setkey_float(rs, "DATAMEAN", l0l1->datamean);
906 drms_setkey_float(rs, "DATARMS", l0l1->data_rms);
907 drms_setkey_float(rs, "DATASKEW", l0l1->dataskew);
908 drms_setkey_float(rs, "DATAKURT", l0l1->datakurt);
909 drms_setkey_int(rs, "DATAVALS", l0l1->datavals);
910 drms_setkey_int(rs, "MISSVALS", l0l1->missvals);
911 missvals = (uint32_t)l0l1->missvals;
912 totvals = (uint32_t)l0l1->datavals + missvals;
913 drms_setkey_int(rs, "TOTVALS", (int)totvals);
914 percentd = (float)((100.0 * (float)l0l1->datavals)/(float)totvals);
915 drms_setkey_float(rs, "PERCENTD", percentd);
916 if(missvals > 0) missflg[i] = missflg[i] | Q_1_MISS0;
917 if(missvals > (uint32_t)(totvals * 0.01))
918 missflg[i] = missflg[i] | Q_1_MISS1;
919 if(missvals > (uint32_t)(totvals * 0.05))
920 missflg[i] = missflg[i] | Q_1_MISS2;
921 if(missvals > (uint32_t)(totvals * 0.25))
922 missflg[i] = missflg[i] | Q_1_MISS3;
923 if(l0l1->datavals == 0)
924 missflg[i] = missflg[i] | Q_MISSALL; //high bit, no data
925 if(hmiaiaflg && nspikes) {
926 if (spikeseg = drms_segment_lookup(rs,"spikes") ) {
927 nbytes = nspikes*sizeof(int);
928 axes[0] = nspikes;
929 axes[1] = 3;
930 spikedata = (int *) malloc (3*nspikes*sizeof(int));
931 ArraySpike = drms_array_create(DRMS_TYPE_INT, 2, axes,
932 (void *) spikedata, &status);
933 memcpy((void *)spikedata, (void *)spikelocs, nbytes);
934 memcpy((void *)(spikedata+nspikes), (void *)oldvalues, nbytes);
935 memcpy((void *)(spikedata+2*nspikes), (void *)newvalues, nbytes);
936 status = drms_segment_write(spikeseg, ArraySpike, 0);
937 drms_free_array(ArraySpike);
938 } // else { printf("spikes segment not found\n"); }
939 }
940 if(!hmiaiaflg) { //only call for hmi
941 //StartTimer(1); //!!TEMP
942 dstatus = cosmic_rays(rs, l0l1->dat1.adata1, l0l1->adatabad, nbad,
943 array_cosmic, &n_cosmic, 4096, 4096);
944 //ftmp = StopTimer(1);
945 //printf( "\nTime sec for cosmic_rays() = %f\n\n", ftmp );
946 if(dstatus) {
947 printk("ERROR: cosmic_rays() error=%d\n", dstatus);
948 noimage[i] = 1;
949 }
950 else {
951 drms_setkey_int(rs, "NBADPERM", nbad);
952 drms_setkey_int(rs, "NBADTOT", n_cosmic);
953 }
954 if(nbad != n_cosmic) { //pretend the bad pix list is like spike
955 nbytes = n_cosmic*sizeof(int);
956 axes[0] = n_cosmic;
957 spikedata = (int *)malloc(nbytes);
958 ArraySpike = drms_array_create(DRMS_TYPE_INT, 1, axes,
959 (void *)spikedata, &status);
960 memcpy((void *)spikedata, (void *)array_cosmic, nbytes);
961 status = drms_segment_write(badoutpixseg, ArraySpike, 0);
962 drms_free_array(ArraySpike);
963 }
964 else {
965 dstatus = drms_segment_write(badoutpixseg, ArrayBad, 0);
966 }
967 if (dstatus) {
968 printk("ERROR: drms_segment_write error=%d for lev1 bad_pixel_list\n",
969 dstatus);
970 noimage[i] = 1;
971 }
972 }
973 FLATERR:
974 drms_close_records(rsetff, DRMS_FREE_RECORD);
975 free(ArrayDark->data);
976 free(Arrayff->data);
977 free(ArrayBad->data);
978
979 TEMPSKIP:
980 if(!hmiaiaflg) {
981 segArray->type = DRMS_TYPE_FLOAT;
982 segArray->bscale = 1.0;
983 segArray->bzero = 0.0;
984 }
985 dstatus = drms_segment_write(segment, segArray, 0);
986 if (dstatus) {
987 printk("ERROR: drms_segment_write error=%d for fsn=%u\n", dstatus,fsnx);
988 noimage[i] = 1;
989 }
990 recnum1 = rs->recnum;
991 printk("*1 %u %u\n", recnum1, fsnx);
992 free(Array0->data);
993
994 x0_lf = DRMS_MISSING_DOUBLE;
995 y0_lf = DRMS_MISSING_DOUBLE;
996 rsun_lf = DRMS_MISSING_DOUBLE;
997 //Don't send calmode or darks to limb_fit()
998 //calmode: HCFTID=17
999 //dark: HSHIEXP=0 and HCAMID=0 or 1
1000 int skiplimb = 0;
1001 //int skiplimb = 1; //!!!TEMP for Thanksgiving comet
1002 hcftid = drms_getkey_int(rs, "HCFTID", &rstatus);
1003 if(hcftid == CAL_HCFTID) { //don't call limb_fit()
1004 printk("Cal mode image fsn=%u\n", fsnx);
1005 skiplimb = 1;
1006 }
1007 else {
1008 if(hshiexp == 0) {
1009 if(hcamid == 0 || hcamid == 1) {
1010 printk("Dark image fsn=%u\n", fsnx);
1011 skiplimb = 1;
1012 }
1013 }
1014 }
1015 lstatus = 1; //default bad limb_fit
1016 //WCS calculations for condition 1
1017 //(For WCS conditions see mail from Rock 10Aug2010 16:40)
1018 if(!skiplimb && !hmiaiaflg) { //only call for HMI
1019 //StartTimer(1); //!!TEMP
1020 //lstatus = limb_fit(rs,l0l1->dat1.adata1,&rsun_lf,&x0_lf,&y0_lf,4096,4096,0, dpath);
1021 lstatus = limb_fit(rs,l0l1->dat1.adata1,&rsun_lf,&x0_lf,&y0_lf,4096,4096,0);
1022 if(lstatus) {
1023 printk("ERROR: limb_fit() %d error for fsn=%u\n", lstatus, fsnx);
1024 //noimage[i] = 1;
1025 limbmiss[i] = 1;
1026 }
1027 else { //limb_fit() returns good values
1028 //ftmp = StopTimer(1);
1029 //printf( "\nTime sec for limb_fit() = %f\n\n", ftmp );
1030 //printf("Calling WCS condition 1\n"); //!!TEMP
1031 drms_setkey_float(rs, "RSUN_LF", (float)rsun_lf);
1032 drms_setkey_float(rs, "R_SUN", (float)rsun_lf);
1033 drms_setkey_float(rs, "X0_LF", (float)x0_lf);
1034 drms_setkey_float(rs, "Y0_LF", (float)y0_lf);
1035 drms_setkey_float(rs, "CRVAL1", 0.0);
1036 drms_setkey_float(rs, "CRVAL2", 0.0);
1037 drms_setkey_string(rs, "CTYPE1", "HPLN-TAN");
1038 drms_setkey_string(rs, "CTYPE2", "HPLT-TAN");
1039 drms_setkey_string(rs, "CUNIT1", "arcsec");
1040 drms_setkey_string(rs, "CUNIT2", "arcsec");
1041 cdelt1 = (float)IOdata.rsun_obs/rsun_lf;
1042 drms_setkey_float(rs, "CDELT1", cdelt1);
1043 drms_setkey_float(rs, "CDELT2", cdelt1);
1044 crpix1 = (float)x0_lf + 1; //set up for correction by heightformation()
1045 crpix2 = (float)y0_lf + 1;
1046 drms_setkey_float(rs, "CRPIX1", crpix1);
1047 drms_setkey_float(rs, "CRPIX2", crpix2);
1048 rsun = (float)rsun_lf;
1049 crota2 = imageloc[i].instrot + ptdata.sat_rot;
1050 drms_setkey_float(rs, "CROTA2", crota2);
1051 goto WCSEND;
1052 }
1053 }
1054 //WCS calculations for condition 2
1055 //(As a first pass, this is following Rock's notes and can be simplified)
1056 int cond2 = 0;
1057 if(!strcmp(ptdata.acs_mode, "SCIENCE") || !strcmp(ptdata.acs_mode, DRMS_MISSING_STRING)) {
1058 //if(quicklook) {
1059 if(!skiplimb && !hmiaiaflg) { //hmi, not dark or cal
1060 cond2 = 1;
1061 }
1062 if(hmiaiaflg) {
1063 cond2 = 1;
1064 if(hshiexp == 0) {
1065 if(hcamid == 0 || hcamid == 1) { //it's aia dark
1066 cond2 = 0;
1067 }
1068 }
1069 }
1070 if(cond2) { //this is WCS condition 2
1071 /**********************************
1072 if(!hmiaiaflg) {
1073 cond2 = 0;
1074 if((rsun_lf == DRMS_MISSING_DOUBLE) || (x0_lf == DRMS_MISSING_DOUBLE) || (y0_lf == DRMS_MISSING_DOUBLE)) {
1075 cond2 = 1;
1076 }
1077 }
1078 **************************************/
1079 if(cond2) {
1080 printf("Calling WCS condition 2\n"); //!!TEMP
1081 rsun = (float)IOdata.rsun_obs/imageloc[i].imscale;
1082 drms_setkey_float(rs, "R_SUN", rsun);
1083 drms_setkey_string(rs, "CTYPE1", "HPLN-TAN");
1084 drms_setkey_string(rs, "CTYPE2", "HPLT-TAN");
1085 drms_setkey_float(rs, "CRVAL1", 0.0);
1086 drms_setkey_float(rs, "CRVAL2", 0.0);
1087 drms_setkey_string(rs, "CUNIT1", "arcsec");
1088 drms_setkey_string(rs, "CUNIT2", "arcsec");
1089 cdelt1 = (float)imageloc[i].imscale;
1090 drms_setkey_float(rs, "CDELT1", cdelt1);
1091 drms_setkey_float(rs, "CDELT2", cdelt1);
1092 crpix1 = imageloc[i].x + 1;
1093 crpix2 = imageloc[i].y + 1;
1094 drms_setkey_float(rs, "CRPIX1", crpix1);
1095 drms_setkey_float(rs, "CRPIX2", crpix2);
1096 crota2 = imageloc[i].instrot + ptdata.sat_rot;
1097 drms_setkey_float(rs, "CROTA2", crota2);
1098 goto WCSEND;
1099 }
1100 }
1101 //}
1102 }
1103 //WCS calculations for condition 3
1104 int cond3 = 0;
1105 if(strcmp(ptdata.acs_mode, "SCIENCE")) { //not is sci pointing mode
1106 if(!hmiaiaflg && !skiplimb && lstatus) { //hmi,limb fit NG, not dark or cal
1107 cond3 = 1;
1108 }
1109 else {
1110 if(hmiaiaflg) {
1111 cond3 = 1;
1112 if(hshiexp == 0) {
1113 if(hcamid == 0 || hcamid == 1) { //it's aia dark
1114 cond3 = 0;
1115 }
1116 }
1117 }
1118 }
1119 }
1120 /******extraneous already lstatus means limb fit ng
1121 if(cond3) {
1122 if(!hmiaiaflg) {
1123 cond3 = 0;
1124 if((rsun_lf == DRMS_MISSING_DOUBLE) || (x0_lf == DRMS_MISSING_DOUBLE) || (y0_lf == DRMS_MISSING_DOUBLE)) {
1125 cond3 = 1;
1126 }
1127 }
1128 }
1129 **************************************/
1130 if(cond3) {
1131 printf("Calling WCS condition 3\n"); //!!TEMP
1132 rsun = (float)IOdata.rsun_obs/imageloc[i].imscale;
1133 drms_setkey_float(rs, "R_SUN", rsun);
1134 drms_setkey_string(rs, "CTYPE1", "HPLN-TAN");
1135 drms_setkey_string(rs, "CTYPE2", "HPLT-TAN");
1136 drms_setkey_string(rs, "CUNIT1", "arcsec");
1137 drms_setkey_string(rs, "CUNIT2", "arcsec");
1138 drms_setkey_float(rs, "CRVAL1", 0.0);
1139 drms_setkey_float(rs, "CRVAL2", 0.0);
1140 cdelt1 = (float)imageloc[i].imscale;
1141 drms_setkey_float(rs, "CDELT1", cdelt1);
1142 drms_setkey_float(rs, "CDELT2", cdelt1);
1143 //below missing - SC_Y/Z_INRT_BIAS/IMSCL_MP + 1
1144 if(hmiaiaflg) { //aia
1145 crpix1 = imageloc[i].x + (ptdata.sat_y0 - imageloc[i].yinrtb)/imageloc[i].imscale + 1;
1146 crpix2 = imageloc[i].y + (ptdata.sat_z0 - imageloc[i].zinrtb)/imageloc[i].imscale + 1;
1147 }
1148 else { //hmi
1149 crpix1 = imageloc[i].x - (ptdata.sat_y0 - imageloc[i].yinrtb)/imageloc[i].imscale + 1;
1150 crpix2 = imageloc[i].y - (ptdata.sat_z0 - imageloc[i].zinrtb)/imageloc[i].imscale + 1;
1151 }
1152 drms_setkey_float(rs, "CRPIX1", crpix1);
1153 drms_setkey_float(rs, "CRPIX2", crpix2);
1154 crota2 = imageloc[i].instrot + ptdata.sat_rot;
1155 drms_setkey_float(rs, "CROTA2", crota2);
1156 goto WCSEND;
1157 }
1158 //WCS calculations for condition 4
1159 int cond4 = 0;
1160 if(!hmiaiaflg && skiplimb) { //hmi, dark or cal
1161 cond4 = 1;
1162 }
1163 else {
1164 if(hshiexp == 0) {
1165 if(hcamid == 0 || hcamid == 1) { //it's aia dark
1166 cond4 = 1;
1167 }
1168 }
1169 }
1170 if(cond4) {
1171 printf("Calling WCS condition 4\n"); //!!TEMP
1172 rsun = DRMS_MISSING_FLOAT;
1173 drms_setkey_float(rs, "R_SUN", rsun);
1174 drms_setkey_string(rs, "CTYPE1", "RAW");
1175 drms_setkey_string(rs, "CTYPE2", "RAW");
1176 drms_setkey_string(rs, "CUNIT1", DRMS_MISSING_STRING);
1177 drms_setkey_string(rs, "CUNIT2", DRMS_MISSING_STRING);
1178 drms_setkey_float(rs, "CRVAL1", DRMS_MISSING_FLOAT);
1179 drms_setkey_float(rs, "CRVAL2", DRMS_MISSING_FLOAT);
1180 cdelt1 = DRMS_MISSING_FLOAT;
1181 drms_setkey_float(rs, "CDELT1", cdelt1);
1182 drms_setkey_float(rs, "CDELT2", cdelt1);
1183 crpix1 = DRMS_MISSING_FLOAT;
1184 crpix2 = DRMS_MISSING_FLOAT;
1185 drms_setkey_float(rs, "CRPIX1", crpix1);
1186 drms_setkey_float(rs, "CRPIX2", crpix2);
1187 crota2 = imageloc[i].instrot + ptdata.sat_rot;
1188 drms_setkey_float(rs, "CROTA2", crota2);
1189 goto WCSEND;
1190 }
1191
1192 WCSEND:
1193 if(!hmiaiaflg && !lstatus) { //only do for hmi and good limb fit
1194 //Now call Sebastien's heightformation() fuction (email 08/09/10 17:50)
1195 //21Sep2010 change -crota2 to crota2 and HFCORRVR to 2
1196 if(!(dstatus = heightformation(fid, IOdata.obs_vr, &cdelt1, &rsun, &crpix1, &crpix2, crota2))) {
1197 drms_setkey_float(rs, "CDELT1", cdelt1);
1198 drms_setkey_float(rs, "CDELT2", cdelt1);
1199 drms_setkey_float(rs, "R_SUN", rsun);
1200 drms_setkey_float(rs, "CRPIX1", crpix1);
1201 drms_setkey_float(rs, "CRPIX2", crpix2);
1202 //drms_setkey_int(rs, "HFCORRVR", 2);
1203 }
1204 else {
1205 //drms_setkey_int(rs, "HFCORRVR", 0);
1206 printk("ERROR: heightformation() returned error for FID=%d\n", fid);
1207 }
1208 }
1209 drms_setkey_int(rs, "CALVER32", 0x12); //Put back in on 09Sep2013
1210
1211 //if(hmiaiaflg) { //aia
1212 // int wl = drms_getkey_int(rs, "WAVELNTH", &rstatus);
1213 //}
1214
1215 do_quallev1(rs0, rs, i, fsnx);
1216 //ftmp = StopTimer(2); //!!!TEMP
1217 //printf( "\nTime sec in inside loop for fsn=%u : %f\n\n", fsnx, ftmp );
1218
1219 } //END do for all the sorted lev0 records
1220
1221 drms_close_records(rset0, DRMS_FREE_RECORD); //close lev0 records
1222 drms_close_records(rset1, DRMS_INSERT_RECORD); //close lev1 records
1223 return(0);
1224 }
1225
1226 int compare_rptr(const void *a, const void *b)
1227 {
1228 TIME t1, t2;
1229 int rstatus;
1230 DRMS_Record_t *x=(DRMS_Record_t *)a, *y=(DRMS_Record_t *)b;
1231
1232 t1 = drms_getkey_time(x, "t_obs", &rstatus);
1233 if(rstatus) t1 = DRMS_MISSING_TIME; //treat error as missing t_obs
1234 t2 = drms_getkey_time(y, "t_obs", &rstatus);
1235 if(rstatus) t2 = DRMS_MISSING_TIME;
1236 if(t1 < t2) return(-1);
1237 if(t1 > t2) return(1);
1238 return(0);
1239 }
1240
1241 // Initial setup stuff called when main is first entered.
1242 void setup()
1243 {
1244 FILE *fin;
1245 char string[128], cwdbuf[128], idstr[256], lfile[128];
1246 int tpid;
1247
1248 sdo_epoch = sscan_time("1958.01.01_00:00:00_TAI");
1249 do_datestr();
1250 printk_set(h1log, h1log); // set for printk calls
1251 printk("%s\n", datestr);
1252 gethostname(idstr, 256);
1253 printf("Host: %s\n", idstr);
1254 printk("Host: %s\n", idstr);
1255 getcwd(cwdbuf, 126);
1256 sprintf(idstr, "Cwd: %s\nCall: ", cwdbuf);
1257 sprintf(string, "build_lev1_hmi started as pid=%d ppid=%d user=%s\n",
1258 getpid(), getppid(), username);
1259 strcat(idstr, string);
1260 printk("%s", idstr);
1261 printf("%s", idstr);
1262 if(restartflg) printk("-r ");
1263 sprintf(argmode, "mode=%s", mode);
1264 sprintf(arginstru, "instru=%s", instru);
1265 sprintf(argdsin, "dsin=%s", dsin);
1266 sprintf(argdsout, "dsout=%s", dsout);
1267 if(modeflg) {
1268 sprintf(argbx, "brec=%lld", brec);
1269 sprintf(argex, "erec=%lld", erec);
1270 }
1271 else {
1272 sprintf(argbx, "bfsn=%lld", bfsn);
1273 sprintf(argex, "efsn=%lld", efsn);
1274 }
1275 sprintf(argquick, "quicklook=%d", quicklook);
1276 sprintf(arglogfile, "logfile=%s", logname);
1277 printk("%s %s %s %s %s %s %s %s %s\n",
1278 argmode, arginstru, argdsin, argdsout, argbx, argex, argquick, arglogfile, argdsaiabad);
1279 printf("%s %s %s %s %s %s %s %s %s\n",
1280 argmode, arginstru, argdsin, argdsout, argbx, argex, argquick, arglogfile, argdsaiabad);
1281 if(!restartflg) {
1282 //printk("tlmseriesname=%s\nlev0seriesname=%s\n",
1283 // tlmseriesname, lev0seriesname);
1284 }
1285 sprintf(bld_vers, "%s", jsoc_version);
1286 sprintf(idstr, "ps -ef | grep %s", LEV1VIEWERNAME);
1287 fin = popen(idstr, "r");
1288 while(fgets(string, sizeof string, fin)) { //get ps line
1289 if(!(strstr(string, "perl"))) continue;
1290 sscanf(string, "%s %d", idstr, &tpid); /* get user name & process id */
1291 sprintf(lfile, "%s/build_lev1_hmi_restart_%d.touch", LEV1LOG_BASEDIR, tpid);
1292 sprintf(idstr, "/bin/touch %s", lfile);
1293 printk("%s\n", idstr);
1294 system(idstr);
1295 }
1296 umask(002); // allow group write
1297 //Image.initialized = 0; // init the two image structures
1298 //ImageOld.initialized = 0;
1299 //Img = &Image;
1300 //ImgO = &ImageOld;
1301 //Img->initialized = 0;
1302 //ImgO->initialized = 0;
1303 }
1304
1305 // Module main function.
1306 int DoIt(void)
1307 {
1308 long long numofrecs, frec, lrec;
1309 int numrec, numofchunks, i;
1310 char line[80];
1311 const char *dpath = NULL;
1312
1313 if (nice_intro())
1314 return (0);
1315 if(!(username = (char *)getenv("USER"))) username = "nouser";
1316 instru = cmdparams_get_str(&cmdparams, "instru", NULL);
1317 if(strcmp(instru, "hmi") && strcmp(instru, "aia")) {
1318 printf("Error: instru=%s must be given as 'hmi' or 'aia'\n", instru);
1319 printk("Error: instru=%s must be given as 'hmi' or 'aia'\n", instru);
1320 return(0);
1321 }
1322 mode = cmdparams_get_str(&cmdparams, "mode", NULL);
1323 if(strcmp(mode, "recnum") && strcmp(mode, "fsn")) {
1324 printf("Error: mode= must be given as 'recnum' or 'fsn'\n");
1325 return(0);
1326 }
1327 if(!strcmp(instru, "aia")) hmiaiaflg = 1;
1328 if(!strcmp(mode, "recnum")) modeflg = 1;
1329 dsin = cmdparams_get_str(&cmdparams, "dsin", NULL);
1330 dsout = cmdparams_get_str(&cmdparams, "dsout", NULL);
1331 dsff = cmdparams_get_str(&cmdparams, "dsff", NULL);
1332 dsaiabad = cmdparams_get_str(&cmdparams, "dsaiabad", NULL);
1333 if (strcmp(dsaiabad, NOTSPECIFIED)) sprintf(argdsaiabad, "dsaiabad=%s", dsaiabad);
1334 else sprintf(argdsaiabad, " ");
1335 brec = cmdparams_get_int(&cmdparams, "brec", NULL);
1336 erec = cmdparams_get_int(&cmdparams, "erec", NULL);
1337 bfsn = cmdparams_get_int(&cmdparams, "bfsn", NULL);
1338 efsn = cmdparams_get_int(&cmdparams, "efsn", NULL);
1339 quicklook = cmdparams_get_int(&cmdparams, "quicklook", NULL);
1340
1341 dpath = cmdparams_get_str(&cmdparams, kDpath, NULL);
1342
1343 //quicklook = 1; //!!TEMP for test
1344 if(modeflg) { //recnum mode
1345 if(brec == 0 || erec == 0) {
1346 fprintf(stderr, "brec and erec must be given for recnum mode. 0 not allowed\n");
1347 return(0);
1348 }
1349 if(brec > erec) {
1350 fprintf(stderr, "brec must be <= erec\n");
1351 return(0);
1352 }
1353 bnumx = brec;
1354 enumx = erec;
1355 }
1356 else { //fsn mode
1357 if(bfsn == 0 || efsn == 0) {
1358 fprintf(stderr, "bfsn and efsn must be given for fsn mode. 0 not allowed\n");
1359 return(0);
1360 }
1361 if(bfsn > efsn) {
1362 fprintf(stderr, "bfsn must be <= efsn\n");
1363 return(0);
1364 }
1365 bnumx = bfsn;
1366 enumx = efsn;
1367 }
1368 logfile = cmdparams_get_str(&cmdparams, "logfile", NULL);
1369 if (strcmp(dsin, NOTSPECIFIED) == 0) {
1370 if(hmiaiaflg == 0) dsin = LEV0SERIESNAMEHMI;
1371 else dsin = LEV0SERIESNAMEAIA;
1372 }
1373 if (strcmp(dsout, NOTSPECIFIED) == 0) {
1374 if(hmiaiaflg == 0) dsout = LEV1SERIESNAMEHMI;
1375 else dsout = LEV1SERIESNAMEAIA;
1376 }
1377 if(hmiaiaflg) {
1378 if (strcmp(dsff, NOTSPECIFIED)) sprintf(dsffname, "%s", dsff);
1379 else sprintf(dsffname, "%s", DSFFNAMEAIA);
1380 if(strstr(dsin, "hmi") || strstr(dsout, "hmi")) {
1381 printf("Warning: You said instru=aia but have 'hmi' in ds name?\n");
1382 printf("Do you want to abort this [y/n]? ");
1383 if(gets(line) == NULL) { return(0); }
1384 if(strcmp(line, "n")) { return(0); }
1385 }
1386 }
1387 else {
1388 sprintf(dsffname, "%s", DSFFNAMEHMI);
1389 if(strstr(dsin, "aia") || strstr(dsout, "aia")) {
1390 printf("Warning: You said instru=hmi but have 'aia' in ds name?\n");
1391 printf("Do you want to abort this [y/n]? ");
1392 if(gets(line) == NULL) { return(0); }
1393 if(strcmp(line, "n")) { return(0); }
1394 }
1395 }
1396 if (strcmp(logfile, NOTSPECIFIED) == 0) {
1397 sprintf(logname, H1LOGFILE, gettimetag());
1398 }
1399 else {
1400 sprintf(logname, "%s", logfile);
1401 }
1402 if(restartflg || LOGTEST) {
1403 if((h1logfp=fopen(logname, "a")) == NULL)
1404 fprintf(stderr, "**Can't open for append the log file %s\n", logname);
1405 }
1406 else {
1407 if((h1logfp=fopen(logname, "w")) == NULL)
1408 fprintf(stderr, "**Can't open the log file %s\n", logname);
1409 }
1410 setup();
1411 numofrecs = (enumx - bnumx) + 1;
1412 numrec = NUMRECLEV1; //# of records to do at a time
1413 numofchunks = numofrecs/numrec;
1414 if((numofrecs % numrec) != 0) numofchunks++; //extra loop for partial chunk
1415 lrec = bnumx-1;
1416 for(i = 0; i < numofchunks; i++) {
1417 frec = lrec+1; lrec = (frec + numrec)-1;
1418 if(lrec > enumx) lrec=enumx;
1419 if(do_ingest(frec, lrec, dpath)) { //do a chunk to get files from the lev0
1420 printf("build_lev1_hmi abort\nSee log: %s\n", logname);
1421 send_mail("build_lev1_hmi abort\nSee log: %s\n", logname);
1422 return(0);
1423 }
1424 if(modeflg) { //only do for recnum mode
1425 if(stat(STOP_FILE, &stbuf) == 0) {
1426 printf("Stop file %s seen. Exit build_lev1_hmi.\n", STOP_FILE);
1427 break;
1428 }
1429 }
1430 }
1431 printf("build_lev1_hmi done last fsn=%u\n", fsnx);
1432 return(0);
1433 }