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