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