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