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 xudong 1.13 // float psize = drms_getkey_float(inRec, "SIZE", &status); 932 // if (psize != psize) { 933 934 if (minlon != minlon \|\| maxlon != maxlon) { // check lons instead of SIZE
935 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
936 xudong 1.1 double omega = drms_getkey_double(inRec, "OMEGA_DT", &status); if (status) return 1; 937 char firstRecQuery[100], t0_str[100]; 938 sprint_time(t0_str, t0, "TAI", 0); 939 snprintf(firstRecQuery, 100, "%s[%d][%s]", inRec->seriesinfo->seriesname, harpnum, t0_str); 940 DRMS_RecordSet_t tmpRS = drms_open_records(drms_env, firstRecQuery, &status); 941 if (status \|\| tmpRS->n != 1) return 1; 942 DRMS_Record_t tmpRec = tmpRS->records[0]; 943 double minlon0 = drms_getkey_double(tmpRec, "LONDTMIN", &status); if (status) return 1; 944 double maxlon0 = drms_getkey_double(tmpRec, "LONDTMAX", &status); if (status) return 1; 945 minlon = minlon0 + (trec - t0) * omega / SECINDAY; 946 maxlon = maxlon0 + (trec - t0) * omega / SECINDAY; 947 printf("%s, %f, %f\n", firstRecQuery, minlon, maxlon); 948 } 949 950 mInfo->xc = (maxlon + minlon) / 2. + disk_lonc; 951 mInfo->yc = (maxlat + minlat) / 2.; 952 953 /* Size / 954 955 mInfo->ncol = round((maxlon - minlon) / mInfo->xscale); 956 mInfo->nrow = round((maxlat - minlat) / mInfo->yscale); 957 xudong 1.1 958 return 0; 959 960 } 961 962 963 / 964 * Fetch ephemeris info from a DRMS record 965 * No error checking for now 966 * 967 / 968 969 int getEphemeris(DRMS_Record_t inRec, struct ephemeris ephem) 970 { 971 972 int status = 0; 973 974 float crota2 = drms_getkey_float(inRec, "CROTA2", &status); // rotation 975 double sina = sin(crota2 RADSINDEG); 976 double cosa = cos(crota2 * RADSINDEG); 977 978 xudong 1.1 ephem->pa = - crota2 * RADSINDEG; 979 ephem->disk_latc = drms_getkey_float(inRec, "CRLT_OBS", &status) * RADSINDEG; 980 ephem->disk_lonc = drms_getkey_float(inRec, "CRLN_OBS", &status) * RADSINDEG; 981 982 float crvalx = 0.0; 983 float crvaly = 0.0; 984 float crpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status); 985 float crpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status); 986 float cdelt = drms_getkey_float(inRec, "CDELT1", &status); // in arcsec, assumimg dx=dy
987 xudong 1.6 printf("cdelt=%f\n",cdelt);
988 xudong 1.1 ephem->disk_xc = PIX_X(0.0,0.0) - 1.0; // Center of disk in pixel, starting at 0 989 ephem->disk_yc = PIX_Y(0.0,0.0) - 1.0; 990 991 float dSun = drms_getkey_float(inRec, "DSUN_OBS", &status); 992 float rSun_ref = drms_getkey_float(inRec, "RSUN_REF", &status); 993 if (status) rSun_ref = 6.96e8; 994 995 ephem->asd = asin(rSun_ref/dSun); 996 ephem->rSun = asin(rSun_ref / dSun) * RAD2ARCSEC / cdelt; 997 998 return 0; 999 1000 } 1001 1002 1003 /* 1004 * Compute the coordinates to be sampled on full disk image 1005 * mInfo->xi_out & mInfo->zeta_out 1006 * This is oversampled, its size is ncol0 & nrow0 as shown below 1007 * 1008 * 1009 xudong 1.1 / 1010 1011 void findCoord(struct mapInfo mInfo) 1012 { 1013 1014 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN 1015 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2; 1016 1017 float xscale0 = mInfo->xscale / mInfo->nbin * RADSINDEG; // oversampling resolution 1018 float yscale0 = mInfo->yscale / mInfo->nbin * RADSINDEG; // in rad 1019 1020 double lonc = mInfo->xc * RADSINDEG; // in rad 1021 double latc = mInfo->yc * RADSINDEG; 1022 1023 double disk_lonc = (mInfo->ephem).disk_lonc; 1024 double disk_latc = (mInfo->ephem).disk_latc; 1025 1026 double rSun = (mInfo->ephem).rSun; 1027 double disk_xc = (mInfo->ephem).disk_xc / rSun; 1028 double disk_yc = (mInfo->ephem).disk_yc / rSun; 1029 double pa = (mInfo->ephem).pa; 1030 xudong 1.1 1031 // Temp pointers 1032 1033 float xi_out = mInfo->xi_out; 1034 float zeta_out = mInfo->zeta_out; 1035 1036 // start 1037 1038 double x, y; // map coord 1039 double lat, lon; // helio coord 1040 double xi, zeta; // image coord (for one point) 1041 1042 int ind_map; 1043 1044 for (int row0 = 0; row0 < nrow0; row0++) { 1045 for (int col0 = 0; col0 < ncol0; col0++) { 1046 1047 ind_map = row0 * ncol0 + col0; 1048 1049 x = (col0 + 0.5 - ncol0/2.) * xscale0; // in rad 1050 y = (row0 + 0.5 - nrow0/2.) * yscale0; 1051 xudong 1.1 1052 /* map grid [x, y] corresponds to the point [lon, lat] in the heliographic coordinates. 1053 * the [x, y] are in radians with respect of the center of the map [xcMap, ycMap]. 1054 * projection methods could be Mercator, Lambert, and many others. [maplonc, mapLatc] 1055 * is the heliographic longitude and latitude of the map center. Both are in degree. 1056 / 1057 1058 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) { 1059 xi_out[ind_map] = -1; 1060 zeta_out[ind_map] = -1; 1061 continue; 1062 } 1063 1064 / map the grid [lon, lat] in the heliographic coordinates to [xi, zeta], a point in the 1065 * image coordinates. The image properties, xCenter, yCenter, rSun, pa, ecc and chi are given. 1066 / 1067 1068 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta, 1069 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) { 1070 xi_out[ind_map] = -1; 1071 zeta_out[ind_map] = -1; 1072 xudong 1.1 continue; 1073 } 1074 1075 xi_out[ind_map] = xi rSun; 1076 zeta_out[ind_map] = zeta * rSun; 1077 1078 } 1079 } 1080 1081 } 1082 1083 1084 /* 1085 * Sampling function 1086 * oversampling by nbin, then binning using a Gaussian 1087 * save results in outData, always of float type 1088 * 1089 / 1090 1091 int performSampling(float outData, float inData, struct mapInfo mInfo) 1092 { 1093 xudong 1.1 1094 int status = 0; 1095 1096 int ncol0 = mInfo->ncol * mInfo->nbin + (mInfo->nbin / 2) * 2; // pad with nbin/2 on edge to avoid NAN 1097 int nrow0 = mInfo->nrow * mInfo->nbin + (mInfo->nbin / 2) * 2; 1098 1099 float outData0 = (float ) (malloc(ncol0 * nrow0 * sizeof(float))); 1100 1101 float xi_out = mInfo->xi_out; 1102 float zeta_out = mInfo->zeta_out;
1103 xudong 1.6
1104 xudong 1.1 // Interpolation 1105 1106 struct fint_struct pars; 1107 int interpOpt = INTERP; // Use Wiener by default, 6 order, 1 constraint 1108 1109 switch (interpOpt) { 1110 case 0: // Wiener, 6 order, 1 constraint 1111 init_finterpolate_wiener(&pars, 6, 1, 6, 2, 1, 1, NULL, dpath); 1112 break; 1113 case 1: // Cubic convolution 1114 init_finterpolate_cubic_conv(&pars, 1., 3.); 1115 break; 1116 case 2: // Bilinear 1117 init_finterpolate_linear(&pars, 1.); 1118 break; 1119 default: 1120 return 1; 1121 } 1122 1123 finterpolate(&pars, inData, xi_out, zeta_out, outData0, 1124 FOURK, FOURK, FOURK, ncol0, nrow0, ncol0, DRMS_MISSING_FLOAT); 1125 xudong 1.1 1126 // Rebinning, smoothing 1127 1128 frebin(outData0, outData, ncol0, nrow0, mInfo->nbin, 1); // Gaussian
1129 xudong 1.6 free(outData0); // Dec 18 2012
1130 xudong 1.1 1131 // 1132 1133 return 0; 1134 1135 } 1136 1137 1138 /* 1139 * Performing local vector transformation 1140 * xyz: z refers to vertical (radial) component, x EW (phi), y NS 1141 * 1142 / 1143 1144 void vectorTransform(float bx_map, float by_map, float bz_map, struct mapInfo mInfo) 1145 { 1146 1147 int ncol = mInfo->ncol; 1148 int nrow = mInfo->nrow; 1149 1150 float xscale = mInfo->xscale RADSINDEG; // in rad 1151 xudong 1.1 float yscale = mInfo->yscale * RADSINDEG; 1152 1153 double lonc = mInfo->xc * RADSINDEG; // in rad 1154 double latc = mInfo->yc * RADSINDEG; 1155 1156 double disk_lonc = (mInfo->ephem).disk_lonc; 1157 double disk_latc = (mInfo->ephem).disk_latc; 1158 1159 double rSun = (mInfo->ephem).rSun; 1160 double disk_xc = (mInfo->ephem).disk_xc / rSun; 1161 double disk_yc = (mInfo->ephem).disk_yc / rSun; 1162 double pa = (mInfo->ephem).pa; 1163 1164 int ind_map; 1165 double x, y; 1166 double lat, lon; // lat / lon for current point 1167 1168 double bx_tmp, by_tmp, bz_tmp; 1169 1170 // 1171 1172 xudong 1.1 for (int row = 0; row < mInfo->nrow; row++) { 1173 for (int col = 0; col < mInfo->ncol; col++) { 1174 1175 ind_map = row * mInfo->ncol + col; 1176 1177 x = (col + 0.5 - ncol / 2.) * xscale; 1178 y = (row + 0.5 - nrow / 2.) * yscale; 1179 1180 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) { 1181 bx_map[ind_map] = DRMS_MISSING_FLOAT; 1182 by_map[ind_map] = DRMS_MISSING_FLOAT; 1183 bz_map[ind_map] = DRMS_MISSING_FLOAT; 1184 continue; 1185 } 1186 1187 bx_tmp = by_tmp = bz_tmp = 0; 1188 1189 img2helioVector (bx_map[ind_map], by_map[ind_map], bz_map[ind_map], 1190 &bx_tmp, &by_tmp, &bz_tmp, 1191 lon, lat, disk_lonc, disk_latc, pa); 1192 1193 xudong 1.1 bx_map[ind_map] = bx_tmp; 1194 by_map[ind_map] = by_tmp; 1195 bz_map[ind_map] = bz_tmp; 1196 1197 } 1198 }
1199 xudong 1.6
1200 xudong 1.1 } 1201 1202 1203 1204 /* 1205 * Map and propogate vector field errors 1206 * 1207 / 1208 1209 int getBErr(float bx_err, float by_err, float bz_err,
1210 xudong 1.6 DRMS_Record_t inRec, struct mapInfo mInfo)
1211 xudong 1.1 { 1212 1213 int status = 0; 1214 1215 // Get variances and covariances, filling factor assume to be 1 1216 1217 float bT = (float ) (malloc(FOURK2 * sizeof(float))); // field 1218 float bI = (float ) (malloc(FOURK2 * sizeof(float))); // inclination 1219 float bA = (float ) (malloc(FOURK2 * sizeof(float))); // azimuth 1220 1221 float errbT = (float ) (malloc(FOURK2 * sizeof(float))); 1222 float errbI = (float ) (malloc(FOURK2 * sizeof(float))); 1223 float errbA = (float ) (malloc(FOURK2 * sizeof(float))); 1224 1225 float errbTbI = (float ) (malloc(FOURK2 * sizeof(float))); 1226 float errbTbA = (float ) (malloc(FOURK2 * sizeof(float))); 1227 float errbIbA = (float ) (malloc(FOURK2 * sizeof(float))); 1228 1229 if (readVectorBErr(inRec, 1230 bT, bI, bA, 1231 errbT, errbI, errbA, 1232 xudong 1.1 errbTbI, errbTbA, errbIbA)) { 1233 printf("Read full disk variances & covariances error\n"); 1234 free(bT); free(bI); free(bA); 1235 free(errbT); free(errbI); free(errbA); 1236 free(errbTbI); free(errbTbA); free(errbIbA); 1237 return 1; 1238 } 1239 1240 // Size 1241 1242 int ncol = mInfo->ncol; 1243 int nrow = mInfo->nrow; 1244 1245 float xscale = mInfo->xscale * RADSINDEG; // in rad 1246 float yscale = mInfo->yscale * RADSINDEG; 1247 1248 double lonc = mInfo->xc * RADSINDEG; // in rad 1249 double latc = mInfo->yc * RADSINDEG; 1250 1251 double disk_lonc = (mInfo->ephem).disk_lonc; 1252 double disk_latc = (mInfo->ephem).disk_latc; 1253 xudong 1.1 1254 double rSun = (mInfo->ephem).rSun; 1255 double disk_xc = (mInfo->ephem).disk_xc / rSun; 1256 double disk_yc = (mInfo->ephem).disk_yc / rSun; 1257 double pa = (mInfo->ephem).pa; 1258 1259 // Start 1260 1261 double x, y; // map coord 1262 double lat, lon; // spherical coord 1263 double xi, zeta; // image coord, round to full pixel 1264 1265 int ind_map, ind_img; 1266 1267 double bpSigma2, btSigma2, brSigma2; // variances after prop
1268 xudong 1.6
1269 xudong 1.1 for (int row = 0; row < mInfo->nrow; row++) { 1270 for (int col = 0; col < mInfo->ncol; col++) { 1271 1272 ind_map = row * mInfo->ncol + col; 1273 1274 x = (col + 0.5 - ncol / 2.) * xscale; 1275 y = (row + 0.5 - nrow / 2.) * yscale; 1276 1277 if (plane2sphere (x, y, latc, lonc, &lat, &lon, (int) mInfo->proj)) { 1278 bx_err[ind_map] = DRMS_MISSING_FLOAT; 1279 by_err[ind_map] = DRMS_MISSING_FLOAT; 1280 bz_err[ind_map] = DRMS_MISSING_FLOAT; 1281 continue; 1282 } 1283 1284 if (sphere2img (lat, lon, disk_latc, disk_lonc, &xi, &zeta, 1285 disk_xc, disk_yc, 1.0, pa, 0., 0., 0., 0.)) { 1286 bx_err[ind_map] = DRMS_MISSING_FLOAT; 1287 bx_err[ind_map] = DRMS_MISSING_FLOAT; 1288 bx_err[ind_map] = DRMS_MISSING_FLOAT; 1289 continue; 1290 xudong 1.1 } 1291 1292 xi = rSun; xi = round(xi); 1293 zeta = rSun; zeta = round(zeta); // nearest neighbor 1294 1295 ind_img = round(zeta * FOURK + xi); 1296 1297 if (errorprop(bT, bA, bI, 1298 errbT, errbA, errbI, errbTbA, errbTbI, errbIbA, 1299 lon, lat, disk_lonc, disk_latc, pa, FOURK, FOURK, xi, zeta, 1300 &btSigma2, &bpSigma2, &brSigma2)) { 1301 bx_err[ind_map] = DRMS_MISSING_FLOAT; 1302 by_err[ind_map] = DRMS_MISSING_FLOAT; 1303 bz_err[ind_map] = DRMS_MISSING_FLOAT; 1304 continue; 1305 } 1306 1307 bx_err[ind_map] = sqrt(bpSigma2); 1308 by_err[ind_map] = sqrt(btSigma2); 1309 bz_err[ind_map] = sqrt(brSigma2); 1310 1311 xudong 1.1 } 1312 } 1313 1314 // 1315 1316 free(bT); free(bI); free(bA); 1317 free(errbT); free(errbI); free(errbA); 1318 free(errbTbI); free(errbTbA); free(errbIbA); 1319 return 0; 1320 1321 } 1322 1323 1324 1325 /* 1326 * Read full disk vector magnetograms 1327 * Fill factor is 1, use default disambiguity resolution 1328 * 1329 / 1330 1331 int readVectorB(DRMS_Record_t inRec, float bx_img, float by_img, float bz_img) 1332 xudong 1.1 { 1333 1334 int status = 0; 1335 1336 DRMS_Segment_t inSeg; 1337 DRMS_Array_t *inArray_ambig;
1338 xudong 1.6 DRMS_Array_t inArray_bTotal, inArray_bAzim, *inArray_bIncl;
1339 xudong 1.1 1340 char ambig; 1341 float bTotal, bAzim, bIncl; 1342 1343 inSeg = drms_segment_lookup(inRec, "disambig"); 1344 inArray_ambig = drms_segment_read(inSeg, DRMS_TYPE_CHAR, &status); 1345 if (status) return 1; 1346 ambig = (char )inArray_ambig->data; 1347 1348 inSeg = drms_segment_lookup(inRec, "field"); 1349 inArray_bTotal = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status); 1350 if (status) return 1; 1351 bTotal = (float )inArray_bTotal->data; 1352 1353 inSeg = drms_segment_lookup(inRec, "azimuth"); 1354 inArray_bAzim = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status); 1355 if (status) return 1; 1356 bAzim = (float )inArray_bAzim->data; 1357 1358 inSeg = drms_segment_lookup(inRec, "inclination"); 1359 inArray_bIncl = drms_segment_read(inSeg, DRMS_TYPE_FLOAT, &status); 1360 xudong 1.1 if (status) return 1; 1361 bIncl = (float )inArray_bIncl->data; 1362 1363 // Convert CCD xyz 1364 1365 int llx, lly; // lower-left corner 1366 int bmx, bmy; // bitmap size 1367 1368 llx = (int)(drms_getkey_float(inRec, "CRPIX1", &status)) - 1; 1369 lly = (int)(drms_getkey_float(inRec, "CRPIX2", &status)) - 1; 1370 1371 bmx = inArray_ambig->axis[0]; 1372 bmy = inArray_ambig->axis[1]; 1373 1374 int kx, ky, kOff; 1375 int ix = 0, jy = 0, yOff = 0, iData = 0; 1376 int xDim = FOURK, yDim = FOURK; 1377 1378 for (jy = 0; jy < yDim; jy++) 1379 { 1380 ix = 0; 1381 xudong 1.1 yOff = jy * xDim; 1382 ky = jy - lly; 1383 for (ix = 0; ix < xDim; ix++) 1384 { 1385 iData = yOff + ix; 1386 kx = ix - llx; 1387 1388 // zero azi pointing up, zero incl pointing out from sun 1389 bx_img[iData] = - bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * sin(bAzim[iData] * RADSINDEG); 1390 by_img[iData] = bTotal[iData] * sin(bIncl[iData] * RADSINDEG) * cos(bAzim[iData] * RADSINDEG); 1391 bz_img[iData] = bTotal[iData] * cos(bIncl[iData] * RADSINDEG); 1392 1393 // Disambiguation 1394 1395 if (kx < 0 \|\| kx >= bmx \|\| ky < 0 \|\| ky >= bmy) { 1396 continue; 1397 } else { 1398 kOff = ky * bmx + kx; 1399 if (ambig[kOff] % 2) { // 180 1400 bx_img[iData] = -1.; by_img[iData] = -1.; 1401 } 1402 xudong 1.1 } 1403 } 1404 } 1405 1406 // Clean up 1407 1408 drms_free_array(inArray_ambig); 1409 drms_free_array(inArray_bTotal); 1410 drms_free_array(inArray_bAzim); 1411 drms_free_array(inArray_bIncl); 1412 1413 return 0; 1414 1415 } 1416 1417 1418 /* 1419 * Read variances and covariances of vector magnetograms 1420 * 1421 / 1422 1423 xudong 1.1 int readVectorBErr(DRMS_Record_t inRec, 1424 float bT, float bI, float bA, 1425 float errbT, float errbI, float errbA, 1426 float errbTbI, float errbTbA, float errbIbA) 1427 { 1428 1429 int status = 0; 1430 1431 float data_ptr[9]; 1432 char *segName[9] = {"field", "inclination", "azimuth",
1433 xudong 1.6 "field_err", "inclination_err", "azimuth_err", 1434 "field_inclination_err", "field_az_err", "inclin_azimuth_err"};
1435 xudong 1.1 DRMS_Segment_t inSegs[9]; 1436 DRMS_Array_t inArrays[9]; 1437 1438 // Read full disk images 1439 1440 for (int iSeg = 0; iSeg < 9; iSeg++) { 1441 1442 inSegs[iSeg] = drms_segment_lookup(inRec, segName[iSeg]); 1443 inArrays[iSeg] = drms_segment_read(inSegs[iSeg], DRMS_TYPE_FLOAT, &status); 1444 data_ptr[iSeg] = (float ) inArrays[iSeg]->data; 1445 1446 } 1447 1448 float bT0 = data_ptr[0], bI0 = data_ptr[1], bA0 = data_ptr[2]; 1449 float errbT0 = data_ptr[3], errbI0 = data_ptr[4], errbA0 = data_ptr[5]; 1450 float errbTbI0 = data_ptr[6], errbTbA0 = data_ptr[7], errbIbA0 = data_ptr[8]; 1451 1452 // Convert errors to variances, correlation coefficients to covariances 1453 1454 for (int i = 0; i < FOURK2; i++) { 1455 1456 xudong 1.1 if (fabs(errbI0[i]) > 180.) errbI0[i] = 180.; 1457 if (fabs(errbA0[i]) > 180.) errbA0[i] = 180.; 1458 1459 bT[i] = bT0[i]; 1460 bI[i] = bI0[i]; 1461 bA[i] = bA0[i]; 1462 1463 errbT[i] = errbT0[i] * errbT0[i]; 1464 errbI[i] = errbI0[i] * errbI0[i] * RADSINDEG * RADSINDEG; 1465 errbA[i] = errbA0[i] * errbA0[i] * RADSINDEG * RADSINDEG; 1466 1467 errbTbI[i] = errbTbI0[i] * errbT0[i] * errbI0[i] * RADSINDEG; 1468 errbTbA[i] = errbTbA0[i] * errbT0[i] * errbA0[i] * RADSINDEG; 1469 errbIbA[i] = errbIbA0[i] * errbI0[i] * errbA0[i] * RADSINDEG * RADSINDEG;
1470 xudong 1.6
1471 xudong 1.1 } 1472 1473 // 1474 1475 for (int iSeg = 0; iSeg < 9; iSeg++) drms_free_array(inArrays[iSeg]);
1476 xudong 1.6
1477 xudong 1.1 return 0; 1478 1479 } 1480 1481 1482 /* 1483 * Create Cutout record: top level subroutine 1484 * Do the loops on segments and set the keywords here 1485 * Work is done in writeCutout routine below 1486 * 1487 / 1488 1489 int createCutRecord(DRMS_Record_t mharpRec, DRMS_Record_t bharpRec, 1490 DRMS_Record_t dopRec, DRMS_Record_t contRec, 1491 DRMS_Record_t sharpRec, struct swIndex *swKeys_ptr) 1492 { 1493 1494 int status = 0; 1495 1496 int iHarpSeg; 1497 int nMharpSegs = ARRLENGTH(MharpSegs), nBharpSegs = ARRLENGTH(BharpSegs); 1498 xudong 1.1 1499 // Cutout Mharp 1500 1501 for (iHarpSeg = 0; iHarpSeg < nMharpSegs; iHarpSeg++) { 1502 if (writeCutout(sharpRec, mharpRec, mharpRec, MharpSegs[iHarpSeg])) { 1503 printf("Mharp cutout fails for %s\n", MharpSegs[iHarpSeg]); 1504 break; 1505 } 1506 }
1507 xudong 1.3 if (iHarpSeg != nMharpSegs) { 1508 SHOW("Cutout: segment number unmatch\n"); 1509 return 1; // if failed 1510 }
1511 xudong 1.1 printf("Magnetogram cutout done.\n"); 1512 1513 // Cutout Doppler 1514 1515 if (writeCutout(sharpRec, dopRec, mharpRec, "Dopplergram")) { 1516 printf("Doppler cutout failed\n"); 1517 return 1; 1518 } 1519 printf("Dopplergram cutout done.\n"); 1520 1521 // Cutout Continuum 1522 1523 if (writeCutout(sharpRec, contRec, mharpRec, "continuum")) { 1524 printf("Continuum cutout failed\n"); 1525 return 1; 1526 } 1527 printf("Intensitygram cutout done.\n"); 1528 1529 // Coutout Bharp 1530 1531 for (iHarpSeg = 0; iHarpSeg < nBharpSegs; iHarpSeg++) { 1532 xudong 1.1 if (writeCutout(sharpRec, bharpRec, mharpRec, BharpSegs[iHarpSeg])) { 1533 printf("Bharp cutout fails for %s\n", BharpSegs[iHarpSeg]); 1534 break; 1535 } 1536 } 1537 if (iHarpSeg != nBharpSegs) return 1; // if failed 1538 printf("Vector B cutout done.\n"); 1539 1540 // Keywords & Links 1541 1542 drms_copykey(sharpRec, mharpRec, "T_REC"); 1543 drms_copykey(sharpRec, mharpRec, "HARPNUM"); 1544 1545 DRMS_Link_t mHarpLink = hcon_lookup_lower(&sharpRec->links, "MHARP"); 1546 if (mHarpLink) drms_link_set("MHARP", sharpRec, mharpRec); 1547 DRMS_Link_t bHarpLink = hcon_lookup_lower(&sharpRec->links, "BHARP"); 1548 if (bHarpLink) drms_link_set("BHARP", sharpRec, bharpRec); 1549 1550 setSWIndex(sharpRec, swKeys_ptr); // Set space weather indices
1551 xudong 1.10 setKeys(sharpRec, bharpRec, NULL); // Set all other keywords, NULL specifies cutout
1552 xudong 1.6
1553 xudong 1.1 // Stats
1554 xudong 1.6
1555 xudong 1.1 int nCutSegs = ARRLENGTH(CutSegs); 1556 for (int iSeg = 0; iSeg < nCutSegs; iSeg++) { 1557 DRMS_Segment_t outSeg = drms_segment_lookupnum(sharpRec, iSeg); 1558 DRMS_Array_t outArray = drms_segment_read(outSeg, DRMS_TYPE_FLOAT, &status); 1559 set_statistics(outSeg, outArray, 1); 1560 drms_free_array(outArray); 1561 }
1562 xudong 1.6
1563 xudong 1.1 return 0; 1564 1565 } 1566 1567 1568 /* 1569 * Get cutout and write segment 1570 * Change DISAMB_AZI to apply disambiguation to azimuth 1571 * 1572 / 1573 1574 int writeCutout(DRMS_Record_t outRec, DRMS_Record_t inRec, DRMS_Record_t harpRec, char SegName) 1575 { 1576 1577 int status = 0; 1578 1579 DRMS_Segment_t inSeg = NULL, outSeg = NULL; 1580 DRMS_Array_t cutoutArray = NULL; 1581 // DRMS_Type_t arrayType; 1582 1583 int ll[2], ur[2], nx, ny, nxny; // lower-left and upper right coords 1584 xudong 1.1 1585 /* Info / 1586 1587 inSeg = drms_segment_lookup(inRec, SegName); 1588 if (!inSeg) return 1; 1589 1590 nx = (int) drms_getkey_float(harpRec, "CRSIZE1", &status); 1591 ny = (int) drms_getkey_float(harpRec, "CRSIZE2", &status); 1592 nxny = nx ny; 1593 ll[0] = (int) drms_getkey_float(harpRec, "CRPIX1", &status) - 1; if (status) return 1; 1594 ll[1] = (int) drms_getkey_float(harpRec, "CRPIX2", &status) - 1; if (status) return 1; 1595 ur[0] = ll[0] + nx - 1; if (status) return 1; 1596 ur[1] = ll[1] + ny - 1; if (status) return 1; 1597 1598 if (inSeg->axis[0] == nx && inSeg->axis[1] == ny) { // for bitmaps, infomaps, etc. 1599 cutoutArray = drms_segment_read(inSeg, DRMS_TYPE_DOUBLE, &status); 1600 if (status) return 1; 1601 } else if (inSeg->axis[0] == FOURK && inSeg->axis[1] == FOURK) { // for full disk ones 1602 cutoutArray = drms_segment_readslice(inSeg, DRMS_TYPE_DOUBLE, ll, ur, &status); 1603 if (status) return 1; 1604 } else { 1605 xudong 1.1 return 1; 1606 }
1607 xudong 1.6
1608 xudong 1.1 /* Adding disambiguation resolution to cutout azimuth? */
1609 xudong 1.6
1610 xudong 1.1 #if DISAMB_AZI 1611 if (!strcmp(SegName, "azimuth")) {
1612 xudong 1.6 DRMS_Segment_t *disambSeg = NULL; 1613 disambSeg = drms_segment_lookup(inRec, "disambig");
1614 xudong 1.1 if (!disambSeg) {drms_free_array(cutoutArray); return 1;} 1615 DRMS_Array_t *disambArray; 1616 if (disambSeg->axis[0] == nx && disambSeg->axis[1] == ny) { 1617 disambArray = drms_segment_read(disambSeg, DRMS_TYPE_CHAR, &status); 1618 if (status) {drms_free_array(cutoutArray); return 1;} 1619 } else {
1620 xudong 1.6 drms_free_array(cutoutArray);
1621 xudong 1.1 return 1; 1622 } 1623 double azimuth = (double ) cutoutArray->data; 1624 char disamb = (char ) disambArray->data; 1625 for (int n = 0; n < nxny; n++) { 1626 if (disamb[n]) azimuth[n] += 180.; 1627 } 1628 drms_free_array(disambArray); 1629 } 1630 #endif
1631 xudong 1.6
1632 xudong 1.1 /* Write out */ 1633 1634 outSeg = drms_segment_lookup(outRec, SegName); 1635 if (!outSeg) return 1;
1636 xudong 1.7 // drms_array_convert_inplace(outSeg->info->type, 0, 1, cutoutArray); // Jan 02 2013
1637 xudong 1.1 outSeg->axis[0] = cutoutArray->axis[0]; 1638 outSeg->axis[1] = cutoutArray->axis[1];
1639 xudong 1.8 // cutoutArray->parent_segment = outSeg; 1640 cutoutArray->israw = 0; // always compressed
1641 xudong 1.1 cutoutArray->bzero = outSeg->bzero; 1642 cutoutArray->bscale = outSeg->bscale; // Same as inArray's 1643 status = drms_segment_write(outSeg, cutoutArray, 0); 1644 drms_free_array(cutoutArray); 1645 if (status) return 1; 1646 1647 return 0; 1648 1649 } 1650 1651 1652 /* 1653 * Compute space weather indices, no error checking for now 1654 * Based on M. Bobra's swharp_vectorB.c 1655 * No error checking for now 1656 * 1657 / 1658 1659 void computeSWIndex(struct swIndex swKeys_ptr, DRMS_Record_t inRec, struct mapInfo mInfo) 1660 { 1661 1662 xudong 1.1 int status = 0; 1663 int nx = mInfo->ncol, ny = mInfo->nrow; 1664 int nxny = nx * ny; 1665 int dims[2] = {nx, ny}; 1666 // Get bx, by, bz, mask 1667
1668 xudong 1.6 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
1669 mbobra 1.5 DRMS_Segment_t bitmaskSeg = drms_segment_lookup(inRec, "bitmap"); 1670 DRMS_Array_t bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status); 1671 int bitmask = (int ) bitmaskArray->data; // get the previously made mask array
1672 xudong 1.6 1673 //Use conf_disambig map as a threshold on spaceweather quantities
1674 xudong 1.2 DRMS_Segment_t *maskSeg = drms_segment_lookup(inRec, "conf_disambig");
1675 xudong 1.1 DRMS_Array_t maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); 1676 int mask = (int *) maskArray->data; // get the previously made mask array
1677 xudong 1.6
1678 xudong 1.1 DRMS_Segment_t bxSeg = drms_segment_lookup(inRec, BP_SEG_CEA); 1679 DRMS_Array_t bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status); 1680 float bx = (float ) bxArray->data; // bx 1681 1682 DRMS_Segment_t bySeg = drms_segment_lookup(inRec, BT_SEG_CEA); 1683 DRMS_Array_t byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status); 1684 float by = (float ) byArray->data; // by 1685 for (int i = 0; i < nxny; i++) by[i] = -1; 1686 1687 DRMS_Segment_t bzSeg = drms_segment_lookup(inRec, BR_SEG_CEA); 1688 DRMS_Array_t bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status); 1689 float bz = (float *) bzArray->data; // bz 1690 1691 // Get emphemeris 1692 1693 //float cdelt1_orig = drms_getkey_float(inRec, "CDELT1", &status);
1694 xudong 1.9 float cdelt1 = drms_getkey_float(inRec, "CDELT1", &status); 1695 float dsun_obs = drms_getkey_float(inRec, "DSUN_OBS", &status);
1696 xudong 1.1 double rsun_ref = drms_getkey_double(inRec, "RSUN_REF", &status); 1697 double rsun_obs = drms_getkey_double(inRec, "RSUN_OBS", &status); 1698 float imcrpix1 = drms_getkey_float(inRec, "IMCRPIX1", &status); 1699 float imcrpix2 = drms_getkey_float(inRec, "IMCRPIX2", &status); 1700 float crpix1 = drms_getkey_float(inRec, "CRPIX1", &status); 1701 float crpix2 = drms_getkey_float(inRec, "CRPIX2", &status); 1702 1703 //float cdelt1=( (rsun_refcdelt1_origPI/180.) / (dsun_obs) )(180./PI)(3600.); //convert cdelt1 from degrees to arcsec (approximately)
1704 xudong 1.6 1705 printf("cdelt1=%f\n",cdelt1); 1706 printf("rsun_ref=%f\n",rsun_ref); 1707 printf("rsun_obs=%f\n",rsun_obs); 1708 printf("dsun_obs=%f\n",dsun_obs); 1709
1710 xudong 1.2 // Temp arrays
1711 xudong 1.1 1712 float bh = (float ) (malloc(nxny * sizeof(float))); 1713 float bt = (float ) (malloc(nxny * sizeof(float))); 1714 float jz = (float ) (malloc(nxny * sizeof(float)));
1715 xudong 1.9 float jz_smooth = (float ) (malloc(nxny * sizeof(float)));
1716 xudong 1.1 float bpx = (float ) (malloc(nxny * sizeof(float))); 1717 float bpy = (float ) (malloc(nxny * sizeof(float))); 1718 float bpz = (float ) (malloc(nxny * sizeof(float))); 1719 float derx = (float ) (malloc(nxny * sizeof(float))); 1720 float dery = (float ) (malloc(nxny * sizeof(float))); 1721 float derx_bt = (float ) (malloc(nxny * sizeof(float))); 1722 float dery_bt = (float ) (malloc(nxny * sizeof(float))); 1723 float derx_bh = (float ) (malloc(nxny * sizeof(float))); 1724 float dery_bh = (float ) (malloc(nxny * sizeof(float))); 1725 float derx_bz = (float ) (malloc(nxny * sizeof(float))); 1726 float dery_bz = (float ) (malloc(nxny * sizeof(float)));
1727 xudong 1.9 1728 // The Computations
1729 xudong 1.1 1730 if (computeAbsFlux(bz, dims, &(swKeys_ptr->absFlux), &(swKeys_ptr->mean_vf),
1731 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)){
1732 xudong 1.1 swKeys_ptr->absFlux = DRMS_MISSING_FLOAT; // If fail, fill in NaN 1733 swKeys_ptr->mean_vf = DRMS_MISSING_FLOAT; 1734 } 1735 1736 for (int i = 0; i < nxny; i++) bpz[i] = bz[i];
1737 xudong 1.6 greenpot(bpx, bpy, bpz, nx, ny);
1738 xudong 1.1
1739 mbobra 1.5 computeBh(bx, by, bz, bh, dims, &(swKeys_ptr->mean_hf), mask, bitmask);
1740 xudong 1.1
1741 mbobra 1.5 if (computeGamma(bx, by, bz, bh, dims, &(swKeys_ptr->mean_gamma), mask, bitmask))
1742 xudong 1.1 swKeys_ptr->mean_gamma = DRMS_MISSING_FLOAT; 1743
1744 mbobra 1.5 computeB_total(bx, by, bz, bt, dims, mask, bitmask);
1745 xudong 1.1
1746 mbobra 1.5 if (computeBtotalderivative(bt, dims, &(swKeys_ptr->mean_derivative_btotal), mask, bitmask, derx_bt, dery_bt))
1747 xudong 1.1 swKeys_ptr->mean_derivative_btotal = DRMS_MISSING_FLOAT; 1748
1749 mbobra 1.5 if (computeBhderivative(bh, dims, &(swKeys_ptr->mean_derivative_bh), mask, bitmask, derx_bh, dery_bh))
1750 xudong 1.1 swKeys_ptr->mean_derivative_bh = DRMS_MISSING_FLOAT; 1751
1752 mbobra 1.5 if (computeBzderivative(bz, dims, &(swKeys_ptr->mean_derivative_bz), mask, bitmask, derx_bz, dery_bz))
1753 xudong 1.1 swKeys_ptr->mean_derivative_bz = DRMS_MISSING_FLOAT; // If fail, fill in NaN 1754
1755 xudong 1.9 computeJz(bx, by, dims, jz, mask, bitmask, cdelt1, rsun_ref, rsun_obs, derx, dery); 1756 1757 struct fresize_struct convolution_array; 1758 init_fresize_gaussian(&convolution_array, 4.0f, 12, 1); 1759 fresize(&convolution_array, jz, jz_smooth, nx, ny, nx, nx, ny, nx, 0, 0, 0.0f); 1760 1761 if(computeJzsmooth(bx, by, dims, jz_smooth, &(swKeys_ptr->mean_jz), &(swKeys_ptr->us_i), mask, bitmask, 1762 cdelt1, rsun_ref, rsun_obs, derx, dery)) 1763 { 1764 swKeys_ptr->mean_jz = DRMS_MISSING_FLOAT;
1765 xudong 1.1 swKeys_ptr->us_i = DRMS_MISSING_FLOAT; 1766 }
1767 xudong 1.9 printf("swKeys_ptr->mean_jz=%f\n",swKeys_ptr->mean_jz); 1768
1769 xudong 1.1
1770 xudong 1.9 if (computeAlpha(bz, dims, jz_smooth, &(swKeys_ptr->mean_alpha), mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1771 xudong 1.1 swKeys_ptr->mean_alpha = DRMS_MISSING_FLOAT; 1772
1773 xudong 1.9 if (computeHelicity(bz, dims, jz_smooth, &(swKeys_ptr->mean_ih),
1774 xudong 1.1 &(swKeys_ptr->total_us_ih), &(swKeys_ptr->total_abs_ih),
1775 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
1776 xudong 1.1 swKeys_ptr->mean_ih = DRMS_MISSING_FLOAT; 1777 swKeys_ptr->total_us_ih = DRMS_MISSING_FLOAT; 1778 swKeys_ptr->total_abs_ih = DRMS_MISSING_FLOAT; 1779 } 1780
1781 xudong 1.9 if (computeSumAbsPerPolarity(bz, jz_smooth, dims, &(swKeys_ptr->totaljz),
1782 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs))
1783 xudong 1.1 swKeys_ptr->totaljz = DRMS_MISSING_FLOAT;
1784 xudong 1.6
1785 xudong 1.1 1786 if (computeFreeEnergy(bx, by, bpx, bpy, dims, 1787 &(swKeys_ptr->meanpot), &(swKeys_ptr->totpot),
1788 mbobra 1.5 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) {
1789 xudong 1.1 swKeys_ptr->meanpot = DRMS_MISSING_FLOAT; // If fail, fill in NaN 1790 swKeys_ptr->totpot = DRMS_MISSING_FLOAT; 1791 }
1792 xudong 1.9
1793 xudong 1.1 if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, 1794 &(swKeys_ptr->meanshear_angle), &(swKeys_ptr->area_w_shear_gt_45), 1795 &(swKeys_ptr->meanshear_angleh), &(swKeys_ptr->area_w_shear_gt_45h),
1796 mbobra 1.5 mask, bitmask)) {
1797 xudong 1.1 swKeys_ptr->meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN 1798 swKeys_ptr->area_w_shear_gt_45 = DRMS_MISSING_FLOAT; 1799 swKeys_ptr->meanshear_angleh = DRMS_MISSING_FLOAT; // If fail, fill in NaN 1800 swKeys_ptr->area_w_shear_gt_45h = DRMS_MISSING_FLOAT; 1801 }
1802 xudong 1.6
1803 xudong 1.1 // Clean up 1804
1805 xudong 1.6 drms_free_array(bitmaskArray); // Dec 18 2012 Xudong
1806 xudong 1.1 drms_free_array(maskArray);
1807 xudong 1.9 drms_free_array(bxArray);
1808 xudong 1.1 drms_free_array(byArray); 1809 drms_free_array(bzArray); 1810
1811 xudong 1.9 free(bh); free(bt); free(jz); free(jz_smooth);
1812 xudong 1.1 free(bpx); free(bpy); free(bpz); 1813 free(derx); free(dery); 1814 free(derx_bt); free(dery_bt); 1815 free(derx_bz); free(dery_bz); 1816 free(derx_bh); free(dery_bh); 1817
1818 xudong 1.9 free_fresize(&convolution_array);
1819 xudong 1.1 } 1820 1821 1822 /* 1823 * Set space weather indices, no error checking for now 1824 * 1825 / 1826 1827 void setSWIndex(DRMS_Record_t outRec, struct swIndex swKeys_ptr) 1828 { 1829 drms_setkey_float(outRec, "USFLUX", swKeys_ptr->mean_vf); 1830 drms_setkey_float(outRec, "MEANGAM", swKeys_ptr->mean_gamma); 1831 drms_setkey_float(outRec, "MEANGBT", swKeys_ptr->mean_derivative_btotal); 1832 drms_setkey_float(outRec, "MEANGBH", swKeys_ptr->mean_derivative_bh); 1833 drms_setkey_float(outRec, "MEANGBZ", swKeys_ptr->mean_derivative_bz); 1834 drms_setkey_float(outRec, "MEANJZD", swKeys_ptr->mean_jz); 1835 drms_setkey_float(outRec, "TOTUSJZ", swKeys_ptr->us_i); 1836 drms_setkey_float(outRec, "MEANALP", swKeys_ptr->mean_alpha); 1837 drms_setkey_float(outRec, "MEANJZH", swKeys_ptr->mean_ih); 1838 drms_setkey_float(outRec, "TOTUSJH", swKeys_ptr->total_us_ih); 1839 drms_setkey_float(outRec, "ABSNJZH", swKeys_ptr->total_abs_ih); 1840 xudong 1.1 drms_setkey_float(outRec, "SAVNCPP", swKeys_ptr->totaljz); 1841 drms_setkey_float(outRec, "MEANPOT", swKeys_ptr->meanpot); 1842 drms_setkey_float(outRec, "TOTPOT", swKeys_ptr->totpot); 1843 drms_setkey_float(outRec, "MEANSHR", swKeys_ptr->meanshear_angle); 1844 drms_setkey_float(outRec, "SHRGT45", swKeys_ptr->area_w_shear_gt_45); 1845 }; 1846 1847 1848 / 1849 * Set all keywords, no error checking for now 1850 * 1851 */ 1852
1853 xudong 1.10 void setKeys(DRMS_Record_t outRec, DRMS_Record_t inRec, struct mapInfo *mInfo)
1854 xudong 1.1 {
1855 xudong 1.6 copy_me_keys(inRec, outRec); 1856 copy_patch_keys(inRec, outRec); 1857 copy_geo_keys(inRec, outRec); 1858 copy_ambig_keys(inRec, outRec); 1859
1860 xudong 1.10 int status = 0; 1861 1862 // Change a few geometry keywords for CEA records 1863 if (mInfo != NULL) { 1864 1865 drms_setkey_float(outRec, "CRPIX1", mInfo->ncol/2. + 0.5); 1866 drms_setkey_float(outRec, "CRPIX2", mInfo->nrow/2. + 0.5); 1867 1868 drms_setkey_float(outRec, "CRVAL1", mInfo->xc); 1869 drms_setkey_float(outRec, "CRVAL2", mInfo->yc); 1870 drms_setkey_float(outRec, "CDELT1", mInfo->xscale); 1871 drms_setkey_float(outRec, "CDELT2", mInfo->yscale); 1872 drms_setkey_string(outRec, "CUNIT1", "degree"); 1873 drms_setkey_string(outRec, "CUNIT2", "degree"); 1874 1875 char key[64]; 1876 snprintf (key, 64, "CRLN-%s", wcsCode[(int) mInfo->proj]); 1877 drms_setkey_string(outRec, "CTYPE1", key); 1878 snprintf (key, 64, "CRLT-%s", wcsCode[(int) mInfo->proj]); 1879 drms_setkey_string(outRec, "CTYPE2", key); 1880 drms_setkey_float(outRec, "CROTA2", 0.0); 1881 xudong 1.10 1882 } else { 1883 1884 float disk_xc = drms_getkey_float(inRec, "IMCRPIX1", &status); 1885 float disk_yc = drms_getkey_float(inRec, "IMCRPIX2", &status); 1886 float x_ll = drms_getkey_float(inRec, "CRPIX1", &status); 1887 float y_ll = drms_getkey_float(inRec, "CRPIX2", &status); 1888 // Defined as disk center's pixel address wrt lower-left of cutout 1889 drms_setkey_float(outRec, "CRPIX1", disk_xc - x_ll + 1.); 1890 drms_setkey_float(outRec, "CRPIX2", disk_yc - y_ll + 1.); 1891 // Always 0. 1892 drms_setkey_float(outRec, "CRVAL1", 0); 1893 drms_setkey_float(outRec, "CRVAL2", 0); 1894 1895 } 1896
1897 xudong 1.6 char timebuf[1024]; 1898 float UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */ 1899 double val;
1900 xudong 1.10 status = 0; 1901
1902 xudong 1.6 val = drms_getkey_double(inRec, "DATE",&status); 1903 drms_setkey_double(outRec, "DATE_B", val); 1904 sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0); 1905 drms_setkey_string(outRec, "DATE", timebuf); 1906 1907 // set cvs commit version into keyword HEADER
1908 xudong 1.13 char *cvsinfo = strdup("$Id$");
1909 xudong 1.6 // status = drms_setkey_string(outRec, "HEADER", cvsinfo);
1910 xudong 1.4 status = drms_setkey_string(outRec, "CODEVER7", cvsinfo);
1911 xudong 1.6
1912 xudong 1.1 }; 1913 1914 // 1915 // 1916 1917 /* ############# Nearest neighbour interpolation ############### / 1918 1919 float nnb (float f, int nx, int ny, double x, double y) 1920 { 1921 1922 if (x <= -0.5 \|\| y <= -0.5 \|\| x > nx - 0.5 \|\| y > ny - 0.5) 1923 return DRMS_MISSING_FLOAT; 1924 int ilow = floor (x); 1925 int jlow = floor (y); 1926 int i = ((x - ilow) > 0.5) ? ilow + 1 : ilow; 1927 int j = ((y - jlow) > 0.5) ? jlow + 1 : jlow; 1928 return f[j * nx + i]; 1929 1930 } 1931 1932 /* ################## Wrapper for Jesper's rebin code ################## / 1933 xudong 1.1 1934 void frebin (float image_in, float image_out, int nx, int ny, int nbin, int gauss) 1935 { 1936 1937 struct fresize_struct fresizes; 1938 int nxout, nyout, xoff, yoff; 1939 int nlead = nx; 1940 1941 nxout = nx / nbin; nyout = ny / nbin; 1942 if (gauss && nbin != 1) 1943 init_fresize_gaussian(&fresizes, (nbin / 2), (nbin / 2 2), nbin); // for nbin=3, sigma=1, half truncate width=2 1944 else 1945 init_fresize_bin(&fresizes, nbin); 1946 xoff = nbin / 2 + nbin / 2; 1947 yoff = nbin / 2 + nbin / 2; 1948 fresize(&fresizes, image_in, image_out, nx, ny, nlead, nxout, nyout, nxout, xoff, yoff, DRMS_MISSING_FLOAT); 1949 1950 }

Karen Tian