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