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