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