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

Karen Tian