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