1 mbobra 1.1 /*
2 * smarp.c
3 *
4 * This module creates the pipeline for Space Weather MDI Active Region Patches (SMARPs).
5 * It is a modified version of sharp.c, created by Xudong Sun and Monica Bobra.
6 * It takes the mdi.mtarp series to create the following:
7 *
8 * Series 1: mdi.smarp_cea_96m
9 * CEA remapped magnetogram, bitmap, continuum (same size in map coordinate)
10 * Space weather indices based on line-of-sight magnetogram in the cutout series
11 *
12 * Series 2: mdi.smarp_96m
13 * cutouts of magnetogram, bitmap, continuum, (TARP defined, various sizes in CCD pixels)
14 * Space weather indices based on line-of-sight magnetogram in the cutout series
15 *
16 * Author:
17 * Monica Bobra; Xudong Sun
18 *
19 * Version:
20 * v0.0 9 February 2018 Monica Bobra
21 *
22 mbobra 1.1 * Notes:
23 * v0.0 No explicit notes
24 *
25 * Example Call:
26 * > smarp "mharp=mdi.Mtarp[13643][2010.10.14_20:48:00_TAI]" "sharp_cea=su_mbobra.smarp_cea_96m" cont="mdi.fd_Ic_interp[2010.10.14_20:48:00_TAI]" "sharp_cut=su_mbobra.smarp_96m"
27 *
28 */
29
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <time.h>
33 #include <sys/time.h>
34 #include <math.h>
35 #include <string.h>
36 #include "jsoc_main.h"
37 #include "astro.h"
38 #include "fstats.h"
39 #include "cartography.c"
40 #include "fresize.h"
41 #include "finterpolate.h"
42 #include "img2helioVector.c"
43 mbobra 1.1 #include "copy_me_keys.c"
44 #include "errorprop.c"
45 #include "smarp_functions.c"
46
47 //#include <mkl.h> // Comment out mkl.h, which can only run on solar3
48 #include <mkl_blas.h>
49 #include <mkl_service.h>
50 #include <mkl_lapack.h>
51 #include <mkl_vml_functions.h>
52 #include <omp.h>
53
54 #define PI (M_PI)
55 #define RADSINDEG (PI/180.)
56 #define RAD2ARCSEC (648000./M_PI)
57 #define SECINDAY (86400.)
58 #define FOURK (1024)
59 #define FOURK2 (1048576)
60
61 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
62
63 // FOR HMI: Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
64 mbobra 1.1 // FOR HMI: Nyqvist rate at disk center is 0.12 degree. Oversample above 0.06 degree
65 #define NYQVIST (0.06)
66
67 // Maximum variation of LONDTMAX-LONDTMIN
68 #define MAXLONDIFF (1.2e-4)
69
70 // Some other things
71 #ifndef MIN
72 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
73 #endif
74 #ifndef MAX
75 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
76 #endif
77
78 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
79 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
80
81 #define kNotSpecified "Not Specified"
82
83 // Macros for WCS transformations. assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
84 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
85 mbobra 1.1 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
86 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
87 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
88 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
89 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
90 #define XSCALE 0.12
91 #define YSCALE 0.12
92 #define NBIN 3
93 #define INTERP 0
94 #define dpath "/home/jsoc/cvs/Development/JSOC"
95
96 /* ========================================================================================================== */
97
98 // Space weather keywords
99 struct swIndex {
100 float mean_vf;
101 float count_mask;
102 float absFlux;
103 float mean_derivative_bz;
104 float Rparam;
105 };
106 mbobra 1.1
107 // Mapping method
108 enum projection {
109 carree,
110 cassini,
111 mercator,
112 cyleqa,
113 sineqa,
114 gnomonic,
115 postel,
116 stereographic,
117 orthographic,
118 lambert
119 };
120
121 // WSC code
122 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
123 "SIN", "ZEA"};
124
125 // Ephemeris information
126 struct ephemeris {
127 mbobra 1.1 double disk_lonc, disk_latc;
128 double disk_xc, disk_yc;
129 double rSun, asd, pa;
130 };
131
132 // Mapping information
133 struct mapInfo {
134 float xc, yc; // reference point: center
135 int nrow, ncol; // size
136 float xscale, yscale; // scale
137 int nbin;
138 enum projection proj; // projection method
139 struct ephemeris ephem; // ephemeris info
140 float *xi_out, *zeta_out; // coordinate on full disk image to sample at
141 };
142
143 /* ========================================================================================================== */
144
145 /* Get all input data series */
146 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery);
147
148 mbobra 1.1 /* Get other data series */
149 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
150
151 /* Find record from record set with given T_rec */
152 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
153
154 /* Create CEA record */
155 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
156
157 /* Mapping single segment, wrapper */
158 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
159 struct mapInfo *mInfo, char *segName);
160
161 /* Determine reference point coordinate and patch size according to input */
162 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
163
164 /* Get ephemeris information */
165 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
166
167 /* Compute the coordinates at which the full disk image is sampled */
168 void findCoord(struct mapInfo *mInfo);
169 mbobra 1.1
170 /* Mapping function */
171 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt);
172
173 // ===================
174
175 /* Create Cutout record */
176 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
177
178 /* Get cutout and write segment */
179 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
180
181 // ===================
182
183 /* Compute space weather indices */
184 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
185
186 /* Set space weather indices */
187 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
188
189 /* Set all keywords */
190 mbobra 1.1 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo);
191
192 // ===================
193
194 /* Nearest neighbor interpolation */
195 float nnb (float *f, int nx, int ny, double x, double y);
196
197 /* Wrapper for Jesper's rebin code */
198 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
199
200 /* ========================================================================================================== */
201
202 /* Cutout segment names, input identical to output */
203 char *MharpSegs[] = {"magnetogram", "bitmap"};
204 char *CutSegs[] = {"magnetogram", "bitmap", "continuum"};
205 char *CEASegs[] = {"magnetogram", "bitmap", "continuum"};
206 // For BUNIT
207 char *CutBunits[] = {"Mx/cm^2", " ", "DN/s"};
208 char *CEABunits[] = {"Mx/cm^2", " ", "DN/s"};
209 /* ========================================================================================================== */
210
211 mbobra 1.1 char *module_name = "smarp";
212
213 ModuleArgs_t module_args[] =
214 {
215 {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
216 {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
217 {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
218 {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
219 {ARG_END}
220 };
221
222 int DoIt(void)
223 {
224 int errbufstat = setvbuf(stderr, NULL, _IONBF, BUFSIZ);
225 int outbufstat = setvbuf(stdout, NULL, _IONBF, BUFSIZ);
226 int status = DRMS_SUCCESS;
227 int nrecs, irec;
228 char *mharpQuery;
229 char *contQuery;
230 char *sharpCeaQuery, *sharpCutQuery;
231 DRMS_RecordSet_t *mharpRS = NULL;
232 mbobra 1.1 DRMS_RecordSet_t *contRS = NULL;
233
234 /* Get parameters */
235
236 mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
237 sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
238 sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
239 contQuery = (char *) params_get_str(&cmdparams, "cont");
240
241 /* Get input data, check everything */
242 if (getInputRS(&mharpRS, mharpQuery))
243 DIE("Input harp data error.");
244 nrecs = mharpRS->n;
245
246 if (getInputRS_aux(&contRS, contQuery, mharpRS))
247 DIE("Input continuum data error.");
248
249 /* Start */
250
251 printf("==============\nStart. %d image(s) in total.\n", nrecs);
252
253 mbobra 1.1 for (irec = 0; irec < nrecs; irec++) {
254
255 /* Records in work */
256
257 DRMS_Record_t *mharpRec = NULL;
258 DRMS_Record_t *contRec = NULL;
259
260 mharpRec = mharpRS->records[irec];
261 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
262
263 struct swIndex swKeys;
264
265 if (getInputRec_aux(&contRec, contRS, trec)) {
266 printf("Fetching Continuum failed, image #%d skipped.\n", irec);
267 continue;
268 }
269
270 printf("Obtained all the data \n");
271
272 /* Create CEA record */
273
274 mbobra 1.1 DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
275 if (status) { // if failed
276 printf("Creating CEA failed, image #%d skipped.\n", irec);
277 continue;
278 }
279 if (createCeaRecord(mharpRec, contRec, sharpCeaRec, &swKeys)) { // do the work
280 printf("Creating CEA failed, image #%d skipped.\n", irec);
281 drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
282 continue;
283 } // swKeys updated here
284
285 drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
286
287 printf("Created CEA record \n");
288
289 /* Create Cutout record */
290
291 DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
292 if (status) { // if failed
293 printf("Creating cutout failed, image #%d skipped.\n", irec);
294 continue;
295 mbobra 1.1 }
296
297 if (createCutRecord(mharpRec, contRec, sharpCutRec, &swKeys)) { // do the work
298 printf("Creating cutout failed, image #%d skipped.\n", irec);
299 drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
300 continue;
301 } // swKeys used here
302 drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
303 printf("Created CUT record \n");
304 /* Done */
305
306 printf("Image #%d done.\n", irec);
307
308 } // irec
309
310
311 drms_close_records(mharpRS, DRMS_FREE_RECORD);
312 drms_close_records(contRS, DRMS_FREE_RECORD);
313
314 return 0;
315
316 mbobra 1.1 } // DoIt
317
318 // ===================================================================
319 // ===================================================================
320 // ===================================================================
321
322 /*
323 * Get input data series, including mHarp and bharp
324 * Need all records to match, otherwise quit
325 *
326 */
327
328 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, char *mharpQuery)
329 {
330 int status = 0;
331 *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
332 if (status || (*mharpRS_ptr)->n == 0) return 1;
333 return 0;
334 }
335
336 /*
337 mbobra 1.1 * Get other data series, check all T_REC are available
338 */
339
340 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
341 {
342
343 int status = 0;
344
345 *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
346 if (status || (*inRS_ptr)->n == 0) return status;
347
348 // Check if all T_rec are available, need to match both ways
349 int n = harpRS->n, n0 = (*inRS_ptr)->n;
350
351 for (int i0 = 0; i0 < n0; i0++) {
352 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
353 TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
354 TIME trec = 0;
355 for (int i = 0; i < n; i++) {
356 DRMS_Record_t *harpRec = harpRS->records[i];
357 trec = drms_getkey_time(harpRec, "T_REC", &status);
358 mbobra 1.1 if (fabs(trec0 - trec) < 10) break;
359 }
360 if (fabs(trec0 - trec) >= 10) return 1;
361 }
362
363 for (int i = 0; i < n; i++) {
364 DRMS_Record_t *harpRec = harpRS->records[i];
365 TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
366 TIME trec0 = 0;
367 for (int i0 = 0; i0 < n0; i0++) {
368 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
369 trec0 = drms_getkey_time(inRec, "T_REC", &status);
370 if (fabs(trec0 - trec) < 10) break;
371 }
372 if (fabs(trec0 - trec) >= 10) return 1;
373 }
374
375 return 0;
376
377 }
378
379 mbobra 1.1 /*
380 * Find record from record set with given T_rec
381 */
382
383 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
384 {
385
386 int status = 0;
387
388 int n = inRS->n;
389 for (int i = 0; i < n; i++) {
390 *inRec_ptr = inRS->records[i];
391 TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
392 if (fabs(trec0 - trec) < 10) return 0;
393 }
394
395 return 1;
396
397 }
398
399
400 mbobra 1.1
401
402 /*
403 * Create CEA record: top level subroutine
404 * Also compute all the space weather keywords here
405 */
406
407 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
408 {
409
410 int status = 0;
411 DRMS_Segment_t *inSeg;
412 DRMS_Array_t *inArray;
413 int val;
414
415 struct mapInfo mInfo;
416 mInfo.proj = (enum projection) cyleqa; // projection method
417 mInfo.xscale = XSCALE;
418 mInfo.yscale = YSCALE;
419
420 int ncol0, nrow0; // oversampled map size
421 mbobra 1.1
422 // Get ephemeris
423
424 if (getEphemeris(mharpRec, &(mInfo.ephem))) {
425 SHOW("CEA: get ephemeris error\n");
426 return 1;
427 }
428
429 // Find position
430
431 if (findPosition(mharpRec, &mInfo)) {
432 SHOW("CEA: find position error\n");
433 return 1;
434 }
435
436 // ========================================
437 // Do this for all bitmaps, Aug 12 2013 XS
438 // ========================================
439
440 mInfo.nbin = 1; // for bitmaps. suppress anti-aliasing
441 ncol0 = mInfo.ncol;
442 mbobra 1.1 nrow0 = mInfo.nrow;
443
444 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
445 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
446
447 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
448
449 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
450 SHOW("CEA: mapping bitmap error\n");
451 return 1;
452 }
453 printf("Bitmap mapping done.\n");
454
455 free(mInfo.xi_out);
456 free(mInfo.zeta_out);
457
458 // ========================================
459 // Do this again for floats, Aug 12 2013 XS
460 // ========================================
461 // Create xi_out, zeta_out array in mInfo:
462 // Coordinates to sample in original full disk image
463 mbobra 1.1
464 mInfo.nbin = NBIN;
465 ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
466 nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
467
468 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
469 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
470
471 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
472
473 // Mapping single segment: Mharp, etc.
474
475 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
476 SHOW("CEA: mapping magnetogram error\n");
477 return 1;
478 }
479 printf("Magnetogram mapping done.\n");
480
481 if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
482 SHOW("CEA: mapping continuum error\n");
483 return 1;
484 mbobra 1.1 }
485 printf("Intensitygram mapping done.\n");
486
487 // Keywords & Links
488 copy_patch_keys(mharpRec, sharpRec);
489 copy_geo_keys(mharpRec, sharpRec);
490 // rename HARPNUM to TARPNUM
491 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
492 drms_setkey_double(sharpRec, "TARPNUM", val);
493 // copy everything else
494 drms_copykey(sharpRec, mharpRec, "T_REC");
495 drms_copykey(sharpRec, mharpRec, "CDELT1");
496 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
497 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
498 drms_copykey(sharpRec, mharpRec, "OBS_VR");
499 drms_copykey(sharpRec, mharpRec, "OBS_VW");
500 drms_copykey(sharpRec, mharpRec, "OBS_VN");
501 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
502 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
503 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
504 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
505 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
506 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
507 drms_copykey(sharpRec, mharpRec, "T_OBS");
508 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
509 drms_copykey(sharpRec, mharpRec, "QUALITY");
510 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
511 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
512 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
513 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
514 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
515 drms_copykey(sharpRec, mharpRec, "T_MID1");
516 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
517
518 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
519 if (mHarpLink) {
520 drms_link_set("MTARP", sharpRec, mharpRec);
521 }
522
523 // set other keywords
524 setKeys(sharpRec, mharpRec, &mInfo);
525
526 mbobra 1.1 // set space weather keywords
527 computeSWIndex(swKeys_ptr, sharpRec, &mInfo);
528 printf("Space weather indices done.\n");
529 setSWIndex(sharpRec, swKeys_ptr);
530
531 // set statistical keywords (e.g. DATAMIN, DATAMAX, etc.)
532 //int nCEASegs = ARRLENGTH(CEASegs);
533 int nCEASegs = 3;
534 for (int iSeg = 0; iSeg < 3; iSeg++) {
535 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
536 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
537 int stat = set_statistics(outSeg, outArray, 1);
538 //printf("%d => %d\n", iSeg, stat);
539 drms_free_array(outArray);
540 }
541
542 free(mInfo.xi_out);
543 free(mInfo.zeta_out);
544 return 0;
545
546 }
547 mbobra 1.1
548 /*
549 * Mapping a single segment
550 * Read in full disk image, utilize mapImage for mapping
551 * then write the segment out, segName same in in/out Rec
552 */
553
554 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
555 struct mapInfo *mInfo, char *segName)
556 {
557
558 int status = 0;
559 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
560 int dims[2] = {nx, ny};
561 int interpOpt = INTERP; // Aug 12 XS, default, overridden below for bitmaps and conf_disambig
562
563 // Input full disk array
564
565 DRMS_Segment_t *inSeg = NULL;
566 inSeg = drms_segment_lookup(inRec, segName);
567 if (!inSeg) return 1;
568 mbobra 1.1
569 DRMS_Array_t *inArray = NULL;
570
571 inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
572 if (!inArray) return 1;
573
574 float *inData;
575 int xsz = inArray->axis[0], ysz = inArray->axis[1];
576 if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk
577 float *inData0 = (float *) inArray->data;
578 inData = (float *) (calloc(FOURK2, sizeof(float)));
579 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
580 int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
581 int ind_map;
582 for (int row = 0; row < ysz; row++) {
583 for (int col = 0; col < xsz; col++) {
584 ind_map = (row + y0) * FOURK + (col + x0);
585 inData[ind_map] = inData0[row * xsz + col];
586 }
587 }
588 drms_free_array(inArray); inArray = NULL;
589 mbobra 1.1 } else {
590 inData = (float *) inArray->data;
591 }
592
593 // Mapping
594
595 float *map = (float *) (malloc(nxny * sizeof(float)));
596 if (performSampling(map, inData, mInfo, interpOpt)) // Add interpOpt for different types, Aug 12 XS
597 {if (inArray) drms_free_array(inArray); free(map); return 1;}
598
599 // Write out
600
601 DRMS_Segment_t *outSeg = NULL;
602 outSeg = drms_segment_lookup(sharpRec, segName);
603 if (!outSeg) return 1;
604
605 // DRMS_Type_t arrayType = outSeg->info->type;
606 DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
607 if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
608
609 // convert to needed data type
610 mbobra 1.1
611 // drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013
612
613 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
614 // outArray->parent_segment = outSeg;
615 outArray->israw = 0; // always compressed
616 outArray->bzero = outSeg->bzero;
617 outArray->bscale = outSeg->bscale;
618
619 status = drms_segment_write(outSeg, outArray, 0);
620 if (status) return 0;
621
622 if (inArray) drms_free_array(inArray);
623 if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012
624 if (outArray) drms_free_array(outArray);
625 return 0;
626
627 }
628
629 /*
630 * Determine reference point coordinate and patch size according to keywords
631 mbobra 1.1 * xc, yc are the coordinate of patch center, in degrees
632 * ncol and nrow are the final size
633 */
634
635 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
636 {
637
638 int status = 0;
639 int harpnum = drms_getkey_int(inRec, "TARPNUM", &status);
640 TIME trec = drms_getkey_time(inRec, "T_REC", &status);
641 float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
642
643 /* Center coord */
644 // Changed into double Jun 16 2014 XS
645
646 double minlon = drms_getkey_double(inRec, "LONDTMIN", &status); if (status) return 1; // Stonyhurst lon
647 double maxlon = drms_getkey_double(inRec, "LONDTMAX", &status); if (status) return 1;
648 double minlat = drms_getkey_double(inRec, "LATDTMIN", &status); if (status) return 1;
649 double maxlat = drms_getkey_double(inRec, "LATDTMAX", &status); if (status) return 1;
650
651 // A bug fixer for HARP (per M. Turmon)
652 mbobra 1.1 // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
653 // Also keywords such as "SIZE" will be NaN
654 // We compute minlon & minlat then by
655 // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
656
657 // float psize = drms_getkey_float(inRec, "SIZE", &status);
658 // if (psize != psize) {
659
660 if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE
661 TIME t0 = drms_getkey_time(inRec, "T_FRST1", &status); if (status) return 1; // changed from T_FRST to T_FRST1, T_FRST may not exist
662 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
663 char firstRecQuery[100], t0_str[100];
664 sprint_time(t0_str, t0, "TAI", 0);
665 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
666 DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
667 if (status || tmpRS->n != 1) return 1;
668 DRMS_Record_t *tmpRec = tmpRS->records[0];
669 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
670 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
671 minlon = minlon0 + (trec - t0) * omega / SECINDAY;
672 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
673 mbobra 1.1 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
674 }
675
676 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
677 mInfo->yc = (maxlat + minlat) / 2.;
678
679 /* Size */
680 // Rounded to 1.d3 precision first. Jun 16 2014 XS
681 // The previous fix does not work. LONDTMAX-LONDTMIN varies from frame to frame
682 // Need to find out the maximum possible difference, MAXLONDIFF (1.2e-4)
683 // Now, ncol = (maxlon-minlon)/xscale, if the decimal part is outside 0.5 \pm (MAXLONDIFF/xscale)
684 // proceed as it is. else, all use floor on ncol
685
686 float dpix = (MAXLONDIFF / mInfo->xscale) * 1.5; // "danger zone"
687 float ncol = (maxlon - minlon) / mInfo->xscale;
688 float d_ncol = fabs(ncol - floor(ncol) - 0.5); // distance to 0.5
689 if (d_ncol < dpix) {
690 mInfo->ncol = floor(ncol);
691 } else {
692 mInfo->ncol = round(ncol);
693 }
694 mbobra 1.1
695 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
696
697 return 0;
698
699 }
700
701
702 /*
703 * Fetch ephemeris info from a DRMS record
704 */
705
706 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
707 {
708
709 int status = 0;
710
711 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
712 double sina = sin(crota2 * RADSINDEG);
713 double cosa = cos(crota2 * RADSINDEG);
714
715 mbobra 1.1 ephem->pa = - crota2 * RADSINDEG;
716 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
717 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
718
719 float crvalx = 0.0;
720 float crvaly = 0.0;
721 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
722 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
723 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
724 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0
725 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
726
727 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
728 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
729 if (status) rSun_ref = 6.96e8;
730
731 ephem->asd = asin(rSun_ref/dSun);
732 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
733
734 return 0;
735
736 mbobra 1.1 }
737
738
739 /*
740 * Compute the coordinates to be sampled on full disk image
741 * mInfo->xi_out & mInfo->zeta_out
742 * This is oversampled, its size is ncol0 & nrow0 as shown below
743 */
744
745 void findCoord(struct mapInfo *mInfo)
746 {
747
748 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
749 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
750
751 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution
752 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad
753
754 double lonc = mInfo->xc * RADSINDEG; // in rad
755 double latc = mInfo->yc * RADSINDEG;
756
757 mbobra 1.1 double disk_lonc = (mInfo->ephem).disk_lonc;
758 double disk_latc = (mInfo->ephem).disk_latc;
759
760 double rSun = (mInfo->ephem).rSun;
761 double disk_xc = (mInfo->ephem).disk_xc / rSun;
762 double disk_yc = (mInfo->ephem).disk_yc / rSun;
763 double pa = (mInfo->ephem).pa;
764
765 // Temp pointers
766
767 float *xi_out = mInfo->xi_out;
768 float *zeta_out = mInfo->zeta_out;
769
770 // start
771
772 double x, y; // map coord
773 double lat, lon; // helio coord
774 double xi, zeta; // image coord (for one point)
775
776 int ind_map;
777
778 mbobra 1.1 for (int row0 = 0; row0 < nrow0; row0++) {
779 for (int col0 = 0; col0 < ncol0; col0++) {
780
781 ind_map = row0 * ncol0 + col0;
782
783 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad
784 y = (row0 + 0.5 - nrow0/2.) * yscale0;
785
786 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
787 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
788 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
789 * is the heliographic longitude and latitude of the map center. Both are in degree.
790 */
791
792 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
793 xi_out[ind_map] = -1;
794 zeta_out[ind_map] = -1;
795 continue;
796 }
797
798 /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
799 mbobra 1.1 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
800 */
801
802 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
803 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
804 xi_out[ind_map] = -1;
805 zeta_out[ind_map] = -1;
806 continue;
807 }
808
809 xi_out[ind_map] = xi * rSun;
810 zeta_out[ind_map] = zeta * rSun;
811
812 }
813 }
814
815 }
816
817
818 /*
819 * Sampling function
820 mbobra 1.1 * oversampling by nbin, then binning using a Gaussian
821 * save results in outData, always of float type
822 */
823
824 int performSampling(float *outData, float *inData, struct mapInfo *mInfo, int interpOpt)
825 {
826
827 int status = 0;
828 int ind_map;
829
830 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
831 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
832
833 // Changed Aug 12 2013, XS, for bitmaps
834 float *outData0;
835 if (interpOpt == 3 && mInfo->nbin == 1) {
836 outData0 = outData;
837 } else {
838 outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
839 }
840
841 mbobra 1.1 float *xi_out = mInfo->xi_out;
842 float *zeta_out = mInfo->zeta_out;
843
844 // Interpolation
845
846 struct fint_struct pars;
847 // Aug 12 2013, passed in as argument now
848
849 switch (interpOpt) {
850 case 0: // Wiener, 6 order, 1 constraint
851 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
852 break;
853 case 1: // Cubic convolution
854 init_finterpolate_cubic_conv(&pars, 1., 3.);
855 break;
856 case 2: // Bilinear
857 init_finterpolate_linear(&pars, 1.);
858 break;
859 case 3: // Near neighbor
860 break;
861 default:
862 mbobra 1.1 return 1;
863 }
864
865 //printf("interpOpt = %d, nbin = %d ", interpOpt, mInfo->nbin);
866 if (interpOpt == 3) { // Aug 6 2013, Xudong
867 for (int row0 = 0; row0 < nrow0; row0++) {
868 for (int col0 = 0; col0 < ncol0; col0++) {
869 ind_map = row0 * ncol0 + col0;
870 outData0[ind_map] = nnb(inData, FOURK, FOURK, xi_out[ind_map], zeta_out[ind_map]);
871 }
872 }
873 } else {
874 finterpolate(&pars, inData, xi_out, zeta_out, outData0,
875 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
876 }
877
878 // Rebinning, smoothing
879
880 if (interpOpt == 3 && mInfo->nbin == 1) {
881 return 0;
882 } else {
883 mbobra 1.1 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
884 free(outData0);
885 }
886
887 return 0;
888
889 }
890
891 /*
892 * Create Cutout record: top level subroutine
893 * Do the loops on segments and set the keywords here
894 * Work is done in writeCutout routine below
895 */
896
897 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *contRec, DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
898 {
899
900 int status = 0;
901 int val;
902 int iHarpSeg;
903 int nMharpSegs = ARRLENGTH(MharpSegs);
904 mbobra 1.1
905 // Cutout Mharp
906
907 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
908 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
909 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
910 printf("iHarpSeg nMharpSegs %d %d \n",iHarpSeg,nMharpSegs);
911 break;
912 }
913 }
914 if (iHarpSeg != nMharpSegs) {
915 SHOW("Cutout: segment number mismatch\n");
916 return 1; // if failed
917 }
918 printf("Magnetogram cutout done.\n");
919
920 // Cutout Continuum
921
922 if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
923 printf("Continuum cutout failed\n");
924 return 1;
925 mbobra 1.1 }
926 printf("Intensitygram cutout done.\n");
927
928 // Keywords & Links
929 copy_patch_keys(mharpRec, sharpRec);
930 copy_geo_keys(mharpRec, sharpRec);
931 // rename HARPNUM to TARPNUM
932 val = drms_getkey_double(mharpRec, "HARPNUM", &status);
933 drms_setkey_double(sharpRec, "TARPNUM", val);
934 // copy everything else
935 drms_copykey(sharpRec, mharpRec, "T_REC");
936 drms_copykey(sharpRec, mharpRec, "CDELT1");
937 drms_copykey(sharpRec, mharpRec, "RSUN_OBS");
938 drms_copykey(sharpRec, mharpRec, "DSUN_OBS");
939 drms_copykey(sharpRec, mharpRec, "OBS_VR");
940 drms_copykey(sharpRec, mharpRec, "OBS_VW");
941 drms_copykey(sharpRec, mharpRec, "OBS_VN");
942 drms_copykey(sharpRec, mharpRec, "CRLN_OBS");
943 drms_copykey(sharpRec, mharpRec, "CRLT_OBS");
944 drms_copykey(sharpRec, mharpRec, "CAR_ROT");
945 drms_copykey(sharpRec, mharpRec, "SIZE_SPT");
946 mbobra 1.1 drms_copykey(sharpRec, mharpRec, "AREA_SPT");
947 drms_copykey(sharpRec, mharpRec, "DATE__OBS");
948 drms_copykey(sharpRec, mharpRec, "T_OBS");
949 drms_copykey(sharpRec, mharpRec, "T_MAXPIX");
950 drms_copykey(sharpRec, mharpRec, "QUALITY");
951 drms_copykey(sharpRec, mharpRec, "NPIX_SPT");
952 drms_copykey(sharpRec, mharpRec, "ARS_NCLN");
953 drms_copykey(sharpRec, mharpRec, "ARS_MODL");
954 drms_copykey(sharpRec, mharpRec, "ARS_EDGE");
955 drms_copykey(sharpRec, mharpRec, "ARS_BETA");
956 drms_copykey(sharpRec, mharpRec, "T_MID1");
957 drms_copykey(sharpRec, mharpRec, "T_CMPASS");
958
959 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MTARP");
960 if (mHarpLink) {
961 drms_link_set("MTARP", sharpRec, mharpRec);
962 }
963
964 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
965 setKeys(sharpRec, mharpRec, NULL); // Set all other keywords, NULL specifies cutout
966
967 mbobra 1.1 // Stats
968
969 int nCutSegs = 3;
970 for (int iSeg = 0; iSeg < 3; iSeg++) {
971 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
972 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
973 set_statistics(outSeg, outArray, 1);
974 drms_free_array(outArray);
975 }
976
977 return 0;
978
979 }
980
981 /*
982 * Get cutout and write segment
983 */
984
985 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
986 {
987
988 mbobra 1.1 int status = 0;
989
990 DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
991 DRMS_Array_t *cutoutArray = NULL;
992 // DRMS_Type_t arrayType;
993
994 int ll[2], ur[2], nx, ny, nxny; // lower-left and upper right coords
995
996 /* Info */
997
998 inSeg = drms_segment_lookup(inRec, SegName);
999 if (!inSeg) return 1;
1000 //printf("SegName=%s\n",SegName); fflush(stdout);
1001 nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
1002 ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
1003 nxny = nx * ny;
1004 ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
1005 ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
1006 ur[0] = ll[0] + nx - 1; if (status) return 1;
1007 ur[1] = ll[1] + ny - 1; if (status) return 1;
1008 if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) { // for bitmaps, infomaps, etc.
1009 mbobra 1.1 cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
1010 if (status) return 1;
1011 } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) { // for full disk ones
1012 cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
1013 if (status) return 1;
1014 } else {
1015 return 1;
1016 }
1017 /* Write out */
1018 outSeg = drms_segment_lookup(outRec, SegName);
1019 if (!outSeg) return 1;
1020 outSeg->axis[0] = cutoutArray->axis[0];
1021 outSeg->axis[1] = cutoutArray->axis[1];
1022 cutoutArray->israw = 0; // always compressed
1023 cutoutArray->bzero = outSeg->bzero;
1024 cutoutArray->bscale = outSeg->bscale; // Same as inArray's
1025 status = drms_segment_write(outSeg, cutoutArray, 0);
1026 drms_free_array(cutoutArray);
1027 if (status) return 1;
1028 //printf("line1068\n"); fflush(stdout);
1029 return 0;
1030 mbobra 1.1
1031 }
1032
1033
1034 /*
1035 * Compute space weather indices
1036 */
1037
1038 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1039 {
1040
1041 int status = 0;
1042 int nx = mInfo->ncol, ny = mInfo->nrow;
1043 int nxny = nx * ny;
1044 int dims[2] = {nx, ny};
1045
1046 // Get bx, by, bz, mask
1047
1048 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1049 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1050 DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1051 mbobra 1.1 int *bitmask = (int *) bitmaskArray->data; // get the previously made mask array
1052
1053 //Use magnetogram map to compute R
1054 DRMS_Segment_t *losSeg = drms_segment_lookup(inRec, "magnetogram");
1055 DRMS_Array_t *losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status);
1056 float *los = (float *) losArray->data; // los
1057
1058 // Get emphemeris
1059 float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);
1060 float dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);
1061 double rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
1062 double rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status);
1063 float imcrpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1064 float imcrpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1065 float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);
1066 float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);
1067
1068 // convert cdelt1_orig from degrees to arcsec
1069 float cdelt1 = (atan((rsun_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
1070
1071 // Temp arrays
1072 mbobra 1.1 float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1073 float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1074
1075 // define some values for the R calculation
1076 int scale = round(2.0/cdelt1);
1077 int nx1 = nx/scale;
1078 int ny1 = ny/scale;
1079 int nxp = nx1+40; // same comment as above
1080 int nyp = ny1+40; // why is this a +40 pixel size? is this an MDI pixel?
1081 float *rim = (float *)malloc(nx1*ny1*sizeof(float));
1082 float *p1p0 = (float *)malloc(nx1*ny1*sizeof(float));
1083 float *p1n0 = (float *)malloc(nx1*ny1*sizeof(float));
1084 float *p1p = (float *)malloc(nx1*ny1*sizeof(float));
1085 float *p1n = (float *)malloc(nx1*ny1*sizeof(float));
1086 float *p1 = (float *)malloc(nx1*ny1*sizeof(float));
1087 float *pmap = (float *)malloc(nxp*nyp*sizeof(float));
1088 float *p1pad = (float *)malloc(nxp*nyp*sizeof(float));
1089 float *pmapn = (float *)malloc(nx1*ny1*sizeof(float));
1090
1091
1092 // THREE spaceweather quantities computed: USFLUX, MEANGBZ, R_VALUE
1093 mbobra 1.1 if (computeAbsFlux(los, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1094 &(swKeys_ptr->count_mask), bitmask, cdelt1, rsun_ref, rsun_obs))
1095 {
1096 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1097 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1098 swKeys_ptr->count_mask = DRMS_MISSING_INT;
1099 }
1100
1101
1102 if (computeBzderivative(los, dims, &(swKeys_ptr->mean_derivative_bz),
1103 bitmask, derx_bz, dery_bz))
1104 {
1105 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1106 }
1107
1108
1109 if (computeR(los, dims, &(swKeys_ptr->Rparam), cdelt1, rim, p1p0, p1n0,
1110 p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
1111 {
1112 swKeys_ptr->Rparam = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1113 }
1114 mbobra 1.1
1115
1116 // Clean up the arrays
1117 drms_free_array(bitmaskArray);
1118 //drms_free_array(bzArray);
1119 drms_free_array(losArray);
1120
1121 // free arrays related to Bz derivative
1122 free(derx_bz); free(dery_bz);
1123 // free the arrays that are related to the r calculation
1124 free(rim);
1125 free(p1p0);
1126 free(p1n0);
1127 free(p1p);
1128 free(p1n);
1129 free(p1);
1130 free(pmap);
1131 free(p1pad);
1132 free(pmapn);
1133 }
1134
1135 mbobra 1.1 /*
1136 * Set space weather indices
1137 */
1138
1139 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1140 {
1141 drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1142 drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
1143 drms_setkey_float(outRec, "R_VALUE", swKeys_ptr->Rparam);
1144 drms_setkey_float(outRec, "CMASK", swKeys_ptr->count_mask);
1145 };
1146
1147 /*
1148 * Set all keywords, no error checking for now
1149 */
1150
1151 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *mharpRec, struct mapInfo *mInfo)
1152 {
1153
1154 int status = 0;
1155
1156 mbobra 1.1 // Change a few geometry keywords for CEA & cutout records
1157 if (mInfo != NULL) { // CEA
1158 printf("Calculating CEA keys\n");
1159 drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1160 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1161 drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1162 drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1163 drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1164 drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1165 drms_setkey_string(outRec, "CUNIT1", "degree");
1166 drms_setkey_string(outRec, "CUNIT2", "degree");
1167 char key[64];
1168 snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1169 drms_setkey_string(outRec, "CTYPE1", key);
1170 snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
1171 drms_setkey_string(outRec, "CTYPE2", key);
1172 drms_setkey_float(outRec, "CROTA2", 0.0);
1173 // Set BUNIT for each segment
1174 int nSeg = 3;
1175 for (int iSeg = 0; iSeg < nSeg; iSeg++) {
1176 DRMS_Segment_t *outSeg = NULL;
1177 mbobra 1.1 outSeg = drms_segment_lookup(outRec, CEASegs[iSeg]);
1178 if (!outSeg) continue;
1179 char bunit_xxx[20];
1180 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1181 //printf("%s, %s\n", bunit_xxx, CEABunits[iSeg]);
1182 drms_setkey_string(outRec, bunit_xxx, CEABunits[iSeg]);
1183 }
1184
1185 } else { // Cutout
1186
1187 float disk_xc, disk_yc;
1188 disk_xc = drms_getkey_float(mharpRec, "IMCRPIX1", &status);
1189 disk_yc = drms_getkey_float(mharpRec, "IMCRPIX2", &status);
1190 float x_ll = drms_getkey_float(mharpRec, "CRPIX1", &status);
1191 float y_ll = drms_getkey_float(mharpRec, "CRPIX2", &status);
1192 // Defined as disk center's pixel address wrt lower-left of cutout
1193 drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
1194 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
1195 // Always 0.
1196 drms_setkey_float(outRec, "CRVAL1", 0);
1197 drms_setkey_float(outRec, "CRVAL2", 0);
1198 mbobra 1.1
1199 // Jan 2 2014 XS
1200 int nSeg = ARRLENGTH(CutSegs);
1201 for (int iSeg = 0; iSeg < nSeg; iSeg++)
1202 {
1203 DRMS_Segment_t *outSeg = NULL;
1204 outSeg = drms_segment_lookup(outRec, CutSegs[iSeg]);
1205 if (!outSeg) continue;
1206 // Set Bunit
1207 char bunit_xxx[20];
1208 sprintf(bunit_xxx, "BUNIT_%03d", iSeg);
1209 //printf("%s, %s\n", bunit_xxx, CutBunits[iSeg]);
1210 drms_setkey_string(outRec, bunit_xxx, CutBunits[iSeg]);
1211 }
1212 }
1213
1214
1215 TIME val, trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1216 tnow = (double)time(NULL);
1217 tnow += UNIX_epoch;
1218
1219 mbobra 1.1 val = drms_getkey_time(mharpRec, "DATE", &status);
1220 drms_setkey_time(outRec, "DATE", tnow);
1221
1222 // set cvs commit version into keyword CODEVER7
1223 char *cvsinfo = strdup("$Id");
1224 char *cvsinfo2 = smarp_functions_version();
1225 char cvsinfoall[2048];
1226 strcat(cvsinfoall,cvsinfo);
1227 strcat(cvsinfoall,"\n");
1228 strcat(cvsinfoall,cvsinfo2);
1229 status = drms_setkey_string(outRec, "CODEVER7", cvsinfoall);
1230
1231 };
1232
1233 /* ############# Nearest neighbour interpolation ############### */
1234
1235 float nnb (float *f, int nx, int ny, double x, double y)
1236 {
1237
1238 if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
1239 return DRMS_MISSING_FLOAT;
1240 mbobra 1.1 int ilow = floor (x);
1241 int jlow = floor (y);
1242 int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
1243 int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
1244 return f[j * nx + i];
1245
1246 }
1247
1248 /* ################## Wrapper for Jesper's rebin code ################## */
1249
1250 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
1251 {
1252
1253 struct fresize_struct fresizes;
1254 int nxout, nyout, xoff, yoff;
1255 int nlead = nx;
1256
1257 nxout = nx / nbin; nyout = ny / nbin;
1258 if (gauss && nbin != 1)
1259 init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin); // for nbin=3, sigma=1, half truncate width=2
1260 else
1261 mbobra 1.1 init_fresize_bin(&fresizes, nbin);
1262 xoff = nbin / 2 + nbin / 2;
1263 yoff = nbin / 2 + nbin / 2;
1264 fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
1265
1266 }
|
Return to
smarp.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps