36 xudong 1.1 *
37 * Example:
38 * sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
39 "bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
40 "dop=hmi.V_720s[2012.02.20_10:00]" \
41 "cont=hmi.Ic_720s[2012.02.20_10:00]" \
42 "sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
43 * For comparison:
44 * bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
45 "out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
46 *
47 *
48 */
49
50 #include <stdio.h>
51 #include <stdlib.h>
52 #include <time.h>
53 #include <sys/time.h>
54 #include <math.h>
55 #include <string.h>
56 #include "jsoc_main.h"
57 xudong 1.1 #include "astro.h"
58 #include "fstats.h"
59 #include "cartography.c"
60 #include "fresize.h"
61 #include "finterpolate.h"
62 #include "img2helioVector.c"
63 #include "copy_me_keys.c"
64 #include "errorprop.c"
65 #include "sw_functions.c"
66
67 #include <mkl_blas.h>
68 #include <mkl_service.h>
69 #include <mkl_lapack.h>
70 #include <mkl_vml_functions.h>
71 #include <omp.h>
72
73 #define PI (M_PI)
74 #define RADSINDEG (PI/180.)
75 #define RAD2ARCSEC (648000./M_PI)
76 #define SECINDAY (86400.)
77 #define FOURK (4096)
78 xudong 1.1 #define FOURK2 (16777216)
79
80 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
81
82 // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
83 #define NYQVIST (0.015)
84
85 #ifndef MIN
86 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
87 #endif
88 #ifndef MAX
89 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
90 #endif
91
92 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
93 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
94
95 #define kNotSpecified "Not Specified"
96
97 // Macros for WCS transformations. assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
98 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
99 xudong 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.
100 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
101 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
102 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
103 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
104
105 #define DISAMB_AZI 1
106 #define XSCALE 0.03
107 #define YSCALE 0.03
108 #define NBIN 3
109 #define INTERP 0
110 #define dpath "/home/jsoc/cvs/Development/JSOC"
111
112 /* ========================================================================================================== */
113
114 // Space weather keywords
115 struct swIndex {
116 float mean_vf;
117 float absFlux;
118 float mean_hf;
119 float mean_gamma;
120 xudong 1.1 float mean_derivative_btotal;
121 float mean_derivative_bh;
122 float mean_derivative_bz;
123 float mean_jz;
124 float us_i;
125 float mean_alpha;
126 float mean_ih;
127 float total_us_ih;
128 float total_abs_ih;
129 float totaljz;
130 float totpot;
131 float meanpot;
132 float area_w_shear_gt_45;
133 float meanshear_angle;
134 float area_w_shear_gt_45h;
135 float meanshear_angleh;
136 };
137
138 // Mapping method
139 enum projection {
140 carree,
141 xudong 1.1 cassini,
142 mercator,
143 cyleqa,
144 sineqa,
145 gnomonic,
146 postel,
147 stereographic,
148 orthographic,
149 lambert
150 };
151
152 // Ephemeris
153 struct ephemeris {
154 double disk_lonc, disk_latc;
155 double disk_xc, disk_yc;
156 double rSun, asd, pa;
157 };
158
159 // Mapping information
160 struct mapInfo {
161 float xc, yc; // reference point: center
162 xudong 1.1 int nrow, ncol; // size
163 float xscale, yscale; // scale
164 int nbin;
165 enum projection proj; // projection method
166 struct ephemeris ephem; // ephemeris info
167 float *xi_out, *zeta_out; // coordinate on full disk image to sample at
168 };
169
170
171 /* ========================================================================================================== */
172
173 /* Get all input data series */
174 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
175 char *mharpQuery, char *bharpQuery);
176
177 /* Check if Mharp and Bharp match */
178 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS);
179
180 /* Get other data series */
181 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
182
183 xudong 1.1 /* Find record from record set with given T_rec */
184 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
185
186 // ===================
187
188 /* Create CEA record */
189 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
190 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
191 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
192
193 /* Mapping single segment, wrapper */
194 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
195 struct mapInfo *mInfo, char *segName);
196
197 /* Mapping vector magnetogram, wrapper */
198 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
199
200 /* Mapping errors of vector magnetogram, wraper */
201 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
202
203 /* Determine reference point coordinate and patch size according to input */
204 xudong 1.1 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
205
206 /* Get ephemeris information */
207 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
208
209 /* Compute the coordinates at which the full disk image is sampled */
210 void findCoord(struct mapInfo *mInfo);
211
212 /* Mapping function */
213 int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
214
215 /* Performing local vector transformation */
216 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
217
218 /* Map and propogate errors */
219 int getBErr(float *bx_err, float *by_err, float *bz_err,
220 DRMS_Record_t *inRec, struct mapInfo *mInfo);
221
222 /* Read full disk vector magnetogram */
223 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
224
225 xudong 1.1 /* Read variances and covariances of vector magnetograms */
226 int readVectorBErr(DRMS_Record_t *bharpRec,
227 float *bT, float *bI, float *bA,
228 float *errbT, float *errbI, float *errbA,
229 float *errbTbI, float *errbTbA, float *errbIbA);
230
231 // ===================
232
233 /* Create Cutout record */
234 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
235 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
236 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
237
238 /* Get cutout and write segment */
239 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
240
241 // ===================
242
243 /* Compute space weather indices, no error checking for now */
244 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
245
246 xudong 1.1 /* Set space weather indices, no error checking for now */
247 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
248
249 /* Set all keywords, no error checking for now */
250 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec);
251
252 // ===================
253
254 /* Nearest neighbor interpolation */
255 float nnb (float *f, int nx, int ny, double x, double y);
256
257 /* Wrapper for Jesper's rebin code */
258 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
259
260 /* ========================================================================================================== */
261
262 /* Remap segment names */
263 #define BR_SEG_CEA "Br"
264 #define BT_SEG_CEA "Bt"
265 #define BP_SEG_CEA "Bp"
266 #define BR_ERR_SEG_CEA "Br_err"
267 xudong 1.1 #define BT_ERR_SEG_CEA "Bt_err"
268 #define BP_ERR_SEG_CEA "Bp_err"
269
270 /* Cutout segment names, input identical to output */
271 char *MharpSegs[] = {"magnetogram", "bitmap"};
272 char *BharpSegs[] = {"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
273 "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
274 "conv_flag", // fixed
275 "info_map", "confid_map",
276 "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
277 "field_inclination_err", "field_az_err", "inclin_azimuth_err",
278 "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
279 "disambig", "conf_disambig"};
280 // For stats
281 char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
282 "inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
283 "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
284 "conv_flag", // fixed
285 "info_map", "confid_map",
286 "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
287 "field_inclination_err", "field_az_err", "inclin_azimuth_err",
288 xudong 1.1 "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
289 "disambig", "conf_disambig"};
290 char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
291 BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
292 /* ========================================================================================================== */
293
294
295
296 char *module_name = "sharp";
297 char *version_id = "2012 Jul 02"; /* Version number */
298
299 ModuleArgs_t module_args[] =
300 {
301 {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
302 {ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
303 {ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
304 {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
305 {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
306 {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
307 {ARG_END}
308 };
309 xudong 1.1
310 int DoIt(void)
311 {
312
313 int status = DRMS_SUCCESS;
314 int nrecs, irec;
315
316 char *mharpQuery, *bharpQuery;
317 char *dopQuery, *contQuery;
318 char *sharpCeaQuery, *sharpCutQuery;
319
320 DRMS_RecordSet_t *mharpRS = NULL, *bharpRS = NULL;
321 DRMS_RecordSet_t *dopRS = NULL, *contRS = NULL;
322
323 /* Get parameters */
324
325 mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
326 bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
327 dopQuery = (char *) params_get_str(&cmdparams, "dop");
328 contQuery = (char *) params_get_str(&cmdparams, "cont");
329 sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
330 xudong 1.1 sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
331
332 /* Get input data, check everything */
333
334 if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
335 DIE("Input harp data error.");
336 nrecs = mharpRS->n;
337
338 if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
339 DIE("Input doppler data error.");
340
341 if (getInputRS_aux(&contRS, contQuery, mharpRS))
342 DIE("Input continuum data error.");
343
344 /* Start */
345
346 printf("==============\nStart. %d image(s) in total.\n", nrecs);
347
348 for (irec = 0; irec < nrecs; irec++) {
349
350 /* Records in work */
351 xudong 1.1
352 DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
353 mharpRec = mharpRS->records[irec];
354 bharpRec = bharpRS->records[irec];
355
356 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
357
358 struct swIndex swKeys;
359
360 DRMS_Record_t *dopRec = NULL, *contRec = NULL;
361 if (getInputRec_aux(&dopRec, dopRS, trec)) {
362 printf("Fetching Doppler failed, image #%d skipped.\n", irec);
363 continue;
364 }
365 if (getInputRec_aux(&contRec, contRS, trec)) {
366 printf("Fetching continuum failed, image #%d skipped.\n", irec);
367 continue;
368 }
369
370 /* Create CEA record */
371
372 xudong 1.1 DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
373 if (status) { // if failed
374 printf("Creating CEA failed, image #%d skipped.\n", irec);
375 continue;
376 }
377
378 if (createCeaRecord(mharpRec, bharpRec, dopRec, contRec, sharpCeaRec, &swKeys)) { // do the work
379 printf("Creating CEA failed, image #%d skipped.\n", irec);
380 drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
381 continue;
382 } // swKeys updated here
383
384 drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
385
386 /* Create Cutout record */
387
388 DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
389 if (status) { // if failed
390 printf("Creating cutout failed, image #%d skipped.\n", irec);
391 continue;
392 }
393 xudong 1.1
394 if (createCutRecord(mharpRec, bharpRec, dopRec, contRec, sharpCutRec, &swKeys)) { // do the work
395 printf("Creating cutout failed, image #%d skipped.\n", irec);
396 drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
397 continue;
398 } // swKeys used here
399
400 drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
401
402 /* Done */
403
404 printf("Image #%d done.\n", irec);
405
406 } // irec
407
408 drms_close_records(mharpRS, DRMS_FREE_RECORD);
409 drms_close_records(bharpRS, DRMS_FREE_RECORD);
410
411 return 0;
412
413 } // DoIt
414 xudong 1.1
415
416 // ===================================================================
417 // ===================================================================
418 // ===================================================================
419
420
421 /*
422 * Get input data series, including mHarp and bharp
423 * Need all records to match, otherwise quit
424 *
425 */
426
427 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
428 char *mharpQuery, char *bharpQuery)
429 {
430
431 int status = 0;
432
433 *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
434 if (status || (*mharpRS_ptr)->n == 0) return 1;
435 xudong 1.1
436 *bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
437 if (status || (*bharpRS_ptr)->n == 0) return 1;
438
439 if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
440
441 return 0;
442
443 }
444
445 /*
446 * Check if Mharp and Bharp match
447 *
448 */
449
450 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS)
451 {
452
453 int status = 0;
454 int nrecs = mharpRS->n;
455
456 xudong 1.1 DRMS_Record_t *mharpRec_t = NULL, *bharpRec_t = NULL; // temporary recs for utility
457
458 if (bharpRS->n != nrecs) {
459 return 1; // return 1 if different
460 }
461
462 for (int i = 0; i < nrecs; i++) {
463 mharpRec_t = mharpRS->records[i];
464 bharpRec_t = bharpRS->records[i];
465 if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
466 drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
467 (drms_getkey_time(mharpRec_t, "T_REC", &status) !=
468 drms_getkey_time(bharpRec_t, "T_REC", &status)))
469 {
470 return 1;
471 }
472 }
473
474 return 0;
475
476 }
477 xudong 1.1
478 /*
479 * Get other data series, check all T_REC are available
480 *
481 */
482
483 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
484 {
485
486 int status = 0;
487
488 *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
489 if (status || (*inRS_ptr)->n == 0) return status;
490
491 // Check if all T_rec are available, need to match both ways
492 int n = harpRS->n, n0 = (*inRS_ptr)->n;
493
494 for (int i0 = 0; i0 < n0; i0++) {
495 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
496 TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
497 TIME trec = 0;
498 xudong 1.1 for (int i = 0; i < n; i++) {
499 DRMS_Record_t *harpRec = harpRS->records[i];
500 trec = drms_getkey_time(harpRec, "T_REC", &status);
501 if (fabs(trec0 - trec) < 10) break;
502 }
503 if (fabs(trec0 - trec) >= 10) return 1;
504 }
505
506 for (int i = 0; i < n; i++) {
507 DRMS_Record_t *harpRec = harpRS->records[i];
508 TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
509 TIME trec0 = 0;
510 for (int i0 = 0; i0 < n0; i0++) {
511 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
512 trec0 = drms_getkey_time(inRec, "T_REC", &status);
513 if (fabs(trec0 - trec) < 10) break;
514 }
515 if (fabs(trec0 - trec) >= 10) return 1;
516 }
517
518 return 0;
519 xudong 1.1
520 }
521
522 /*
523 * Find record from record set with given T_rec
524 *
525 */
526
527 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
528 {
529
530 int status = 0;
531
532 int n = inRS->n;
533 for (int i = 0; i < n; i++) {
534 *inRec_ptr = inRS->records[i];
535 TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
536 if (fabs(trec0 - trec) < 10) return 0;
537 }
538
539 return 1;
540 xudong 1.1
541 }
542
543
544
545
546 /*
547 * Create CEA record: top level subroutine
548 * Also compute all the space weather keywords here
549 *
550 */
551
552 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
553 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
554 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
555 {
556
557 int status = 0;
558 DRMS_Segment_t *inSeg;
559 DRMS_Array_t *inArray;
560
561 xudong 1.1 struct mapInfo mInfo;
562 mInfo.proj = (enum projection) cyleqa; // projection method
563 mInfo.xscale = XSCALE;
564 mInfo.yscale = YSCALE;
565 mInfo.nbin = NBIN;
566
567 // Get ephemeris
568
|
651 xudong 1.1
652 // Space weather
653
654 computeSWIndex(swKeys_ptr, sharpRec, &mInfo); // compute it!
655 printf("Space weather indices done.\n");
656
657 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
658
659 // Stats
660
661 int nCEASegs = ARRLENGTH(CEASegs);
662 for (int iSeg = 0; iSeg < nCEASegs; iSeg++) {
663 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
664 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
665 int stat = set_statistics(outSeg, outArray, 1);
666 // printf("%d => %d\n", iSeg, stat);
667 drms_free_array(outArray);
668 }
669
670 free(mInfo.xi_out);
671 free(mInfo.zeta_out);
672 xudong 1.1 return 0;
673
674 }
675
676
677 /*
678 * Mapping a single segment
679 * Read in full disk image, utilize mapImage for mapping
680 * then write the segment out, segName same in in/out Rec
681 *
682 */
683
684 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
685 struct mapInfo *mInfo, char *segName)
686 {
687
688 int status = 0;
689 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
690 int dims[2] = {nx, ny};
691
692 // Input full disk array
693 xudong 1.1
694 DRMS_Segment_t *inSeg = NULL;
695 inSeg = drms_segment_lookup(inRec, segName);
696 if (!inSeg) return 1;
697
698 DRMS_Array_t *inArray = NULL;
699 inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
700 if (!inArray) return 1;
701
702 float *inData;
703 int xsz = inArray->axis[0], ysz = inArray->axis[1];
704 if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk
705 float *inData0 = (float *) inArray->data;
706 inData = (float *) (calloc(FOURK2, sizeof(float)));
707 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
708 int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
709 int ind_map;
710 for (int row = 0; row < ysz; row++) {
711 for (int col = 0; col < xsz; col++) {
712 ind_map = (row + y0) * FOURK + (col + x0);
713 inData[ind_map] = inData0[row * xsz + col];
714 xudong 1.1 }
715 }
716 drms_free_array(inArray); inArray = NULL;
717 } else {
718 inData = (float *) inArray->data;
719 }
720
721 // Mapping
722
723 float *map = (float *) (malloc(nxny * sizeof(float)));
724 if (performSampling(map, inData, mInfo))
725 {if (inArray) drms_free_array(inArray); free(map); return 1;}
726
727 // Write out
728
729 DRMS_Segment_t *outSeg = NULL;
730 outSeg = drms_segment_lookup(sharpRec, segName);
731 if (!outSeg) return 1;
732
733 DRMS_Type_t arrayType = outSeg->info->type;
734 DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
735 xudong 1.1 if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
736
737 // convert to needed data type
738
739 drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray);
740
741 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
742 outArray->parent_segment = outSeg;
743 outArray->israw = 0; // always compressed
744 outArray->bzero = outSeg->bzero;
745 outArray->bscale = outSeg->bscale;
746
747 status = drms_segment_write(outSeg, outArray, 0);
748 if (status) return 0;
749
750 if (inArray) drms_free_array(inArray);
751 if (outArray) drms_free_array(outArray);
752 return 0;
753
754 }
755
756 xudong 1.1
757 /*
758 * Mapping vector magnetogram
759 *
760 */
761
762 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
763 {
764
765 int status = 0;
766 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
767 int dims[2] = {nx, ny};
768
769 // Read in segments, filling factor assume to be 1
770
771 float *bx_img = (float *) (malloc(FOURK2 * sizeof(float)));
772 float *by_img = (float *) (malloc(FOURK2 * sizeof(float)));
773 float *bz_img = (float *) (malloc(FOURK2 * sizeof(float)));
774
775 if (readVectorB(bharpRec, bx_img, by_img, bz_img)) {
776 printf("Read full disk image error\n");
777 xudong 1.1 free(bx_img); free(by_img); free(bz_img);
778 return 1;
779 }
780
781 // Mapping
782
783 float *bx_map = NULL, *by_map = NULL, *bz_map = NULL; // intermediate maps, in CCD bxyz representation
784
785 bx_map = (float *) (malloc(nxny * sizeof(float)));
786 if (performSampling(bx_map, bx_img, mInfo))
787 {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
788
789 by_map = (float *) (malloc(nxny * sizeof(float)));
790 if (performSampling(by_map, by_img, mInfo))
791 {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
792
793 bz_map = (float *) (malloc(nxny * sizeof(float)));
794 if (performSampling(bz_map, bz_img, mInfo))
795 {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
796
797 free(bx_img); free(by_img); free(bz_img);
798 xudong 1.1
799 // Vector transform
800
801 vectorTransform(bx_map, by_map, bz_map, mInfo);
802
803 for (int i = 0; i < nxny; i++) by_map[i] *= -1; // positive theta pointing south
804
805 // Write out
806
807 DRMS_Segment_t *outSeg;
808 DRMS_Array_t *outArray;
809
810 float *data_prt[3] = {bx_map, by_map, bz_map};
811 char *segName[3] = {BP_SEG_CEA, BT_SEG_CEA, BR_SEG_CEA};
812
813 for (int iSeg = 0; iSeg < 3; iSeg++) {
814 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
815 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
816 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
817 outArray->parent_segment = outSeg;
818 outArray->israw = 0;
819 xudong 1.1 outArray->bzero = outSeg->bzero;
820 outArray->bscale = outSeg->bscale;
821 status = drms_segment_write(outSeg, outArray, 0);
822 if (status) return 1;
823 drms_free_array(outArray);
824 }
825
826 //
827
828 return 0;
829
830 }
831
832
833 /*
834 * Mapping vector magnetogram errors
835 *
836 */
837
838 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
839 {
840 xudong 1.1
841 int status = 0;
842
843 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
844 int dims[2] = {nx, ny};
845
846 // Compute propogated errors, using nearest neighbour interpolation
847
848 float *bx_err = (float *) (malloc(nxny * sizeof(float)));
849 float *by_err = (float *) (malloc(nxny * sizeof(float)));
850 float *bz_err = (float *) (malloc(nxny * sizeof(float)));
851
852 if (getBErr(bx_err, by_err, bz_err, bharpRec, mInfo)) {
853 free(bx_err); free(by_err); free(bz_err);
854 return 1;
855 }
856
857 // Write out
858
859 DRMS_Segment_t *outSeg;
860 DRMS_Array_t *outArray;
861 xudong 1.1
862 float *data_prt[3] = {bx_err, by_err, bz_err};
863 char *segName[3] = {BP_ERR_SEG_CEA, BT_ERR_SEG_CEA, BR_ERR_SEG_CEA};
864
865 for (int iSeg = 0; iSeg < 3; iSeg++) {
866 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
867 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
868 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
869 outArray->parent_segment = outSeg;
870 outArray->israw = 0;
871 outArray->bzero = outSeg->bzero;
872 outArray->bscale = outSeg->bscale;
873 status = drms_segment_write(outSeg, outArray, 0);
874 if (status) return 1;
875 drms_free_array(outArray);
876 }
877
878 //
879
880 return 0;
881
882 xudong 1.1 }
883
884
885
886 /*
887 * Determine reference point coordinate and patch size according to keywords
888 * xc, yc are the coordinate of patch center, in degrees
889 * ncol and nrow are the final size
890 *
891 */
892
893 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
894 {
895
896 int status = 0;
897 int harpnum = drms_getkey_int(inRec, "HARPNUM", &status);
898 TIME trec = drms_getkey_time(inRec, "T_REC", &status);
899 float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
900
901 /* Center coord */
902
903 xudong 1.1 float minlon = drms_getkey_float(inRec, "LONDTMIN", &status); if (status) return 1; // Stonyhurst lon
904 float maxlon = drms_getkey_float(inRec, "LONDTMAX", &status); if (status) return 1;
905 float minlat = drms_getkey_float(inRec, "LATDTMIN", &status); if (status) return 1;
906 float maxlat = drms_getkey_float(inRec, "LATDTMAX", &status); if (status) return 1;
907
908 // A bug fixer for HARP (per M. Turmon)
909 // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
910 // Also keywords such as "SIZE" will be NaN
911 // We compute minlon & minlat then by
912 // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
913
914 float psize = drms_getkey_float(inRec, "SIZE", &status);
915 if (psize != psize) {
916 TIME t0 = drms_getkey_time(inRec, "T_FRST", &status); if (status) return 1;
917 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
918 char firstRecQuery[100], t0_str[100];
919 sprint_time(t0_str, t0, "TAI", 0);
920 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
921 DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
922 if (status || tmpRS->n != 1) return 1;
923 DRMS_Record_t *tmpRec = tmpRS->records[0];
924 xudong 1.1 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
925 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
926 minlon = minlon0 + (trec - t0) * omega / SECINDAY;
927 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
928 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
929 }
930
931 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
932 mInfo->yc = (maxlat + minlat) / 2.;
933
934 /* Size */
935
936 mInfo->ncol = round((maxlon - minlon) / mInfo->xscale);
937 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
938
939 return 0;
940
941 }
942
943
944 /*
945 xudong 1.1 * Fetch ephemeris info from a DRMS record
946 * No error checking for now
947 *
948 */
949
950 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
951 {
952
953 int status = 0;
954
955 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
956 double sina = sin(crota2 * RADSINDEG);
957 double cosa = cos(crota2 * RADSINDEG);
958
959 ephem->pa = - crota2 * RADSINDEG;
960 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
961 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
962
963 float crvalx = 0.0;
964 float crvaly = 0.0;
965 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
966 xudong 1.1 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
967 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
968 printf("cdelt=%f\n",cdelt);
969 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0
970 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
971
972 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
973 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
974 if (status) rSun_ref = 6.96e8;
975
976 ephem->asd = asin(rSun_ref/dSun);
977 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
978
979 return 0;
980
981 }
982
983
984 /*
985 * Compute the coordinates to be sampled on full disk image
986 * mInfo->xi_out & mInfo->zeta_out
987 xudong 1.1 * This is oversampled, its size is ncol0 & nrow0 as shown below
988 *
989 *
990 */
991
992 void findCoord(struct mapInfo *mInfo)
993 {
994
995 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
996 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
997
998 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution
999 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad
1000
1001 double lonc = mInfo->xc * RADSINDEG; // in rad
1002 double latc = mInfo->yc * RADSINDEG;
1003
1004 double disk_lonc = (mInfo->ephem).disk_lonc;
1005 double disk_latc = (mInfo->ephem).disk_latc;
1006
1007 double rSun = (mInfo->ephem).rSun;
1008 xudong 1.1 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1009 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1010 double pa = (mInfo->ephem).pa;
1011
1012 // Temp pointers
1013
1014 float *xi_out = mInfo->xi_out;
1015 float *zeta_out = mInfo->zeta_out;
1016
1017 // start
1018
1019 double x, y; // map coord
1020 double lat, lon; // helio coord
1021 double xi, zeta; // image coord (for one point)
1022
1023 int ind_map;
1024
1025 for (int row0 = 0; row0 < nrow0; row0++) {
1026 for (int col0 = 0; col0 < ncol0; col0++) {
1027
1028 ind_map = row0 * ncol0 + col0;
1029 xudong 1.1
1030 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad
1031 y = (row0 + 0.5 - nrow0/2.) * yscale0;
1032
1033 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1034 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1035 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1036 * is the heliographic longitude and latitude of the map center. Both are in degree.
1037 */
1038
1039 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1040 xi_out[ind_map] = -1;
1041 zeta_out[ind_map] = -1;
1042 continue;
1043 }
1044
1045 /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
1046 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1047 */
1048
1049 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1050 xudong 1.1 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1051 xi_out[ind_map] = -1;
1052 zeta_out[ind_map] = -1;
1053 continue;
1054 }
1055
1056 xi_out[ind_map] = xi * rSun;
1057 zeta_out[ind_map] = zeta * rSun;
1058
1059 }
1060 }
1061
1062 }
1063
1064
1065 /*
1066 * Sampling function
1067 * oversampling by nbin, then binning using a Gaussian
1068 * save results in outData, always of float type
1069 *
1070 */
1071 xudong 1.1
1072 int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1073 {
1074
1075 int status = 0;
1076
1077 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
1078 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1079
1080 float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1081
1082 float *xi_out = mInfo->xi_out;
1083 float *zeta_out = mInfo->zeta_out;
1084
1085 // Interpolation
1086
1087 struct fint_struct pars;
1088 int interpOpt = INTERP; // Use Wiener by default, 6 order, 1 constraint
1089
1090 switch (interpOpt) {
1091 case 0: // Wiener, 6 order, 1 constraint
1092 xudong 1.1 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
1093 break;
1094 case 1: // Cubic convolution
1095 init_finterpolate_cubic_conv(&pars, 1., 3.);
1096 break;
1097 case 2: // Bilinear
1098 init_finterpolate_linear(&pars, 1.);
1099 break;
1100 default:
1101 return 1;
1102 }
1103
1104 finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1105 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1106
1107 // Rebinning, smoothing
1108
1109 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
1110
1111 //
1112
1113 xudong 1.1 return 0;
1114
1115 }
1116
1117
1118 /*
1119 * Performing local vector transformation
1120 * xyz: z refers to vertical (radial) component, x EW (phi), y NS
1121 *
1122 */
1123
1124 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo)
1125 {
1126
1127 int ncol = mInfo->ncol;
1128 int nrow = mInfo->nrow;
1129
1130 float xscale = mInfo->xscale * RADSINDEG; // in rad
1131 float yscale = mInfo->yscale * RADSINDEG;
1132
1133 double lonc = mInfo->xc * RADSINDEG; // in rad
1134 xudong 1.1 double latc = mInfo->yc * RADSINDEG;
1135
1136 double disk_lonc = (mInfo->ephem).disk_lonc;
1137 double disk_latc = (mInfo->ephem).disk_latc;
1138
1139 double rSun = (mInfo->ephem).rSun;
1140 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1141 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1142 double pa = (mInfo->ephem).pa;
1143
1144 int ind_map;
1145 double x, y;
1146 double lat, lon; // lat / lon for current point
1147
1148 double bx_tmp, by_tmp, bz_tmp;
1149
1150 //
1151
1152 for (int row = 0; row < mInfo->nrow; row++) {
1153 for (int col = 0; col < mInfo->ncol; col++) {
1154
1155 xudong 1.1 ind_map = row * mInfo->ncol + col;
1156
1157 x = (col + 0.5 - ncol / 2.) * xscale;
1158 y = (row + 0.5 - nrow / 2.) * yscale;
1159
1160 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1161 bx_map[ind_map] = DRMS_MISSING_FLOAT;
1162 by_map[ind_map] = DRMS_MISSING_FLOAT;
1163 bz_map[ind_map] = DRMS_MISSING_FLOAT;
1164 continue;
1165 }
1166
1167 bx_tmp = by_tmp = bz_tmp = 0;
1168
1169 img2helioVector (bx_map[ind_map], by_map[ind_map], bz_map[ind_map],
1170 &bx_tmp, &by_tmp, &bz_tmp,
1171 lon, lat, disk_lonc, disk_latc, pa);
1172
1173 bx_map[ind_map] = bx_tmp;
1174 by_map[ind_map] = by_tmp;
1175 bz_map[ind_map] = bz_tmp;
1176 xudong 1.1
1177 }
1178 }
1179
1180 }
1181
1182
1183
1184 /*
1185 * Map and propogate vector field errors
1186 *
1187 */
1188
1189 int getBErr(float *bx_err, float *by_err, float *bz_err,
1190 DRMS_Record_t *inRec, struct mapInfo *mInfo)
1191 {
1192
1193 int status = 0;
1194
1195 // Get variances and covariances, filling factor assume to be 1
1196
1197 xudong 1.1 float *bT = (float *) (malloc(FOURK2 * sizeof(float))); // field
1198 float *bI = (float *) (malloc(FOURK2 * sizeof(float))); // inclination
1199 float *bA = (float *) (malloc(FOURK2 * sizeof(float))); // azimuth
1200
1201 float *errbT = (float *) (malloc(FOURK2 * sizeof(float)));
1202 float *errbI = (float *) (malloc(FOURK2 * sizeof(float)));
1203 float *errbA = (float *) (malloc(FOURK2 * sizeof(float)));
1204
1205 float *errbTbI = (float *) (malloc(FOURK2 * sizeof(float)));
1206 float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1207 float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1208
1209 if (readVectorBErr(inRec,
1210 bT, bI, bA,
1211 errbT, errbI, errbA,
1212 errbTbI, errbTbA, errbIbA)) {
1213 printf("Read full disk variances & covariances error\n");
1214 free(bT); free(bI); free(bA);
1215 free(errbT); free(errbI); free(errbA);
1216 free(errbTbI); free(errbTbA); free(errbIbA);
1217 return 1;
1218 xudong 1.1 }
1219
1220 // Size
1221
1222 int ncol = mInfo->ncol;
1223 int nrow = mInfo->nrow;
1224
1225 float xscale = mInfo->xscale * RADSINDEG; // in rad
1226 float yscale = mInfo->yscale * RADSINDEG;
1227
1228 double lonc = mInfo->xc * RADSINDEG; // in rad
1229 double latc = mInfo->yc * RADSINDEG;
1230
1231 double disk_lonc = (mInfo->ephem).disk_lonc;
1232 double disk_latc = (mInfo->ephem).disk_latc;
1233
1234 double rSun = (mInfo->ephem).rSun;
1235 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1236 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1237 double pa = (mInfo->ephem).pa;
1238
1239 xudong 1.1 // Start
1240
1241 double x, y; // map coord
1242 double lat, lon; // spherical coord
1243 double xi, zeta; // image coord, round to full pixel
1244
1245 int ind_map, ind_img;
1246
1247 double bpSigma2, btSigma2, brSigma2; // variances after prop
1248
1249 for (int row = 0; row < mInfo->nrow; row++) {
1250 for (int col = 0; col < mInfo->ncol; col++) {
1251
1252 ind_map = row * mInfo->ncol + col;
1253
1254 x = (col + 0.5 - ncol / 2.) * xscale;
1255 y = (row + 0.5 - nrow / 2.) * yscale;
1256
1257 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1258 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1259 by_err[ind_map] = DRMS_MISSING_FLOAT;
1260 xudong 1.1 bz_err[ind_map] = DRMS_MISSING_FLOAT;
1261 continue;
1262 }
1263
1264 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1265 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1266 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1267 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1268 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1269 continue;
1270 }
1271
1272 xi *= rSun; xi = round(xi);
1273 zeta *= rSun; zeta = round(zeta); // nearest neighbor
1274
1275 ind_img = round(zeta * FOURK + xi);
1276
1277 if (errorprop(bT, bA, bI,
1278 errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1279 lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1280 &btSigma2, &bpSigma2, &brSigma2)) {
1281 xudong 1.1 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1282 by_err[ind_map] = DRMS_MISSING_FLOAT;
1283 bz_err[ind_map] = DRMS_MISSING_FLOAT;
1284 continue;
1285 }
1286
1287 bx_err[ind_map] = sqrt(bpSigma2);
1288 by_err[ind_map] = sqrt(btSigma2);
1289 bz_err[ind_map] = sqrt(brSigma2);
1290
1291 }
1292 }
1293
1294 //
1295
1296 free(bT); free(bI); free(bA);
1297 free(errbT); free(errbI); free(errbA);
1298 free(errbTbI); free(errbTbA); free(errbIbA);
1299 return 0;
1300
1301 }
1302 xudong 1.1
1303
1304
1305 /*
1306 * Read full disk vector magnetograms
1307 * Fill factor is 1, use default disambiguity resolution
1308 *
1309 */
1310
1311 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img)
1312 {
1313
1314 int status = 0;
1315
1316 DRMS_Segment_t *inSeg;
1317 DRMS_Array_t *inArray_ambig;
1318 DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
1319
1320 char *ambig;
1321 float *bTotal, *bAzim, *bIncl;
1322
1323 xudong 1.1 inSeg = drms_segment_lookup(inRec, "disambig");
1324 inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
1325 if (status) return 1;
1326 ambig = (char *)inArray_ambig->data;
1327
1328 inSeg = drms_segment_lookup(inRec, "field");
1329 inArray_bTotal = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1330 if (status) return 1;
1331 bTotal = (float *)inArray_bTotal->data;
1332
1333 inSeg = drms_segment_lookup(inRec, "azimuth");
1334 inArray_bAzim = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1335 if (status) return 1;
1336 bAzim = (float *)inArray_bAzim->data;
1337
1338 inSeg = drms_segment_lookup(inRec, "inclination");
1339 inArray_bIncl = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1340 if (status) return 1;
1341 bIncl = (float *)inArray_bIncl->data;
1342
1343 // Convert CCD xyz
1344 xudong 1.1
1345 int llx, lly; // lower-left corner
1346 int bmx, bmy; // bitmap size
1347
1348 llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1349 lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1350
1351 bmx = inArray_ambig->axis[0];
1352 bmy = inArray_ambig->axis[1];
1353
1354 int kx, ky, kOff;
1355 int ix = 0, jy = 0, yOff = 0, iData = 0;
1356 int xDim = FOURK, yDim = FOURK;
1357
1358 for (jy = 0; jy < yDim; jy++)
1359 {
1360 ix = 0;
1361 yOff = jy * xDim;
1362 ky = jy - lly;
1363 for (ix = 0; ix < xDim; ix++)
1364 {
1365 xudong 1.1 iData = yOff + ix;
1366 kx = ix - llx;
1367
1368 // zero azi pointing up, zero incl pointing out from sun
1369 bx_img[iData] = - bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * sin(bAzim[iData] * RADSINDEG);
1370 by_img[iData] = bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * cos(bAzim[iData] * RADSINDEG);
1371 bz_img[iData] = bTotal[iData] * cos(bIncl[iData] * RADSINDEG);
1372
1373 // Disambiguation
1374
1375 if (kx < 0 || kx >= bmx || ky < 0 || ky >= bmy) {
1376 continue;
1377 } else {
1378 kOff = ky * bmx + kx;
1379 if (ambig[kOff] % 2) { // 180
1380 bx_img[iData] *= -1.; by_img[iData] *= -1.;
1381 }
1382 }
1383 }
1384 }
1385
1386 xudong 1.1 // Clean up
1387
1388 drms_free_array(inArray_ambig);
1389 drms_free_array(inArray_bTotal);
1390 drms_free_array(inArray_bAzim);
1391 drms_free_array(inArray_bIncl);
1392
1393 return 0;
1394
1395 }
1396
1397
1398 /*
1399 * Read variances and covariances of vector magnetograms
1400 *
1401 */
1402
1403 int readVectorBErr(DRMS_Record_t *inRec,
1404 float *bT, float *bI, float *bA,
1405 float *errbT, float *errbI, float *errbA,
1406 float *errbTbI, float *errbTbA, float *errbIbA)
1407 xudong 1.1 {
1408
1409 int status = 0;
1410
1411 float *data_ptr[9];
1412 char *segName[9] = {"field", "inclination", "azimuth",
1413 "field_err", "inclination_err", "azimuth_err",
1414 "field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1415 DRMS_Segment_t *inSegs[9];
1416 DRMS_Array_t *inArrays[9];
1417
1418 // Read full disk images
1419
1420 for (int iSeg = 0; iSeg < 9; iSeg++) {
1421
1422 inSegs[iSeg] = drms_segment_lookup(inRec, segName[iSeg]);
1423 inArrays[iSeg] = drms_segment_read(inSegs[iSeg], DRMS_TYPE_FLOAT, &status);
1424 data_ptr[iSeg] = (float *) inArrays[iSeg]->data;
1425
1426 }
1427
1428 xudong 1.1 float *bT0 = data_ptr[0], *bI0 = data_ptr[1], *bA0 = data_ptr[2];
1429 float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
1430 float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
1431
1432 // Convert errors to variances, correlation coefficients to covariances
1433
1434 for (int i = 0; i < FOURK2; i++) {
1435
1436 if (fabs(errbI0[i]) > 180.) errbI0[i] = 180.;
1437 if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
1438
1439 bT[i] = bT0[i];
1440 bI[i] = bI0[i];
1441 bA[i] = bA0[i];
1442
1443 errbT[i] = errbT0[i] * errbT0[i];
1444 errbI[i] = errbI0[i] * errbI0[i] * RADSINDEG * RADSINDEG;
1445 errbA[i] = errbA0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1446
1447 errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1448 errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1449 xudong 1.1 errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1450
1451 }
1452
1453 //
1454
1455 for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1456
1457 return 0;
1458
1459 }
1460
1461
1462 /*
1463 * Create Cutout record: top level subroutine
1464 * Do the loops on segments and set the keywords here
1465 * Work is done in writeCutout routine below
1466 *
1467 */
1468
1469 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1470 xudong 1.1 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1471 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1472 {
1473
1474 int status = 0;
1475
1476 int iHarpSeg;
1477 int nMharpSegs = ARRLENGTH(MharpSegs), nBharpSegs = ARRLENGTH(BharpSegs);
1478
1479 // Cutout Mharp
1480
1481 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
1482 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
1483 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
1484 break;
1485 }
1486 }
|
1689 xudong 1.1
1690 float *bh = (float *) (malloc(nxny * sizeof(float)));
1691 float *bt = (float *) (malloc(nxny * sizeof(float)));
1692 float *jz = (float *) (malloc(nxny * sizeof(float)));
1693 float *bpx = (float *) (malloc(nxny * sizeof(float)));
1694 float *bpy = (float *) (malloc(nxny * sizeof(float)));
1695 float *bpz = (float *) (malloc(nxny * sizeof(float)));
1696 float *derx = (float *) (malloc(nxny * sizeof(float)));
1697 float *dery = (float *) (malloc(nxny * sizeof(float)));
1698 float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1699 float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1700 float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1701 float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1702 float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1703 float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
1704
1705 // Compute
1706
1707 if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1708 mask, cdelt1, rsun_ref, rsun_obs)){
1709 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1710 xudong 1.1 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1711 }
1712
1713 for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1714 greenpot(bpx, bpy, bpz, nx, ny);
1715
1716 computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask);
1717
1718 if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask))
1719 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1720
1721 computeB_total(bx, by, bz, bt, dims, mask);
1722
1723 if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, derx_bt, dery_bt))
1724 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
1725
1726 if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, derx_bh, dery_bh))
1727 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
1728
1729 if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, derx_bz, dery_bz))
1730 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1731 xudong 1.1
1732
1733
1734 if(computeJz(bx, by, dims, jz, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask,
1735 cdelt1, rsun_ref, rsun_obs, derx, dery)) {
1736 swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1737 swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
1738 }
1739
1740 printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1741
1742 if (computeAlpha(bz, dims, jz, &(swKeys_ptr->mean_alpha), mask, cdelt1, rsun_ref, rsun_obs))
1743 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
1744
1745 if (computeHelicity(bz, dims, jz, &(swKeys_ptr->mean_ih),
1746 &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
1747 mask, cdelt1, rsun_ref, rsun_obs)) {
1748 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
1749 swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
1750 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
1751 }
1752 xudong 1.1
1753 if (computeSumAbsPerPolarity(bz, jz, dims, &(swKeys_ptr->totaljz),
1754 mask, cdelt1, rsun_ref, rsun_obs))
1755 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
1756
1757
1758 if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1759 &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1760 mask, cdelt1, rsun_ref, rsun_obs)) {
1761 swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1762 swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
1763 }
1764
1765 if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
1766 &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
1767 &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
1768 mask)) {
1769 swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1770 swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
1771 swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1772 swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
1773 xudong 1.1 }
1774
1775 // Clean up
1776
1777 drms_free_array(maskArray);
1778 drms_free_array(bxArray);
1779 drms_free_array(byArray);
1780 drms_free_array(bzArray);
1781
1782 free(bh); free(bt); free(jz);
1783 free(bpx); free(bpy); free(bpz);
1784 free(derx); free(dery);
1785 free(derx_bt); free(dery_bt);
1786 free(derx_bz); free(dery_bz);
1787 free(derx_bh); free(dery_bh);
1788
1789 }
1790
1791
1792 /*
1793 * Set space weather indices, no error checking for now
1794 xudong 1.1 *
1795 */
1796
1797 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1798 {
1799 drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1800 drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
1801 drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
1802 drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
1803 drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
1804 drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
1805 drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
1806 drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
1807 drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
1808 drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
1809 drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
1810 drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
1811 drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
1812 drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
1813 drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
1814 drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
1815 xudong 1.1 };
1816
1817
1818 /*
1819 * Set all keywords, no error checking for now
1820 *
1821 */
1822
1823 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec)
1824 {
1825 copy_me_keys(inRec, outRec);
1826 copy_patch_keys(inRec, outRec);
1827 copy_geo_keys(inRec, outRec);
1828 copy_ambig_keys(inRec, outRec);
|
Return to
sharp.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps