1 xudong 1.1 /*
2 * sharp.c
3 *
4 * This module creates the pipeline space weather harps
5 * It is a hard-coded strip-down version of bmap.c
6 * It takes the Mharp and Bharp series and crete the following quantities
7 * Series 1: Sharp_CEA
8 * CEA remapped magnetogram, bitmap, continuum, doppler (same size in map coordinate, need manual spec?)
9 * CEA remapped vector field (Br, Bt, Bp) (same as above)
10 * Space weather indices based on vector cutouts (step 2)
11 * Series 2: Sharp_cutout:
12 * cutouts of magnetogram, bitmap, continuum, doppler (HARP defined, various sizes in CCD pixels)
13 * cutouts of all vector data segments (same as above)
14 * Series 3: Other remaps
15 *
16 * Author:
17 * Xudong Sun; Monica Bobra
18 *
19 * Version:
20 * v0.0 Jul 02 2012
21 * v0.1 Jul 23 2012
|
22 xudong 1.7 * v0.2 Sep 04 2012
|
23 xudong 1.6 * v0.3 Dec 18 2012
|
24 xudong 1.7 * v0.4 Jan 02 2013
|
25 xudong 1.10 * v0.5 Jam 23 2012
|
26 xudong 1.1 *
27 * Notes:
28 * v0.0
29 * Mharp & Bharp must be fully specified; other input are series names only
30 * All input records need to match, otherwise quit
31 * Mapping parameters depend on keywords of each record only, not necessarily consistent for now
32 * Cutout doesn't work for char segments yet (drms_segment_readslice bug)
33 * SW indices require ephemeris info which is not passed properly as of now
34 * v0.1
35 * Fixed char I/O thanks to Art
36 * SW indices fixed
37 * Added doppler and continuum
|
38 xudong 1.3 * Added other keywords: HEADER (populated by cvs build version), DATE_B
|
39 xudong 1.6 * v0.3
40 * Fixed memory leakage of 0.15G per rec; denoted with "Dec 18"
|
41 xudong 1.7 * v0.4
42 * Took out convert_inplace(). Was causing all the images to be int
|
43 xudong 1.10 * v0.5
44 * Corrected ephemeris keywords, added argument mInfo for setKeys()
|
45 xudong 1.1 *
46 * Example:
47 * sharp "mharp=hmi.Mharp_720s[1404][2012.02.20_10:00]" \
|
48 xudong 1.6 "bharp=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
49 "dop=hmi.V_720s[2012.02.20_10:00]" \
50 "cont=hmi.Ic_720s[2012.02.20_10:00]" \
51 "sharp_cea=su_xudong.Sharp_CEA" "sharp_cut=su_xudong.Sharp_Cut"
|
52 xudong 1.1 * For comparison:
53 * bmap "in=hmi_test.Bharp_720s_fd10[1404][2012.02.20_10:00]" \
|
54 xudong 1.6 "out=hmi_test.B_720s_CEA" -s -a "map=cyleqa"
|
55 xudong 1.1 *
56 *
57 */
|
58 xudong 1.6
|
59 xudong 1.1 #include <stdio.h>
60 #include <stdlib.h>
61 #include <time.h>
62 #include <sys/time.h>
63 #include <math.h>
64 #include <string.h>
65 #include "jsoc_main.h"
66 #include "astro.h"
67 #include "fstats.h"
68 #include "cartography.c"
69 #include "fresize.h"
70 #include "finterpolate.h"
71 #include "img2helioVector.c"
72 #include "copy_me_keys.c"
73 #include "errorprop.c"
74 #include "sw_functions.c"
75
76 #include <mkl_blas.h>
77 #include <mkl_service.h>
78 #include <mkl_lapack.h>
79 #include <mkl_vml_functions.h>
80 xudong 1.1 #include <omp.h>
|
81 xudong 1.9
|
82 xudong 1.1 #define PI (M_PI)
83 #define RADSINDEG (PI/180.)
84 #define RAD2ARCSEC (648000./M_PI)
85 #define SECINDAY (86400.)
86 #define FOURK (4096)
87 #define FOURK2 (16777216)
88
89 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
90
91 // Nyqvist rate at disk center is 0.03 degree. Oversample above 0.015 degree
92 #define NYQVIST (0.015)
93
94 #ifndef MIN
95 #define MIN(a,b) (((a)<(b)) ? (a) : (b))
96 #endif
97 #ifndef MAX
98 #define MAX(a,b) (((a)>(b)) ? (a) : (b))
99 #endif
100
101 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);}
102 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
103 xudong 1.1
104 #define kNotSpecified "Not Specified"
|
105 xudong 1.6
|
106 xudong 1.1 // Macros for WCS transformations. assume crpix1, crpix2 = CRPIX1, CRPIX2, sina,cosa = sin and cos of CROTA2 resp.
107 // and crvalx and crvaly are CRVAL1 and CRVAL2, cdelt = CDELT1 == CDELT2, then
108 // PIX_X and PIX_Y are CCD pixel addresses, WX and WY are arc-sec W and N on the Sun from disk center.
109 #define PIX_X(wx,wy) ((((wx-crvalx)*cosa + (wy-crvaly)*sina)/cdelt)+crpix1)
110 #define PIX_Y(wx,wy) ((((wy-crvaly)*cosa - (wx-crvalx)*sina)/cdelt)+crpix2)
111 #define WX(pix_x,pix_y) (((pix_x-crpix1)*cosa - (pix_y-crpix2)*sina)*cdelt+crvalx)
112 #define WY(pix_x,pix_y) (((pix_y-crpix2)*cosa + (pix_x-crpix1)*sina)*cdelt+crvaly)
113
114 #define DISAMB_AZI 1
115 #define XSCALE 0.03
116 #define YSCALE 0.03
117 #define NBIN 3
118 #define INTERP 0
119 #define dpath "/home/jsoc/cvs/Development/JSOC"
120
121 /* ========================================================================================================== */
122
123 // Space weather keywords
124 struct swIndex {
125 float mean_vf;
126 float absFlux;
127 xudong 1.1 float mean_hf;
128 float mean_gamma;
129 float mean_derivative_btotal;
130 float mean_derivative_bh;
131 float mean_derivative_bz;
132 float mean_jz;
133 float us_i;
134 float mean_alpha;
135 float mean_ih;
136 float total_us_ih;
137 float total_abs_ih;
138 float totaljz;
139 float totpot;
140 float meanpot;
141 float area_w_shear_gt_45;
142 float meanshear_angle;
143 float area_w_shear_gt_45h;
144 float meanshear_angleh;
145 };
146
147 // Mapping method
148 xudong 1.1 enum projection {
149 carree,
150 cassini,
151 mercator,
152 cyleqa,
153 sineqa,
154 gnomonic,
155 postel,
156 stereographic,
157 orthographic,
158 lambert
159 };
160
|
161 xudong 1.10 // WSC code
162 char *wcsCode[] = {"CAR", "CAS", "MER", "CEA", "GLS", "TAN", "ARC", "STG",
163 "SIN", "ZEA"};
164
|
165 xudong 1.1 // Ephemeris
166 struct ephemeris {
167 double disk_lonc, disk_latc;
168 double disk_xc, disk_yc;
169 double rSun, asd, pa;
170 };
171
172 // Mapping information
173 struct mapInfo {
174 float xc, yc; // reference point: center
175 int nrow, ncol; // size
176 float xscale, yscale; // scale
177 int nbin;
178 enum projection proj; // projection method
179 struct ephemeris ephem; // ephemeris info
180 float *xi_out, *zeta_out; // coordinate on full disk image to sample at
181 };
182
183
184 /* ========================================================================================================== */
185
186 xudong 1.1 /* Get all input data series */
187 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
188 char *mharpQuery, char *bharpQuery);
189
190 /* Check if Mharp and Bharp match */
191 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS);
192
193 /* Get other data series */
194 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS);
195
196 /* Find record from record set with given T_rec */
197 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec);
198
199 // ===================
200
201 /* Create CEA record */
202 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
203 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
204 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
205
206 /* Mapping single segment, wrapper */
207 xudong 1.1 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
208 struct mapInfo *mInfo, char *segName);
209
210 /* Mapping vector magnetogram, wrapper */
211 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
212
213 /* Mapping errors of vector magnetogram, wraper */
214 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo);
215
216 /* Determine reference point coordinate and patch size according to input */
217 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo);
218
219 /* Get ephemeris information */
220 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem);
221
222 /* Compute the coordinates at which the full disk image is sampled */
223 void findCoord(struct mapInfo *mInfo);
224
225 /* Mapping function */
226 int performSampling(float *outData, float *inData, struct mapInfo *mInfo);
227
228 xudong 1.1 /* Performing local vector transformation */
229 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo);
230
231 /* Map and propogate errors */
232 int getBErr(float *bx_err, float *by_err, float *bz_err,
|
233 xudong 1.6 DRMS_Record_t *inRec, struct mapInfo *mInfo);
|
234 xudong 1.1
235 /* Read full disk vector magnetogram */
236 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img);
237
238 /* Read variances and covariances of vector magnetograms */
239 int readVectorBErr(DRMS_Record_t *bharpRec,
240 float *bT, float *bI, float *bA,
241 float *errbT, float *errbI, float *errbA,
242 float *errbTbI, float *errbTbA, float *errbIbA);
243
244 // ===================
245
246 /* Create Cutout record */
247 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
248 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
249 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr);
250
251 /* Get cutout and write segment */
252 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName);
253
254 // ===================
255 xudong 1.1
256 /* Compute space weather indices, no error checking for now */
257 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo);
258
259 /* Set space weather indices, no error checking for now */
260 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr);
261
262 /* Set all keywords, no error checking for now */
|
263 xudong 1.10 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo);
|
264 xudong 1.1
265 // ===================
266
267 /* Nearest neighbor interpolation */
268 float nnb (float *f, int nx, int ny, double x, double y);
269
270 /* Wrapper for Jesper's rebin code */
271 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss);
272
273 /* ========================================================================================================== */
274
275 /* Remap segment names */
276 #define BR_SEG_CEA "Br"
277 #define BT_SEG_CEA "Bt"
278 #define BP_SEG_CEA "Bp"
279 #define BR_ERR_SEG_CEA "Br_err"
280 #define BT_ERR_SEG_CEA "Bt_err"
281 #define BP_ERR_SEG_CEA "Bp_err"
282
283 /* Cutout segment names, input identical to output */
284 char *MharpSegs[] = {"magnetogram", "bitmap"};
285 xudong 1.1 char *BharpSegs[] = {"inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
286 "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
287 "conv_flag", // fixed
288 "info_map", "confid_map",
289 "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
290 "field_inclination_err", "field_az_err", "inclin_azimuth_err",
291 "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
292 "disambig", "conf_disambig"};
293 // For stats
294 char *CutSegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum",
295 "inclination", "azimuth", "field", "vlos_mag", "dop_width", "eta_0",
296 "damping", "src_continuum", "src_grad", "alpha_mag", "chisq",
297 "conv_flag", // fixed
298 "info_map", "confid_map",
299 "inclination_err", "azimuth_err", "field_err", "vlos_err", "alpha_err",
300 "field_inclination_err", "field_az_err", "inclin_azimuth_err",
301 "field_alpha_err","inclination_alpha_err", "azimuth_alpha_err",
302 "disambig", "conf_disambig"};
303 char *CEASegs[] = {"magnetogram", "bitmap", "Dopplergram", "continuum", "disambig",
|
304 xudong 1.6 BR_SEG_CEA, BT_SEG_CEA, BP_SEG_CEA, BR_ERR_SEG_CEA, BT_ERR_SEG_CEA, BP_ERR_SEG_CEA};
|
305 xudong 1.1 /* ========================================================================================================== */
306
307
308
309 char *module_name = "sharp";
|
310 xudong 1.6 char *version_id = "2012 Dec 18"; /* Version number */
|
311 xudong 1.1
312 ModuleArgs_t module_args[] =
313 {
314 {ARG_STRING, "mharp", kNotSpecified, "Input Mharp series."},
315 {ARG_STRING, "bharp", kNotSpecified, "Input Bharp series."},
316 {ARG_STRING, "dop", kNotSpecified, "Input Doppler series."},
317 {ARG_STRING, "cont", kNotSpecified, "Input Continuum series."},
318 {ARG_STRING, "sharp_cea", kNotSpecified, "Output Sharp CEA series."},
319 {ARG_STRING, "sharp_cut", kNotSpecified, "Output Sharp cutout series."},
320 {ARG_END}
321 };
322
323 int DoIt(void)
324 {
325
326 int status = DRMS_SUCCESS;
327 int nrecs, irec;
328
329 char *mharpQuery, *bharpQuery;
330 char *dopQuery, *contQuery;
331 char *sharpCeaQuery, *sharpCutQuery;
332 xudong 1.1
333 DRMS_RecordSet_t *mharpRS = NULL, *bharpRS = NULL;
334 DRMS_RecordSet_t *dopRS = NULL, *contRS = NULL;
335
336 /* Get parameters */
337
338 mharpQuery = (char *) params_get_str(&cmdparams, "mharp");
339 bharpQuery = (char *) params_get_str(&cmdparams, "bharp");
340 dopQuery = (char *) params_get_str(&cmdparams, "dop");
341 contQuery = (char *) params_get_str(&cmdparams, "cont");
342 sharpCeaQuery = (char *) params_get_str(&cmdparams, "sharp_cea");
343 sharpCutQuery = (char *) params_get_str(&cmdparams, "sharp_cut");
344
345 /* Get input data, check everything */
346
347 if (getInputRS(&mharpRS, &bharpRS, mharpQuery, bharpQuery))
348 DIE("Input harp data error.");
349 nrecs = mharpRS->n;
350
351 if (getInputRS_aux(&dopRS, dopQuery, mharpRS))
352 DIE("Input doppler data error.");
353 xudong 1.1
354 if (getInputRS_aux(&contRS, contQuery, mharpRS))
355 DIE("Input continuum data error.");
356
357 /* Start */
358
359 printf("==============\nStart. %d image(s) in total.\n", nrecs);
|
360 xudong 1.6
|
361 xudong 1.1 for (irec = 0; irec < nrecs; irec++) {
362
363 /* Records in work */
364
365 DRMS_Record_t *mharpRec = NULL, *bharpRec = NULL;
366 mharpRec = mharpRS->records[irec];
367 bharpRec = bharpRS->records[irec];
368
369 TIME trec = drms_getkey_time(mharpRec, "T_REC", &status);
|
370 xudong 1.6
|
371 xudong 1.1 struct swIndex swKeys;
372
373 DRMS_Record_t *dopRec = NULL, *contRec = NULL;
374 if (getInputRec_aux(&dopRec, dopRS, trec)) {
375 printf("Fetching Doppler failed, image #%d skipped.\n", irec);
376 continue;
377 }
378 if (getInputRec_aux(&contRec, contRS, trec)) {
379 printf("Fetching continuum failed, image #%d skipped.\n", irec);
380 continue;
381 }
|
382 xudong 1.6
|
383 xudong 1.1 /* Create CEA record */
384
385 DRMS_Record_t *sharpCeaRec = drms_create_record(drms_env, sharpCeaQuery, DRMS_PERMANENT, &status);
386 if (status) { // if failed
387 printf("Creating CEA failed, image #%d skipped.\n", irec);
388 continue;
389 }
390
391 if (createCeaRecord(mharpRec, bharpRec, dopRec, contRec, sharpCeaRec, &swKeys)) { // do the work
392 printf("Creating CEA failed, image #%d skipped.\n", irec);
393 drms_close_record(sharpCeaRec, DRMS_FREE_RECORD);
394 continue;
395 } // swKeys updated here
396
397 drms_close_record(sharpCeaRec, DRMS_INSERT_RECORD);
398
399 /* Create Cutout record */
400
401 DRMS_Record_t *sharpCutRec = drms_create_record(drms_env, sharpCutQuery, DRMS_PERMANENT, &status);
402 if (status) { // if failed
403 printf("Creating cutout failed, image #%d skipped.\n", irec);
404 xudong 1.1 continue;
405 }
406
407 if (createCutRecord(mharpRec, bharpRec, dopRec, contRec, sharpCutRec, &swKeys)) { // do the work
408 printf("Creating cutout failed, image #%d skipped.\n", irec);
409 drms_close_record(sharpCutRec, DRMS_FREE_RECORD);
410 continue;
411 } // swKeys used here
412
413 drms_close_record(sharpCutRec, DRMS_INSERT_RECORD);
414
415 /* Done */
416
417 printf("Image #%d done.\n", irec);
418
419 } // irec
|
420 xudong 1.6
|
421 xudong 1.1
422 drms_close_records(mharpRS, DRMS_FREE_RECORD);
423 drms_close_records(bharpRS, DRMS_FREE_RECORD);
|
424 xudong 1.6 drms_close_records(dopRS, DRMS_FREE_RECORD); // Dec 18 2012
425 drms_close_records(contRS, DRMS_FREE_RECORD); // Dec 18 2012
|
426 xudong 1.1
427 return 0;
|
428 xudong 1.6
|
429 xudong 1.1 } // DoIt
430
431
432 // ===================================================================
433 // ===================================================================
434 // ===================================================================
435
436
437 /*
438 * Get input data series, including mHarp and bharp
439 * Need all records to match, otherwise quit
440 *
441 */
442
443 int getInputRS(DRMS_RecordSet_t **mharpRS_ptr, DRMS_RecordSet_t **bharpRS_ptr,
444 char *mharpQuery, char *bharpQuery)
445 {
446
447 int status = 0;
448
449 *mharpRS_ptr = drms_open_records(drms_env, mharpQuery, &status);
450 xudong 1.1 if (status || (*mharpRS_ptr)->n == 0) return 1;
451
452 *bharpRS_ptr = drms_open_records(drms_env, bharpQuery, &status);
453 if (status || (*bharpRS_ptr)->n == 0) return 1;
|
454 xudong 1.6
|
455 xudong 1.1 if (compareHarp((*mharpRS_ptr), (*bharpRS_ptr))) return 1;
456
457 return 0;
458
459 }
460
461 /*
462 * Check if Mharp and Bharp match
463 *
464 */
465
466 int compareHarp(DRMS_RecordSet_t *mharpRS, DRMS_RecordSet_t *bharpRS)
467 {
468
469 int status = 0;
470 int nrecs = mharpRS->n;
471
472 DRMS_Record_t *mharpRec_t = NULL, *bharpRec_t = NULL; // temporary recs for utility
473
474 if (bharpRS->n != nrecs) {
475 return 1; // return 1 if different
476 xudong 1.1 }
477
478 for (int i = 0; i < nrecs; i++) {
479 mharpRec_t = mharpRS->records[i];
480 bharpRec_t = bharpRS->records[i];
481 if ((drms_getkey_int(mharpRec_t, "HARPNUM", &status) !=
482 drms_getkey_int(bharpRec_t, "HARPNUM", &status)) ||
483 (drms_getkey_time(mharpRec_t, "T_REC", &status) !=
484 drms_getkey_time(bharpRec_t, "T_REC", &status)))
485 {
486 return 1;
487 }
488 }
489
490 return 0;
491
492 }
493
494 /*
495 * Get other data series, check all T_REC are available
496 *
497 xudong 1.1 */
498
499 int getInputRS_aux(DRMS_RecordSet_t **inRS_ptr, char *inQuery, DRMS_RecordSet_t *harpRS)
500 {
501
502 int status = 0;
503
504 *inRS_ptr = drms_open_records(drms_env, inQuery, &status);
505 if (status || (*inRS_ptr)->n == 0) return status;
506
507 // Check if all T_rec are available, need to match both ways
508 int n = harpRS->n, n0 = (*inRS_ptr)->n;
|
509 xudong 1.6
|
510 xudong 1.1 for (int i0 = 0; i0 < n0; i0++) {
511 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
512 TIME trec0 = drms_getkey_time(inRec, "T_REC", &status);
513 TIME trec = 0;
514 for (int i = 0; i < n; i++) {
515 DRMS_Record_t *harpRec = harpRS->records[i];
516 trec = drms_getkey_time(harpRec, "T_REC", &status);
517 if (fabs(trec0 - trec) < 10) break;
518 }
519 if (fabs(trec0 - trec) >= 10) return 1;
520 }
521
522 for (int i = 0; i < n; i++) {
523 DRMS_Record_t *harpRec = harpRS->records[i];
524 TIME trec = drms_getkey_time(harpRec, "T_REC", &status);
525 TIME trec0 = 0;
526 for (int i0 = 0; i0 < n0; i0++) {
527 DRMS_Record_t *inRec = (*inRS_ptr)->records[i0];
528 trec0 = drms_getkey_time(inRec, "T_REC", &status);
529 if (fabs(trec0 - trec) < 10) break;
530 }
531 xudong 1.1 if (fabs(trec0 - trec) >= 10) return 1;
532 }
533
534 return 0;
535
536 }
537
538 /*
539 * Find record from record set with given T_rec
540 *
541 */
542
543 int getInputRec_aux(DRMS_Record_t **inRec_ptr, DRMS_RecordSet_t *inRS, TIME trec)
544 {
545
546 int status = 0;
547
548 int n = inRS->n;
549 for (int i = 0; i < n; i++) {
550 *inRec_ptr = inRS->records[i];
551 TIME trec0 = drms_getkey_time((*inRec_ptr), "T_REC", &status);
552 xudong 1.1 if (fabs(trec0 - trec) < 10) return 0;
553 }
554
555 return 1;
556
557 }
558
559
560
561
562 /*
563 * Create CEA record: top level subroutine
564 * Also compute all the space weather keywords here
565 *
566 */
567
568 int createCeaRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
569 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
570 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
571 {
572
573 xudong 1.1 int status = 0;
574 DRMS_Segment_t *inSeg;
575 DRMS_Array_t *inArray;
576
577 struct mapInfo mInfo;
578 mInfo.proj = (enum projection) cyleqa; // projection method
579 mInfo.xscale = XSCALE;
580 mInfo.yscale = YSCALE;
581 mInfo.nbin = NBIN;
582
583 // Get ephemeris
|
584 xudong 1.6
|
585 xudong 1.3 if (getEphemeris(mharpRec, &(mInfo.ephem))) {
586 SHOW("CEA: get ephemeris error\n");
587 return 1;
588 }
|
589 xudong 1.1
590 // Find position
591
|
592 xudong 1.3 if (findPosition(mharpRec, &mInfo)) {
593 SHOW("CEA: find position error\n");
594 return 1;
595 }
|
596 xudong 1.1
597 // Create xi_out, zeta_out array in mInfo:
598 // Coordinates to sample in original full disk image
599
600 int ncol0, nrow0; // oversampled map size
601 ncol0 = mInfo.ncol * mInfo.nbin + (mInfo.nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
602 nrow0 = mInfo.nrow * mInfo.nbin + (mInfo.nbin / 2) * 2;
603
604 mInfo.xi_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
605 mInfo.zeta_out = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
606
607 findCoord(&mInfo); // compute it here so it could be shared by the following 4 functions
608
609 // Mapping single segment: Mharp, etc.
|
610 xudong 1.6
|
611 xudong 1.3 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "magnetogram")) {
612 SHOW("CEA: mapping magnetogram error\n");
613 return 1;
614 }
615 if (mapScaler(sharpRec, mharpRec, mharpRec, &mInfo, "bitmap")) {
616 SHOW("CEA: mapping bitmap error\n");
617 return 1;
618 }
|
619 xudong 1.1 printf("Magnetogram mapping done.\n");
620
|
621 xudong 1.3 if (mapScaler(sharpRec, dopRec, mharpRec, &mInfo, "Dopplergram")) {
622 SHOW("CEA: mapping dopplergram error\n");
623 return 1;
624 }
|
625 xudong 1.1 printf("Dopplergram mapping done.\n");
626
|
627 xudong 1.3 if (mapScaler(sharpRec, contRec, mharpRec, &mInfo, "continuum")) {
628 SHOW("CEA: mapping continuum error\n");
629 return 1;
630 }
|
631 xudong 1.1 printf("Intensitygram mapping done.\n");
|
632 xudong 1.6
|
633 xudong 1.3 if (mapScaler(sharpRec, bharpRec, mharpRec, &mInfo, "conf_disambig")) {
634 SHOW("CEA: mapping conf_disambig error\n");
635 return 1;
636 }
|
637 xudong 1.2 printf("Conf disambig mapping done.\n");
|
638 xudong 1.1
639 // Mapping vector B
640
|
641 xudong 1.3 if (mapVectorB(sharpRec, bharpRec, &mInfo)) {
642 SHOW("CEA: mapping vector B error\n");
643 return 1;
644 }
|
645 xudong 1.1 printf("Vector B mapping done.\n");
646
647 // Mapping vector B errors
648
|
649 xudong 1.3 if (mapVectorBErr(sharpRec, bharpRec, &mInfo)) {
650 SHOW("CEA: mapping vector B uncertainty error\n");
651 return 1;
652 }
|
653 xudong 1.1 printf("Vector B error done.\n");
654
655 // Keywords & Links
656
657 drms_copykey(sharpRec, mharpRec, "T_REC");
658 drms_copykey(sharpRec, mharpRec, "HARPNUM");
659
660 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
661 if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
662 DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
663 if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
664
|
665 xudong 1.10 setKeys(sharpRec, bharpRec, &mInfo); // Set all other keywords
|
666 xudong 1.3 drms_copykey(sharpRec, mharpRec, "QUALITY"); // copied from los records
|
667 xudong 1.6
|
668 xudong 1.1 // Space weather
669
670 computeSWIndex(swKeys_ptr, sharpRec, &mInfo); // compute it!
671 printf("Space weather indices done.\n");
672
673 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
|
674 xudong 1.6
|
675 xudong 1.1 // Stats
676
677 int nCEASegs = ARRLENGTH(CEASegs);
678 for (int iSeg = 0; iSeg < nCEASegs; iSeg++) {
679 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
680 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
681 int stat = set_statistics(outSeg, outArray, 1);
|
682 xudong 1.6 // printf("%d => %d\n", iSeg, stat);
|
683 xudong 1.1 drms_free_array(outArray);
684 }
685
686 free(mInfo.xi_out);
687 free(mInfo.zeta_out);
688 return 0;
689
690 }
691
692
693 /*
694 * Mapping a single segment
695 * Read in full disk image, utilize mapImage for mapping
696 * then write the segment out, segName same in in/out Rec
697 *
698 */
699
700 int mapScaler(DRMS_Record_t *sharpRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec,
701 struct mapInfo *mInfo, char *segName)
702 {
703
704 xudong 1.1 int status = 0;
705 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
706 int dims[2] = {nx, ny};
707
708 // Input full disk array
709
710 DRMS_Segment_t *inSeg = NULL;
711 inSeg = drms_segment_lookup(inRec, segName);
712 if (!inSeg) return 1;
713
714 DRMS_Array_t *inArray = NULL;
715 inArray = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
716 if (!inArray) return 1;
717
718 float *inData;
719 int xsz = inArray->axis[0], ysz = inArray->axis[1];
720 if ((xsz != FOURK) || (ysz != FOURK)) { // for bitmap, make tmp full disk
721 float *inData0 = (float *) inArray->data;
722 inData = (float *) (calloc(FOURK2, sizeof(float)));
723 int x0 = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1;
724 int y0 = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1;
725 xudong 1.1 int ind_map;
726 for (int row = 0; row < ysz; row++) {
727 for (int col = 0; col < xsz; col++) {
728 ind_map = (row + y0) * FOURK + (col + x0);
729 inData[ind_map] = inData0[row * xsz + col];
730 }
731 }
732 drms_free_array(inArray); inArray = NULL;
733 } else {
734 inData = (float *) inArray->data;
735 }
736
737 // Mapping
738
739 float *map = (float *) (malloc(nxny * sizeof(float)));
740 if (performSampling(map, inData, mInfo))
|
741 xudong 1.6 {if (inArray) drms_free_array(inArray); free(map); return 1;}
|
742 xudong 1.1
743 // Write out
744
745 DRMS_Segment_t *outSeg = NULL;
746 outSeg = drms_segment_lookup(sharpRec, segName);
747 if (!outSeg) return 1;
748
|
749 xudong 1.7 // DRMS_Type_t arrayType = outSeg->info->type;
|
750 xudong 1.1 DRMS_Array_t *outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, map, &status);
751 if (status) {if (inArray) drms_free_array(inArray); free(map); return 1;}
752
753 // convert to needed data type
754
|
755 xudong 1.7 // drms_array_convert_inplace(outSeg->info->type, 0, 1, outArray); // Jan 02 2013
|
756 xudong 1.1
757 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
|
758 xudong 1.8 // outArray->parent_segment = outSeg;
759 outArray->israw = 0; // always compressed
|
760 xudong 1.1 outArray->bzero = outSeg->bzero;
761 outArray->bscale = outSeg->bscale;
|
762 xudong 1.6
|
763 xudong 1.1 status = drms_segment_write(outSeg, outArray, 0);
764 if (status) return 0;
765
766 if (inArray) drms_free_array(inArray);
|
767 xudong 1.6 if ((xsz != FOURK) || (ysz != FOURK)) free(inData); // Dec 18 2012
|
768 xudong 1.1 if (outArray) drms_free_array(outArray);
769 return 0;
770
771 }
772
773
774 /*
775 * Mapping vector magnetogram
776 *
777 */
778
779 int mapVectorB(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
780 {
781
782 int status = 0;
783 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
784 int dims[2] = {nx, ny};
785
786 // Read in segments, filling factor assume to be 1
787
788 float *bx_img = (float *) (malloc(FOURK2 * sizeof(float)));
789 xudong 1.1 float *by_img = (float *) (malloc(FOURK2 * sizeof(float)));
790 float *bz_img = (float *) (malloc(FOURK2 * sizeof(float)));
791
792 if (readVectorB(bharpRec, bx_img, by_img, bz_img)) {
793 printf("Read full disk image error\n");
794 free(bx_img); free(by_img); free(bz_img);
795 return 1;
796 }
797
798 // Mapping
799
800 float *bx_map = NULL, *by_map = NULL, *bz_map = NULL; // intermediate maps, in CCD bxyz representation
|
801 xudong 1.6
|
802 xudong 1.1 bx_map = (float *) (malloc(nxny * sizeof(float)));
803 if (performSampling(bx_map, bx_img, mInfo))
|
804 xudong 1.6 {free(bx_img); free(by_img); free(bz_img); free(bx_map); return 1;}
805
|
806 xudong 1.1 by_map = (float *) (malloc(nxny * sizeof(float)));
807 if (performSampling(by_map, by_img, mInfo))
|
808 xudong 1.6 {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
809
|
810 xudong 1.1 bz_map = (float *) (malloc(nxny * sizeof(float)));
811 if (performSampling(bz_map, bz_img, mInfo))
|
812 xudong 1.6 {free(bx_img); free(by_img); free(bz_img); free(bz_map); return 1;}
|
813 xudong 1.1
814 free(bx_img); free(by_img); free(bz_img);
815
816 // Vector transform
817
818 vectorTransform(bx_map, by_map, bz_map, mInfo);
819
820 for (int i = 0; i < nxny; i++) by_map[i] *= -1; // positive theta pointing south
821
822 // Write out
823
824 DRMS_Segment_t *outSeg;
825 DRMS_Array_t *outArray;
826
827 float *data_prt[3] = {bx_map, by_map, bz_map};
828 char *segName[3] = {BP_SEG_CEA, BT_SEG_CEA, BR_SEG_CEA};
829
830 for (int iSeg = 0; iSeg < 3; iSeg++) {
831 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
832 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
833 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
|
834 xudong 1.8 // outArray->parent_segment = outSeg;
835 outArray->israw = 0;
|
836 xudong 1.1 outArray->bzero = outSeg->bzero;
837 outArray->bscale = outSeg->bscale;
838 status = drms_segment_write(outSeg, outArray, 0);
839 if (status) return 1;
840 drms_free_array(outArray);
841 }
842
843 //
844
845 return 0;
846
847 }
848
849
850 /*
851 * Mapping vector magnetogram errors
852 *
853 */
854
855 int mapVectorBErr(DRMS_Record_t *sharpRec, DRMS_Record_t *bharpRec, struct mapInfo *mInfo)
856 {
857 xudong 1.1
858 int status = 0;
859
860 int nx = mInfo->ncol, ny = mInfo->nrow, nxny = nx * ny;
861 int dims[2] = {nx, ny};
862
863 // Compute propogated errors, using nearest neighbour interpolation
864
865 float *bx_err = (float *) (malloc(nxny * sizeof(float)));
866 float *by_err = (float *) (malloc(nxny * sizeof(float)));
867 float *bz_err = (float *) (malloc(nxny * sizeof(float)));
868
869 if (getBErr(bx_err, by_err, bz_err, bharpRec, mInfo)) {
870 free(bx_err); free(by_err); free(bz_err);
871 return 1;
872 }
|
873 xudong 1.6
|
874 xudong 1.1 // Write out
875
876 DRMS_Segment_t *outSeg;
877 DRMS_Array_t *outArray;
878
879 float *data_prt[3] = {bx_err, by_err, bz_err};
880 char *segName[3] = {BP_ERR_SEG_CEA, BT_ERR_SEG_CEA, BR_ERR_SEG_CEA};
881
882 for (int iSeg = 0; iSeg < 3; iSeg++) {
883 outSeg = drms_segment_lookup(sharpRec, segName[iSeg]);
884 outArray = drms_array_create(DRMS_TYPE_FLOAT, 2, dims, data_prt[iSeg], &status);
885 outSeg->axis[0] = outArray->axis[0]; outSeg->axis[1] = outArray->axis[1];
|
886 xudong 1.8 // outArray->parent_segment = outSeg;
887 outArray->israw = 0;
|
888 xudong 1.1 outArray->bzero = outSeg->bzero;
889 outArray->bscale = outSeg->bscale;
890 status = drms_segment_write(outSeg, outArray, 0);
891 if (status) return 1;
892 drms_free_array(outArray);
893 }
894
895 //
896
897 return 0;
898
899 }
900
901
902
903 /*
904 * Determine reference point coordinate and patch size according to keywords
905 * xc, yc are the coordinate of patch center, in degrees
906 * ncol and nrow are the final size
907 *
908 */
909 xudong 1.1
910 int findPosition(DRMS_Record_t *inRec, struct mapInfo *mInfo)
911 {
912
913 int status = 0;
914 int harpnum = drms_getkey_int(inRec, "HARPNUM", &status);
915 TIME trec = drms_getkey_time(inRec, "T_REC", &status);
916 float disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status);
917
918 /* Center coord */
919
920 float minlon = drms_getkey_float(inRec, "LONDTMIN", &status); if (status) return 1; // Stonyhurst lon
921 float maxlon = drms_getkey_float(inRec, "LONDTMAX", &status); if (status) return 1;
922 float minlat = drms_getkey_float(inRec, "LATDTMIN", &status); if (status) return 1;
923 float maxlat = drms_getkey_float(inRec, "LATDTMAX", &status); if (status) return 1;
924
925 // A bug fixer for HARP (per M. Turmon)
926 // When AR is below threshold, "LONDTMIN", "LONDTMAX" will be wrong
927 // Also keywords such as "SIZE" will be NaN
928 // We compute minlon & minlat then by
929 // LONDTMIN(t) = LONDTMIN(t0) + (t - t0) * OMEGA_DT
930 xudong 1.1
|
931 xudong 1.13 // float psize = drms_getkey_float(inRec, "SIZE", &status);
932 // if (psize != psize) {
933
934 if (minlon != minlon || maxlon != maxlon) { // check lons instead of SIZE
|
935 xudong 1.12 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
|
936 xudong 1.1 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1;
937 char firstRecQuery[100], t0_str[100];
938 sprint_time(t0_str, t0, "TAI", 0);
939 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str);
940 DRMS_RecordSet_t *tmpRS = drms_open_records(drms_env, firstRecQuery, &status);
941 if (status || tmpRS->n != 1) return 1;
942 DRMS_Record_t *tmpRec = tmpRS->records[0];
943 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1;
944 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1;
945 minlon = minlon0 + (trec - t0) * omega / SECINDAY;
946 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY;
947 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon);
948 }
949
950 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc;
951 mInfo->yc = (maxlat + minlat) / 2.;
952
953 /* Size */
954
955 mInfo->ncol = round((maxlon - minlon) / mInfo->xscale);
956 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale);
957 xudong 1.1
958 return 0;
959
960 }
961
962
963 /*
964 * Fetch ephemeris info from a DRMS record
965 * No error checking for now
966 *
967 */
968
969 int getEphemeris(DRMS_Record_t *inRec, struct ephemeris *ephem)
970 {
971
972 int status = 0;
973
974 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation
975 double sina = sin(crota2 * RADSINDEG);
976 double cosa = cos(crota2 * RADSINDEG);
977
978 xudong 1.1 ephem->pa = - crota2 * RADSINDEG;
979 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG;
980 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG;
981
982 float crvalx = 0.0;
983 float crvaly = 0.0;
984 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
985 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
986 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
|
987 xudong 1.6 printf("cdelt=%f\n",cdelt);
|
988 xudong 1.1 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0
989 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0;
990
991 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status);
992 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status);
993 if (status) rSun_ref = 6.96e8;
994
995 ephem->asd = asin(rSun_ref/dSun);
996 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt;
997
998 return 0;
999
1000 }
1001
1002
1003 /*
1004 * Compute the coordinates to be sampled on full disk image
1005 * mInfo->xi_out & mInfo->zeta_out
1006 * This is oversampled, its size is ncol0 & nrow0 as shown below
1007 *
1008 *
1009 xudong 1.1 */
1010
1011 void findCoord(struct mapInfo *mInfo)
1012 {
1013
1014 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
1015 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1016
1017 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution
1018 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad
1019
1020 double lonc = mInfo->xc * RADSINDEG; // in rad
1021 double latc = mInfo->yc * RADSINDEG;
1022
1023 double disk_lonc = (mInfo->ephem).disk_lonc;
1024 double disk_latc = (mInfo->ephem).disk_latc;
1025
1026 double rSun = (mInfo->ephem).rSun;
1027 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1028 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1029 double pa = (mInfo->ephem).pa;
1030 xudong 1.1
1031 // Temp pointers
1032
1033 float *xi_out = mInfo->xi_out;
1034 float *zeta_out = mInfo->zeta_out;
1035
1036 // start
1037
1038 double x, y; // map coord
1039 double lat, lon; // helio coord
1040 double xi, zeta; // image coord (for one point)
1041
1042 int ind_map;
1043
1044 for (int row0 = 0; row0 < nrow0; row0++) {
1045 for (int col0 = 0; col0 < ncol0; col0++) {
1046
1047 ind_map = row0 * ncol0 + col0;
1048
1049 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad
1050 y = (row0 + 0.5 - nrow0/2.) * yscale0;
1051 xudong 1.1
1052 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates.
1053 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap].
1054 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc]
1055 * is the heliographic longitude and latitude of the map center. Both are in degree.
1056 */
1057
1058 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1059 xi_out[ind_map] = -1;
1060 zeta_out[ind_map] = -1;
1061 continue;
1062 }
1063
1064 /* map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the
1065 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given.
1066 */
1067
1068 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1069 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1070 xi_out[ind_map] = -1;
1071 zeta_out[ind_map] = -1;
1072 xudong 1.1 continue;
1073 }
1074
1075 xi_out[ind_map] = xi * rSun;
1076 zeta_out[ind_map] = zeta * rSun;
1077
1078 }
1079 }
1080
1081 }
1082
1083
1084 /*
1085 * Sampling function
1086 * oversampling by nbin, then binning using a Gaussian
1087 * save results in outData, always of float type
1088 *
1089 */
1090
1091 int performSampling(float *outData, float *inData, struct mapInfo *mInfo)
1092 {
1093 xudong 1.1
1094 int status = 0;
1095
1096 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN
1097 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2;
1098
1099 float *outData0 = (float *) (malloc(ncol0 * nrow0 * sizeof(float)));
1100
1101 float *xi_out = mInfo->xi_out;
1102 float *zeta_out = mInfo->zeta_out;
|
1103 xudong 1.6
|
1104 xudong 1.1 // Interpolation
1105
1106 struct fint_struct pars;
1107 int interpOpt = INTERP; // Use Wiener by default, 6 order, 1 constraint
1108
1109 switch (interpOpt) {
1110 case 0: // Wiener, 6 order, 1 constraint
1111 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath);
1112 break;
1113 case 1: // Cubic convolution
1114 init_finterpolate_cubic_conv(&pars, 1., 3.);
1115 break;
1116 case 2: // Bilinear
1117 init_finterpolate_linear(&pars, 1.);
1118 break;
1119 default:
1120 return 1;
1121 }
1122
1123 finterpolate(&pars, inData, xi_out, zeta_out, outData0,
1124 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT);
1125 xudong 1.1
1126 // Rebinning, smoothing
1127
1128 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
|
1129 xudong 1.6 free(outData0); // Dec 18 2012
|
1130 xudong 1.1
1131 //
1132
1133 return 0;
1134
1135 }
1136
1137
1138 /*
1139 * Performing local vector transformation
1140 * xyz: z refers to vertical (radial) component, x EW (phi), y NS
1141 *
1142 */
1143
1144 void vectorTransform(float *bx_map, float *by_map, float *bz_map, struct mapInfo *mInfo)
1145 {
1146
1147 int ncol = mInfo->ncol;
1148 int nrow = mInfo->nrow;
1149
1150 float xscale = mInfo->xscale * RADSINDEG; // in rad
1151 xudong 1.1 float yscale = mInfo->yscale * RADSINDEG;
1152
1153 double lonc = mInfo->xc * RADSINDEG; // in rad
1154 double latc = mInfo->yc * RADSINDEG;
1155
1156 double disk_lonc = (mInfo->ephem).disk_lonc;
1157 double disk_latc = (mInfo->ephem).disk_latc;
1158
1159 double rSun = (mInfo->ephem).rSun;
1160 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1161 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1162 double pa = (mInfo->ephem).pa;
1163
1164 int ind_map;
1165 double x, y;
1166 double lat, lon; // lat / lon for current point
1167
1168 double bx_tmp, by_tmp, bz_tmp;
1169
1170 //
1171
1172 xudong 1.1 for (int row = 0; row < mInfo->nrow; row++) {
1173 for (int col = 0; col < mInfo->ncol; col++) {
1174
1175 ind_map = row * mInfo->ncol + col;
1176
1177 x = (col + 0.5 - ncol / 2.) * xscale;
1178 y = (row + 0.5 - nrow / 2.) * yscale;
1179
1180 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1181 bx_map[ind_map] = DRMS_MISSING_FLOAT;
1182 by_map[ind_map] = DRMS_MISSING_FLOAT;
1183 bz_map[ind_map] = DRMS_MISSING_FLOAT;
1184 continue;
1185 }
1186
1187 bx_tmp = by_tmp = bz_tmp = 0;
1188
1189 img2helioVector (bx_map[ind_map], by_map[ind_map], bz_map[ind_map],
1190 &bx_tmp, &by_tmp, &bz_tmp,
1191 lon, lat, disk_lonc, disk_latc, pa);
1192
1193 xudong 1.1 bx_map[ind_map] = bx_tmp;
1194 by_map[ind_map] = by_tmp;
1195 bz_map[ind_map] = bz_tmp;
1196
1197 }
1198 }
|
1199 xudong 1.6
|
1200 xudong 1.1 }
1201
1202
1203
1204 /*
1205 * Map and propogate vector field errors
1206 *
1207 */
1208
1209 int getBErr(float *bx_err, float *by_err, float *bz_err,
|
1210 xudong 1.6 DRMS_Record_t *inRec, struct mapInfo *mInfo)
|
1211 xudong 1.1 {
1212
1213 int status = 0;
1214
1215 // Get variances and covariances, filling factor assume to be 1
1216
1217 float *bT = (float *) (malloc(FOURK2 * sizeof(float))); // field
1218 float *bI = (float *) (malloc(FOURK2 * sizeof(float))); // inclination
1219 float *bA = (float *) (malloc(FOURK2 * sizeof(float))); // azimuth
1220
1221 float *errbT = (float *) (malloc(FOURK2 * sizeof(float)));
1222 float *errbI = (float *) (malloc(FOURK2 * sizeof(float)));
1223 float *errbA = (float *) (malloc(FOURK2 * sizeof(float)));
1224
1225 float *errbTbI = (float *) (malloc(FOURK2 * sizeof(float)));
1226 float *errbTbA = (float *) (malloc(FOURK2 * sizeof(float)));
1227 float *errbIbA = (float *) (malloc(FOURK2 * sizeof(float)));
1228
1229 if (readVectorBErr(inRec,
1230 bT, bI, bA,
1231 errbT, errbI, errbA,
1232 xudong 1.1 errbTbI, errbTbA, errbIbA)) {
1233 printf("Read full disk variances & covariances error\n");
1234 free(bT); free(bI); free(bA);
1235 free(errbT); free(errbI); free(errbA);
1236 free(errbTbI); free(errbTbA); free(errbIbA);
1237 return 1;
1238 }
1239
1240 // Size
1241
1242 int ncol = mInfo->ncol;
1243 int nrow = mInfo->nrow;
1244
1245 float xscale = mInfo->xscale * RADSINDEG; // in rad
1246 float yscale = mInfo->yscale * RADSINDEG;
1247
1248 double lonc = mInfo->xc * RADSINDEG; // in rad
1249 double latc = mInfo->yc * RADSINDEG;
1250
1251 double disk_lonc = (mInfo->ephem).disk_lonc;
1252 double disk_latc = (mInfo->ephem).disk_latc;
1253 xudong 1.1
1254 double rSun = (mInfo->ephem).rSun;
1255 double disk_xc = (mInfo->ephem).disk_xc / rSun;
1256 double disk_yc = (mInfo->ephem).disk_yc / rSun;
1257 double pa = (mInfo->ephem).pa;
1258
1259 // Start
1260
1261 double x, y; // map coord
1262 double lat, lon; // spherical coord
1263 double xi, zeta; // image coord, round to full pixel
1264
1265 int ind_map, ind_img;
1266
1267 double bpSigma2, btSigma2, brSigma2; // variances after prop
|
1268 xudong 1.6
|
1269 xudong 1.1 for (int row = 0; row < mInfo->nrow; row++) {
1270 for (int col = 0; col < mInfo->ncol; col++) {
1271
1272 ind_map = row * mInfo->ncol + col;
1273
1274 x = (col + 0.5 - ncol / 2.) * xscale;
1275 y = (row + 0.5 - nrow / 2.) * yscale;
1276
1277 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) {
1278 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1279 by_err[ind_map] = DRMS_MISSING_FLOAT;
1280 bz_err[ind_map] = DRMS_MISSING_FLOAT;
1281 continue;
1282 }
1283
1284 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta,
1285 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) {
1286 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1287 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1288 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1289 continue;
1290 xudong 1.1 }
1291
1292 xi *= rSun; xi = round(xi);
1293 zeta *= rSun; zeta = round(zeta); // nearest neighbor
1294
1295 ind_img = round(zeta * FOURK + xi);
1296
1297 if (errorprop(bT, bA, bI,
1298 errbT, errbA, errbI, errbTbA, errbTbI, errbIbA,
1299 lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta,
1300 &btSigma2, &bpSigma2, &brSigma2)) {
1301 bx_err[ind_map] = DRMS_MISSING_FLOAT;
1302 by_err[ind_map] = DRMS_MISSING_FLOAT;
1303 bz_err[ind_map] = DRMS_MISSING_FLOAT;
1304 continue;
1305 }
1306
1307 bx_err[ind_map] = sqrt(bpSigma2);
1308 by_err[ind_map] = sqrt(btSigma2);
1309 bz_err[ind_map] = sqrt(brSigma2);
1310
1311 xudong 1.1 }
1312 }
1313
1314 //
1315
1316 free(bT); free(bI); free(bA);
1317 free(errbT); free(errbI); free(errbA);
1318 free(errbTbI); free(errbTbA); free(errbIbA);
1319 return 0;
1320
1321 }
1322
1323
1324
1325 /*
1326 * Read full disk vector magnetograms
1327 * Fill factor is 1, use default disambiguity resolution
1328 *
1329 */
1330
1331 int readVectorB(DRMS_Record_t *inRec, float *bx_img, float *by_img, float *bz_img)
1332 xudong 1.1 {
1333
1334 int status = 0;
1335
1336 DRMS_Segment_t *inSeg;
1337 DRMS_Array_t *inArray_ambig;
|
1338 xudong 1.6 DRMS_Array_t *inArray_bTotal, *inArray_bAzim, *inArray_bIncl;
|
1339 xudong 1.1
1340 char *ambig;
1341 float *bTotal, *bAzim, *bIncl;
1342
1343 inSeg = drms_segment_lookup(inRec, "disambig");
1344 inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status);
1345 if (status) return 1;
1346 ambig = (char *)inArray_ambig->data;
1347
1348 inSeg = drms_segment_lookup(inRec, "field");
1349 inArray_bTotal = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1350 if (status) return 1;
1351 bTotal = (float *)inArray_bTotal->data;
1352
1353 inSeg = drms_segment_lookup(inRec, "azimuth");
1354 inArray_bAzim = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1355 if (status) return 1;
1356 bAzim = (float *)inArray_bAzim->data;
1357
1358 inSeg = drms_segment_lookup(inRec, "inclination");
1359 inArray_bIncl = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status);
1360 xudong 1.1 if (status) return 1;
1361 bIncl = (float *)inArray_bIncl->data;
1362
1363 // Convert CCD xyz
1364
1365 int llx, lly; // lower-left corner
1366 int bmx, bmy; // bitmap size
1367
1368 llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1;
1369 lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1;
1370
1371 bmx = inArray_ambig->axis[0];
1372 bmy = inArray_ambig->axis[1];
1373
1374 int kx, ky, kOff;
1375 int ix = 0, jy = 0, yOff = 0, iData = 0;
1376 int xDim = FOURK, yDim = FOURK;
1377
1378 for (jy = 0; jy < yDim; jy++)
1379 {
1380 ix = 0;
1381 xudong 1.1 yOff = jy * xDim;
1382 ky = jy - lly;
1383 for (ix = 0; ix < xDim; ix++)
1384 {
1385 iData = yOff + ix;
1386 kx = ix - llx;
1387
1388 // zero azi pointing up, zero incl pointing out from sun
1389 bx_img[iData] = - bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * sin(bAzim[iData] * RADSINDEG);
1390 by_img[iData] = bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * cos(bAzim[iData] * RADSINDEG);
1391 bz_img[iData] = bTotal[iData] * cos(bIncl[iData] * RADSINDEG);
1392
1393 // Disambiguation
1394
1395 if (kx < 0 || kx >= bmx || ky < 0 || ky >= bmy) {
1396 continue;
1397 } else {
1398 kOff = ky * bmx + kx;
1399 if (ambig[kOff] % 2) { // 180
1400 bx_img[iData] *= -1.; by_img[iData] *= -1.;
1401 }
1402 xudong 1.1 }
1403 }
1404 }
1405
1406 // Clean up
1407
1408 drms_free_array(inArray_ambig);
1409 drms_free_array(inArray_bTotal);
1410 drms_free_array(inArray_bAzim);
1411 drms_free_array(inArray_bIncl);
1412
1413 return 0;
1414
1415 }
1416
1417
1418 /*
1419 * Read variances and covariances of vector magnetograms
1420 *
1421 */
1422
1423 xudong 1.1 int readVectorBErr(DRMS_Record_t *inRec,
1424 float *bT, float *bI, float *bA,
1425 float *errbT, float *errbI, float *errbA,
1426 float *errbTbI, float *errbTbA, float *errbIbA)
1427 {
1428
1429 int status = 0;
1430
1431 float *data_ptr[9];
1432 char *segName[9] = {"field", "inclination", "azimuth",
|
1433 xudong 1.6 "field_err", "inclination_err", "azimuth_err",
1434 "field_inclination_err", "field_az_err", "inclin_azimuth_err"};
|
1435 xudong 1.1 DRMS_Segment_t *inSegs[9];
1436 DRMS_Array_t *inArrays[9];
1437
1438 // Read full disk images
1439
1440 for (int iSeg = 0; iSeg < 9; iSeg++) {
1441
1442 inSegs[iSeg] = drms_segment_lookup(inRec, segName[iSeg]);
1443 inArrays[iSeg] = drms_segment_read(inSegs[iSeg], DRMS_TYPE_FLOAT, &status);
1444 data_ptr[iSeg] = (float *) inArrays[iSeg]->data;
1445
1446 }
1447
1448 float *bT0 = data_ptr[0], *bI0 = data_ptr[1], *bA0 = data_ptr[2];
1449 float *errbT0 = data_ptr[3], *errbI0 = data_ptr[4], *errbA0 = data_ptr[5];
1450 float *errbTbI0 = data_ptr[6], *errbTbA0 = data_ptr[7], *errbIbA0 = data_ptr[8];
1451
1452 // Convert errors to variances, correlation coefficients to covariances
1453
1454 for (int i = 0; i < FOURK2; i++) {
1455
1456 xudong 1.1 if (fabs(errbI0[i]) > 180.) errbI0[i] = 180.;
1457 if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.;
1458
1459 bT[i] = bT0[i];
1460 bI[i] = bI0[i];
1461 bA[i] = bA0[i];
1462
1463 errbT[i] = errbT0[i] * errbT0[i];
1464 errbI[i] = errbI0[i] * errbI0[i] * RADSINDEG * RADSINDEG;
1465 errbA[i] = errbA0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1466
1467 errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG;
1468 errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG;
1469 errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
|
1470 xudong 1.6
|
1471 xudong 1.1 }
1472
1473 //
1474
1475 for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
|
1476 xudong 1.6
|
1477 xudong 1.1 return 0;
1478
1479 }
1480
1481
1482 /*
1483 * Create Cutout record: top level subroutine
1484 * Do the loops on segments and set the keywords here
1485 * Work is done in writeCutout routine below
1486 *
1487 */
1488
1489 int createCutRecord(DRMS_Record_t *mharpRec, DRMS_Record_t *bharpRec,
1490 DRMS_Record_t *dopRec, DRMS_Record_t *contRec,
1491 DRMS_Record_t *sharpRec, struct swIndex *swKeys_ptr)
1492 {
1493
1494 int status = 0;
1495
1496 int iHarpSeg;
1497 int nMharpSegs = ARRLENGTH(MharpSegs), nBharpSegs = ARRLENGTH(BharpSegs);
1498 xudong 1.1
1499 // Cutout Mharp
1500
1501 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) {
1502 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) {
1503 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]);
1504 break;
1505 }
1506 }
|
1507 xudong 1.3 if (iHarpSeg != nMharpSegs) {
1508 SHOW("Cutout: segment number unmatch\n");
1509 return 1; // if failed
1510 }
|
1511 xudong 1.1 printf("Magnetogram cutout done.\n");
1512
1513 // Cutout Doppler
1514
1515 if (writeCutout(sharpRec, dopRec, mharpRec, "Dopplergram")) {
1516 printf("Doppler cutout failed\n");
1517 return 1;
1518 }
1519 printf("Dopplergram cutout done.\n");
1520
1521 // Cutout Continuum
1522
1523 if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) {
1524 printf("Continuum cutout failed\n");
1525 return 1;
1526 }
1527 printf("Intensitygram cutout done.\n");
1528
1529 // Coutout Bharp
1530
1531 for (iHarpSeg = 0; iHarpSeg < nBharpSegs; iHarpSeg++) {
1532 xudong 1.1 if (writeCutout(sharpRec, bharpRec, mharpRec, BharpSegs[iHarpSeg])) {
1533 printf("Bharp cutout fails for %s\n", BharpSegs[iHarpSeg]);
1534 break;
1535 }
1536 }
1537 if (iHarpSeg != nBharpSegs) return 1; // if failed
1538 printf("Vector B cutout done.\n");
1539
1540 // Keywords & Links
1541
1542 drms_copykey(sharpRec, mharpRec, "T_REC");
1543 drms_copykey(sharpRec, mharpRec, "HARPNUM");
1544
1545 DRMS_Link_t *mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP");
1546 if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec);
1547 DRMS_Link_t *bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP");
1548 if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec);
1549
1550 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
|
1551 xudong 1.10 setKeys(sharpRec, bharpRec, NULL); // Set all other keywords, NULL specifies cutout
|
1552 xudong 1.6
|
1553 xudong 1.1 // Stats
|
1554 xudong 1.6
|
1555 xudong 1.1 int nCutSegs = ARRLENGTH(CutSegs);
1556 for (int iSeg = 0; iSeg < nCutSegs; iSeg++) {
1557 DRMS_Segment_t *outSeg = drms_segment_lookupnum(sharpRec, iSeg);
1558 DRMS_Array_t *outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status);
1559 set_statistics(outSeg, outArray, 1);
1560 drms_free_array(outArray);
1561 }
|
1562 xudong 1.6
|
1563 xudong 1.1 return 0;
1564
1565 }
1566
1567
1568 /*
1569 * Get cutout and write segment
1570 * Change DISAMB_AZI to apply disambiguation to azimuth
1571 *
1572 */
1573
1574 int writeCutout(DRMS_Record_t *outRec, DRMS_Record_t *inRec, DRMS_Record_t *harpRec, char *SegName)
1575 {
1576
1577 int status = 0;
1578
1579 DRMS_Segment_t *inSeg = NULL, *outSeg = NULL;
1580 DRMS_Array_t *cutoutArray = NULL;
1581 // DRMS_Type_t arrayType;
1582
1583 int ll[2], ur[2], nx, ny, nxny; // lower-left and upper right coords
1584 xudong 1.1
1585 /* Info */
1586
1587 inSeg = drms_segment_lookup(inRec, SegName);
1588 if (!inSeg) return 1;
1589
1590 nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status);
1591 ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status);
1592 nxny = nx * ny;
1593 ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1;
1594 ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1;
1595 ur[0] = ll[0] + nx - 1; if (status) return 1;
1596 ur[1] = ll[1] + ny - 1; if (status) return 1;
1597
1598 if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) { // for bitmaps, infomaps, etc.
1599 cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status);
1600 if (status) return 1;
1601 } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) { // for full disk ones
1602 cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status);
1603 if (status) return 1;
1604 } else {
1605 xudong 1.1 return 1;
1606 }
|
1607 xudong 1.6
|
1608 xudong 1.1 /* Adding disambiguation resolution to cutout azimuth? */
|
1609 xudong 1.6
|
1610 xudong 1.1 #if DISAMB_AZI
1611 if (!strcmp(SegName, "azimuth")) {
|
1612 xudong 1.6 DRMS_Segment_t *disambSeg = NULL;
1613 disambSeg = drms_segment_lookup(inRec, "disambig");
|
1614 xudong 1.1 if (!disambSeg) {drms_free_array(cutoutArray); return 1;}
1615 DRMS_Array_t *disambArray;
1616 if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) {
1617 disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status);
1618 if (status) {drms_free_array(cutoutArray); return 1;}
1619 } else {
|
1620 xudong 1.6 drms_free_array(cutoutArray);
|
1621 xudong 1.1 return 1;
1622 }
1623 double *azimuth = (double *) cutoutArray->data;
1624 char *disamb = (char *) disambArray->data;
1625 for (int n = 0; n < nxny; n++) {
1626 if (disamb[n]) azimuth[n] += 180.;
1627 }
1628 drms_free_array(disambArray);
1629 }
1630 #endif
|
1631 xudong 1.6
|
1632 xudong 1.1 /* Write out */
1633
1634 outSeg = drms_segment_lookup(outRec, SegName);
1635 if (!outSeg) return 1;
|
1636 xudong 1.7 // drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray); // Jan 02 2013
|
1637 xudong 1.1 outSeg->axis[0] = cutoutArray->axis[0];
1638 outSeg->axis[1] = cutoutArray->axis[1];
|
1639 xudong 1.8 // cutoutArray->parent_segment = outSeg;
1640 cutoutArray->israw = 0; // always compressed
|
1641 xudong 1.1 cutoutArray->bzero = outSeg->bzero;
1642 cutoutArray->bscale = outSeg->bscale; // Same as inArray's
1643 status = drms_segment_write(outSeg, cutoutArray, 0);
1644 drms_free_array(cutoutArray);
1645 if (status) return 1;
1646
1647 return 0;
1648
1649 }
1650
1651
1652 /*
1653 * Compute space weather indices, no error checking for now
1654 * Based on M. Bobra's swharp_vectorB.c
1655 * No error checking for now
1656 *
1657 */
1658
1659 void computeSWIndex(struct swIndex *swKeys_ptr, DRMS_Record_t *inRec, struct mapInfo *mInfo)
1660 {
1661
1662 xudong 1.1 int status = 0;
1663 int nx = mInfo->ncol, ny = mInfo->nrow;
1664 int nxny = nx * ny;
1665 int dims[2] = {nx, ny};
1666 // Get bx, by, bz, mask
1667
|
1668 xudong 1.6 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
|
1669 mbobra 1.5 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(inRec, "bitmap");
1670 DRMS_Array_t *bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status);
1671 int *bitmask = (int *) bitmaskArray->data; // get the previously made mask array
|
1672 xudong 1.6
1673 //Use conf_disambig map as a threshold on spaceweather quantities
|
1674 xudong 1.2 DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
|
1675 xudong 1.1 DRMS_Array_t *maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
1676 int *mask = (int *) maskArray->data; // get the previously made mask array
|
1677 xudong 1.6
|
1678 xudong 1.1 DRMS_Segment_t *bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA);
1679 DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
1680 float *bx = (float *) bxArray->data; // bx
1681
1682 DRMS_Segment_t *bySeg = drms_segment_lookup(inRec, BT_SEG_CEA);
1683 DRMS_Array_t *byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status);
1684 float *by = (float *) byArray->data; // by
1685 for (int i = 0; i < nxny; i++) by[i] *= -1;
1686
1687 DRMS_Segment_t *bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA);
1688 DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
1689 float *bz = (float *) bzArray->data; // bz
1690
1691 // Get emphemeris
1692
1693 //float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);
|
1694 xudong 1.9 float cdelt1 = drms_getkey_float(inRec, "CDELT1", &status);
1695 float dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);
|
1696 xudong 1.1 double rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status);
1697 double rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status);
1698 float imcrpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status);
1699 float imcrpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status);
1700 float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status);
1701 float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status);
1702
1703 //float cdelt1=( (rsun_ref*cdelt1_orig*PI/180.) / (dsun_obs) )*(180./PI)*(3600.); //convert cdelt1 from degrees to arcsec (approximately)
|
1704 xudong 1.6
1705 printf("cdelt1=%f\n",cdelt1);
1706 printf("rsun_ref=%f\n",rsun_ref);
1707 printf("rsun_obs=%f\n",rsun_obs);
1708 printf("dsun_obs=%f\n",dsun_obs);
1709
|
1710 xudong 1.2 // Temp arrays
|
1711 xudong 1.1
1712 float *bh = (float *) (malloc(nxny * sizeof(float)));
1713 float *bt = (float *) (malloc(nxny * sizeof(float)));
1714 float *jz = (float *) (malloc(nxny * sizeof(float)));
|
1715 xudong 1.9 float *jz_smooth = (float *) (malloc(nxny * sizeof(float)));
|
1716 xudong 1.1 float *bpx = (float *) (malloc(nxny * sizeof(float)));
1717 float *bpy = (float *) (malloc(nxny * sizeof(float)));
1718 float *bpz = (float *) (malloc(nxny * sizeof(float)));
1719 float *derx = (float *) (malloc(nxny * sizeof(float)));
1720 float *dery = (float *) (malloc(nxny * sizeof(float)));
1721 float *derx_bt = (float *) (malloc(nxny * sizeof(float)));
1722 float *dery_bt = (float *) (malloc(nxny * sizeof(float)));
1723 float *derx_bh = (float *) (malloc(nxny * sizeof(float)));
1724 float *dery_bh = (float *) (malloc(nxny * sizeof(float)));
1725 float *derx_bz = (float *) (malloc(nxny * sizeof(float)));
1726 float *dery_bz = (float *) (malloc(nxny * sizeof(float)));
|
1727 xudong 1.9
1728 // The Computations
|
1729 xudong 1.1
1730 if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
|
1731 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)){
|
1732 xudong 1.1 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1733 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT;
1734 }
1735
1736 for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
|
1737 xudong 1.6 greenpot(bpx, bpy, bpz, nx, ny);
|
1738 xudong 1.1
|
1739 mbobra 1.5 computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
|
1740 xudong 1.1
|
1741 mbobra 1.5 if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))
|
1742 xudong 1.1 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT;
1743
|
1744 mbobra 1.5 computeB_total(bx, by, bz, bt, dims, mask, bitmask);
|
1745 xudong 1.1
|
1746 mbobra 1.5 if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))
|
1747 xudong 1.1 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT;
1748
|
1749 mbobra 1.5 if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))
|
1750 xudong 1.1 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT;
1751
|
1752 mbobra 1.5 if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))
|
1753 xudong 1.1 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1754
|
1755 xudong 1.9 computeJz(bx, by, dims, jz, mask, bitmask, cdelt1, rsun_ref, rsun_obs, derx, dery);
1756
1757 struct fresize_struct convolution_array;
1758 init_fresize_gaussian(&convolution_array, 4.0f, 12, 1);
1759 fresize(&convolution_array, jz, jz_smooth, nx, ny, nx, nx, ny, nx, 0, 0, 0.0f);
1760
1761 if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask,
1762 cdelt1, rsun_ref, rsun_obs, derx, dery))
1763 {
1764 swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
|
1765 xudong 1.1 swKeys_ptr->us_i = DRMS_MISSING_FLOAT;
1766 }
|
1767 xudong 1.9 printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz);
1768
|
1769 xudong 1.1
|
1770 xudong 1.9 if (computeAlpha(bz, dims, jz_smooth, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
|
1771 xudong 1.1 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT;
1772
|
1773 xudong 1.9 if (computeHelicity(bz, dims, jz_smooth, &(swKeys_ptr->mean_ih),
|
1774 xudong 1.1 &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
|
1775 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
|
1776 xudong 1.1 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT;
1777 swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT;
1778 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT;
1779 }
1780
|
1781 xudong 1.9 if (computeSumAbsPerPolarity(bz, jz_smooth, dims, &(swKeys_ptr->totaljz),
|
1782 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs))
|
1783 xudong 1.1 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
|
1784 xudong 1.6
|
1785 xudong 1.1
1786 if (computeFreeEnergy(bx, by, bpx, bpy, dims,
1787 &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
|
1788 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
|
1789 xudong 1.1 swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1790 swKeys_ptr->totpot = DRMS_MISSING_FLOAT;
1791 }
|
1792 xudong 1.9
|
1793 xudong 1.1 if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims,
1794 &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45),
1795 &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
|
1796 mbobra 1.5 mask, bitmask)) {
|
1797 xudong 1.1 swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1798 swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT;
1799 swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN
1800 swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT;
1801 }
|
1802 xudong 1.6
|
1803 xudong 1.1 // Clean up
1804
|
1805 xudong 1.6 drms_free_array(bitmaskArray); // Dec 18 2012 Xudong
|
1806 xudong 1.1 drms_free_array(maskArray);
|
1807 xudong 1.9 drms_free_array(bxArray);
|
1808 xudong 1.1 drms_free_array(byArray);
1809 drms_free_array(bzArray);
1810
|
1811 xudong 1.9 free(bh); free(bt); free(jz); free(jz_smooth);
|
1812 xudong 1.1 free(bpx); free(bpy); free(bpz);
1813 free(derx); free(dery);
1814 free(derx_bt); free(dery_bt);
1815 free(derx_bz); free(dery_bz);
1816 free(derx_bh); free(dery_bh);
1817
|
1818 xudong 1.9 free_fresize(&convolution_array);
|
1819 xudong 1.1 }
1820
1821
1822 /*
1823 * Set space weather indices, no error checking for now
1824 *
1825 */
1826
1827 void setSWIndex(DRMS_Record_t *outRec, struct swIndex *swKeys_ptr)
1828 {
1829 drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf);
1830 drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma);
1831 drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal);
1832 drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh);
1833 drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz);
1834 drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz);
1835 drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i);
1836 drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha);
1837 drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih);
1838 drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih);
1839 drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih);
1840 xudong 1.1 drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz);
1841 drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot);
1842 drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot);
1843 drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle);
1844 drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45);
1845 };
1846
1847
1848 /*
1849 * Set all keywords, no error checking for now
1850 *
1851 */
1852
|
1853 xudong 1.10 void setKeys(DRMS_Record_t *outRec, DRMS_Record_t *inRec, struct mapInfo *mInfo)
|
1854 xudong 1.1 {
|
1855 xudong 1.6 copy_me_keys(inRec, outRec);
1856 copy_patch_keys(inRec, outRec);
1857 copy_geo_keys(inRec, outRec);
1858 copy_ambig_keys(inRec, outRec);
1859
|
1860 xudong 1.10 int status = 0;
1861
1862 // Change a few geometry keywords for CEA records
1863 if (mInfo != NULL) {
1864
1865 drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5);
1866 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5);
1867
1868 drms_setkey_float(outRec, "CRVAL1", mInfo->xc);
1869 drms_setkey_float(outRec, "CRVAL2", mInfo->yc);
1870 drms_setkey_float(outRec, "CDELT1", mInfo->xscale);
1871 drms_setkey_float(outRec, "CDELT2", mInfo->yscale);
1872 drms_setkey_string(outRec, "CUNIT1", "degree");
1873 drms_setkey_string(outRec, "CUNIT2", "degree");
1874
1875 char key[64];
1876 snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]);
1877 drms_setkey_string(outRec, "CTYPE1", key);
1878 snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]);
1879 drms_setkey_string(outRec, "CTYPE2", key);
1880 drms_setkey_float(outRec, "CROTA2", 0.0);
1881 xudong 1.10
1882 } else {
1883
1884 float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status);
1885 float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status);
1886 float x_ll = drms_getkey_float(inRec, "CRPIX1", &status);
1887 float y_ll = drms_getkey_float(inRec, "CRPIX2", &status);
1888 // Defined as disk center's pixel address wrt lower-left of cutout
1889 drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.);
1890 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.);
1891 // Always 0.
1892 drms_setkey_float(outRec, "CRVAL1", 0);
1893 drms_setkey_float(outRec, "CRVAL2", 0);
1894
1895 }
1896
|
1897 xudong 1.6 char timebuf[1024];
1898 float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
1899 double val;
|
1900 xudong 1.10 status = 0;
1901
|
1902 xudong 1.6 val = drms_getkey_double(inRec, "DATE",&status);
1903 drms_setkey_double(outRec, "DATE_B", val);
1904 sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0);
1905 drms_setkey_string(outRec, "DATE", timebuf);
1906
1907 // set cvs commit version into keyword HEADER
|
1908 xudong 1.13 char *cvsinfo = strdup("$Id$");
|
1909 xudong 1.6 // status = drms_setkey_string(outRec, "HEADER", cvsinfo);
|
1910 xudong 1.4 status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);
|
1911 xudong 1.6
|
1912 xudong 1.1 };
1913
1914 //
1915 //
1916
1917 /* ############# Nearest neighbour interpolation ############### */
1918
1919 float nnb (float *f, int nx, int ny, double x, double y)
1920 {
1921
1922 if (x <= -0.5 || y <= -0.5 || x > nx - 0.5 || y > ny - 0.5)
1923 return DRMS_MISSING_FLOAT;
1924 int ilow = floor (x);
1925 int jlow = floor (y);
1926 int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow;
1927 int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow;
1928 return f[j * nx + i];
1929
1930 }
1931
1932 /* ################## Wrapper for Jesper's rebin code ################## */
1933 xudong 1.1
1934 void frebin (float *image_in, float *image_out, int nx, int ny, int nbin, int gauss)
1935 {
1936
1937 struct fresize_struct fresizes;
1938 int nxout, nyout, xoff, yoff;
1939 int nlead = nx;
1940
1941 nxout = nx / nbin; nyout = ny / nbin;
1942 if (gauss && nbin != 1)
1943 init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 * 2), nbin); // for nbin=3, sigma=1, half truncate width=2
1944 else
1945 init_fresize_bin(&fresizes, nbin);
1946 xoff = nbin / 2 + nbin / 2;
1947 yoff = nbin / 2 + nbin / 2;
1948 fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT);
1949
1950 }
|