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

Karen Tian