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