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