JSOC/proj/sharp/apps/sharp.c - annotate

Return to sharp.c CVS log

Up to [Development] / JSOC / proj / sharp / apps

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_refcdelt1_origRADSINDEG)/(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

Karen Tian