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

Return to update_sharp_keys.c CVS log

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

1 mbobra 1.1 /* 2 * MODULE NAME: update_sharp_keys.c 3 * 4 * DESCRIPTION: This module recalculates SHARP keywords. This is accomplished by 5 * cloning a record, recalculating keywords of choice (i.e., user input), 6 * and pointing to the same segments as the old record. Associated error keys 7 * are computed by default. 8 * 9 * This module accounts for versioning by appending to the keyword CODEVER7 and HISTORY: 10 * CODEVER7 will contain multiple lines: the production build of sharp.c, the 11 * production build of include file sw_functions.c, and the production build 12 * of update_sharp_keys.c. HISTORY will include a human-readable sentence about which 13 * keywords were updated. 14 * 15 * This module does not produce segments. 16 * 17 * INPUTS : -- DRMS SHARP series 18 * -- DRMS SHARP CEA series 19 * -- HARPNUM 20 * -- comma separated list of keywords to recalculate 21 * 22 mbobra 1.1 * AUTHOR : Monica Bobra 23 * 24 * Version : v0.0 Jun 14 2013 25 * 26 * EXAMPLE :
27 arta 1.5 * update_sharp_keys sharpseries=hmi.sharp_720s sharpceaseries=hmi.sharp_cea_720s // 28 * HARPNUM=1 keylist=USFLUX,TOTPOT
29 mbobra 1.1 * 30 / 31 32 / Include files / 33 #include <stdio.h> 34 #include <stdlib.h> 35 #include <time.h> 36 #include <sys/time.h> 37 #include <math.h> 38 #include <string.h> 39 #include "jsoc_main.h" 40 #include "astro.h" 41 #include "fstats.h" 42 #include "cartography.c" 43 #include "fresize.h" 44 #include "finterpolate.h" 45 #include "img2helioVector.c" 46 #include "copy_me_keys.c" 47 #include "errorprop.c" 48 #include "sw_functions.c" 49 #include <mkl_blas.h> 50 mbobra 1.1 #include <mkl_service.h> 51 #include <mkl_lapack.h> 52 #include <mkl_vml_functions.h> 53 #include <omp.h> 54 55 / Define values / 56 #define PI (M_PI) 57 #define RADSINDEG (PI/180.) 58 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0])) 59 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);} 60 #define SHOW(msg) {printf("%s", msg); fflush(stdout);} 61 62 char module_name = "update_sharp_keys"; /* Module name */
63 arta 1.5 char version_id = "2013 Jun 14"; / Version number */
64 mbobra 1.1 65 ModuleArgs_t module_args[] = 66 {
67 arta 1.5 {ARG_STRING, "sharpseries", NULL, "Input/output Sharp dataseries"}, 68 {ARG_STRING, "sharpceaseries", NULL, "Input/output Sharp CEA dataseries"},
69 mbobra 1.1 {ARG_INT, "HARPNUM", NULL, "HARP number"}, 70 {ARG_STRING, "keylist", NULL, "comma separated list of keywords to update"}, 71 {ARG_END} 72 }; 73 74 75 int DoIt(void) 76 { 77 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); 78 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); 79 80 int status = DRMS_SUCCESS; 81 int nrecs, irec; 82 83 /* Keywords */ 84 float mean_vf; 85 float absFlux; 86 float count_mask; 87 float mean_hf; 88 float mean_gamma; 89 float mean_derivative_btotal; 90 mbobra 1.1 float mean_derivative_bh; 91 float mean_derivative_bz; 92 float mean_jz; 93 float us_i; 94 float mean_alpha; 95 float mean_ih; 96 float total_us_ih; 97 float total_abs_ih; 98 float totaljz; 99 float totpot; 100 float meanpot; 101 float area_w_shear_gt_45; 102 float meanshear_angle; 103 float area_w_shear_gt_45h; 104 float meanshear_angleh; 105 float mean_derivative_btotal_err; 106 float mean_vf_err; 107 float mean_gamma_err; 108 float mean_derivative_bh_err; 109 float mean_derivative_bz_err; 110 float mean_jz_err; 111 mbobra 1.1 float us_i_err; 112 float mean_alpha_err; 113 float mean_ih_err; 114 float total_us_ih_err; 115 float total_abs_ih_err; 116 float totaljz_err; 117 float meanpot_err; 118 float totpot_err; 119 float meanshear_angle_err;
120 arta 1.5 char sharpseries = (char ) params_get_str(&cmdparams, "sharpseries"); 121 char sharpceaseries = (char ) params_get_str(&cmdparams, "sharpceaseries");
122 mbobra 1.1 123 int harpnum = params_get_int(&cmdparams, "HARPNUM"); 124 char sharpquery[256]; 125 char sharpceaquery[256]; 126
127 arta 1.5 sprintf(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum); 128 printf(sharpceaquery, "%s[%d]\n", sharpceaseries,harpnum); 129 130 sprintf(sharpquery, "%s[%d]\n", sharpseries,harpnum); 131 printf(sharpquery, "%s[%d]\n", sharpseries,harpnum);
132 mbobra 1.1
133 arta 1.5 DRMS_RecordSet_t *sharpinrecset = drms_open_records(drms_env, sharpquery, &status);
134 mbobra 1.1 if (status != DRMS_SUCCESS) 135 DIE("Problem opening sharp recordset.");
136 arta 1.5 if (sharpinrecset->n == 0)
137 mbobra 1.1 DIE("Sharp recordset contains no records.");
138 arta 1.5 DRMS_RecordSet_t *sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status); 139 if (status != DRMS_SUCCESS) 140 DIE("Problem cloning sharp records.");
141 mbobra 1.1 142 DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status); 143 if (status != DRMS_SUCCESS) 144 DIE("Problem opening sharp cea recordset.");
145 arta 1.5 if (sharpceainrecset->n == 0) 146 DIE("Sharp cea recordset contains no records."); 147 DRMS_RecordSet_t *sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
148 mbobra 1.1 if (status != DRMS_SUCCESS) 149 DIE("Problem cloning sharp cea records."); 150 151 char keylist = (char ) params_get_str(&cmdparams, "keylist"); 152 153 // Flags to indicate which keyword will be recalculated 154 int meanvfflag = (strstr(keylist,"USFLUX") != NULL); // generalize so that lowercase is acceptable 155 int totpotflag = (strstr(keylist,"TOTPOT") != NULL); 156 int meangamflag = (strstr(keylist,"MEANGAM") != NULL); 157 int meangbtflag = (strstr(keylist,"MEANGBT") != NULL); 158 int meangbzflag = (strstr(keylist,"MEANGBZ") != NULL); 159 int meangbhflag = (strstr(keylist,"MEANGBH") != NULL); 160 int meanjzdflag = (strstr(keylist,"MEANJZD") != NULL); 161 int totusjzflag = (strstr(keylist,"TOTUSJZ") != NULL); 162 int meanalpflag = (strstr(keylist,"MEANALP") != NULL); 163 int meanjzhflag = (strstr(keylist,"MEANJZH") != NULL); 164 int totusjhflag = (strstr(keylist,"TOTUSJH") != NULL); 165 int absnjzhflag = (strstr(keylist,"ABSNJZH") != NULL); 166 int savncppflag = (strstr(keylist,"SAVNCPP") != NULL); 167 int meanpotflag = (strstr(keylist,"MEANPOT") != NULL); 168 int meanshrflag = (strstr(keylist,"MEANSHR") != NULL); 169 mbobra 1.1 int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL); 170 DRMS_Record_t sharpinrec = sharpinrecset->records[0]; 171 DRMS_Record_t sharpceainrec = sharpceainrecset->records[0]; 172 DRMS_Segment_t inseg = drms_segment_lookup(sharpceainrec, "Br"); 173 int nx = inseg->axis[0]; 174 int ny = inseg->axis[1]; 175 int nxny = nx ny; 176 int dims[2] = {nx, ny};
177 arta 1.5
178 mbobra 1.1 // Temp arrays 179 float bh = (float ) (malloc(nxny * sizeof(float))); 180 float bt = (float ) (malloc(nxny * sizeof(float))); 181 float jz = (float ) (malloc(nxny * sizeof(float))); 182 float bpx = (float ) (malloc(nxny * sizeof(float))); 183 float bpy = (float ) (malloc(nxny * sizeof(float))); 184 float bpz = (float ) (malloc(nxny * sizeof(float))); 185 float derx = (float ) (malloc(nxny * sizeof(float))); 186 float dery = (float ) (malloc(nxny * sizeof(float))); 187 float derx_bt = (float ) (malloc(nxny * sizeof(float))); 188 float dery_bt = (float ) (malloc(nxny * sizeof(float))); 189 float derx_bh = (float ) (malloc(nxny * sizeof(float))); 190 float dery_bh = (float ) (malloc(nxny * sizeof(float))); 191 float derx_bz = (float ) (malloc(nxny * sizeof(float))); 192 float dery_bz = (float ) (malloc(nxny * sizeof(float))); 193 float bt_err = (float ) (malloc(nxny * sizeof(float))); 194 float bh_err = (float ) (malloc(nxny * sizeof(float))); 195 float jz_err = (float ) (malloc(nxny * sizeof(float))); 196 float jz_err_squared = (float ) (malloc(nxny * sizeof(float))); 197 float jz_rms_err = (float ) (malloc(nxny * sizeof(float))); 198 float jz_err_squared_smooth = (float ) (malloc(nxny * sizeof(float))); 199 mbobra 1.1 float jz_smooth = (float ) (malloc(nxny * sizeof(float))); 200 201 // ephemeris variables 202 float cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2; 203 double rsun_ref, rsun_obs; 204 205 // outrecords
206 arta 1.5 DRMS_Record_t sharpoutrec, sharpceaoutrec; 207 nrecs=sharpoutrecset->n; 208 printf("nrecs=%d\n",nrecs);
209 mbobra 1.1 210 // prepare to set CODEVER7 (CVS Version of the SHARP module) 211 char cvsinfo0; 212 char history0;
213 mbobra 1.4 char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.3 2013/08/07 00:04:46 arta Exp $");
214 mbobra 1.1 char cvsinfo2 = sw_functions_version(); 215 char cvsinfoall = (char *)malloc(2048); 216 char historyofthemodule[2048]; 217 cvsinfo0 = drms_getkey_string(sharpinrec, "CODEVER7", &status); 218 strcpy(cvsinfoall,cvsinfo0); 219 strcat(cvsinfoall,"\n"); 220 strcat(cvsinfoall,cvsinfo2); 221 222 // prepare to set HISTORY (History of the data) 223 char timebuf[1024]; 224 float UNIX_epoch = -220924792.000; 225 sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0); 226 sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist); 227 printf("historyofthemodule=%s\n",historyofthemodule); 228 history0 = drms_getkey_string(sharpinrec, "HISTORY", &status); 229 strcat(historyofthemodule,"\n"); 230 strcat(historyofthemodule,history0); 231
232 arta 1.5 // no longer need inrecs 233 drms_close_records(sharpinrecset, DRMS_FREE_RECORD); 234 drms_close_records(sharpceainrecset, DRMS_FREE_RECORD); 235
236 mbobra 1.1 for (irec=0;irec<nrecs;irec++) 237 { 238 // Get emphemeris 239 sharpoutrec = sharpoutrecset->records[irec]; 240 sharpceaoutrec = sharpceaoutrecset->records[irec];
241 arta 1.5 cdelt1_orig = drms_getkey_float(sharpceaoutrec, "CDELT1", &status); 242 dsun_obs = drms_getkey_float(sharpoutrec, "DSUN_OBS", &status); 243 rsun_ref = drms_getkey_double(sharpoutrec, "RSUN_REF", &status); 244 rsun_obs = drms_getkey_double(sharpoutrec, "RSUN_OBS", &status); 245 imcrpix1 = drms_getkey_float(sharpoutrec, "IMCRPIX1", &status); 246 imcrpix2 = drms_getkey_float(sharpoutrec, "IMCRPIX2", &status); 247 crpix1 = drms_getkey_float(sharpoutrec, "CRPIX1", &status); 248 crpix2 = drms_getkey_float(sharpoutrec, "CRPIX2", &status); 249 cdelt1 = (atan((rsun_refcdelt1_origRADSINDEG)/(dsun_obs)))(1/RADSINDEG)(3600.);
250 mbobra 1.1
251 arta 1.5 // set CODEVER7 and HISTORY
252 mbobra 1.1 drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall); 253 drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule); 254 drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall); 255 drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule); 256 257 // Get bx, by, bz, mask 258 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
259 arta 1.5 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceaoutrec, "bitmap");
260 mbobra 1.1 DRMS_Array_t bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status); 261 int bitmask = (int *) bitmaskArray->data; // get the previously made mask array 262 263 //Use conf_disambig map as a threshold on spaceweather quantities
264 arta 1.5 DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceaoutrec, "conf_disambig");
265 mbobra 1.1 DRMS_Array_t maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); 266 if (status != DRMS_SUCCESS) 267 DIE("No CONF_DISAMBIG segment."); 268 int mask = (int *) maskArray->data; // get the previously made mask array 269
270 arta 1.5 DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceaoutrec, "Bp");
271 mbobra 1.1 DRMS_Array_t bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status); 272 float bx = (float *) bxArray->data; // bx 273
274 arta 1.5 DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceaoutrec, "Bt");
275 mbobra 1.1 DRMS_Array_t byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status); 276 float by = (float ) byArray->data; // by 277 for (int i = 0; i < nxny; i++) by[i] = -1; 278
279 arta 1.5 DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceaoutrec, "Br");
280 mbobra 1.1 DRMS_Array_t bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status); 281 float bz = (float *) bzArray->data; // bz 282
283 arta 1.5 DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceaoutrec, "Br_err");
284 mbobra 1.1 DRMS_Array_t bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status); 285 float bz_err = (float *) bz_errArray->data; // bz_err 286
287 arta 1.5 DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceaoutrec, "Bt_err");
288 mbobra 1.1 DRMS_Array_t by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status); 289 float by_err = (float *) by_errArray->data; // by_err 290
291 arta 1.5 DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceaoutrec, "Bp_err");
292 mbobra 1.1 DRMS_Array_t bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status); 293 float bx_err = (float *) bx_errArray->data; // bx_err 294
295 arta 1.5 /*** USFLUX, Example Function 1 *********************************sdfdsf/
296 mbobra 1.1 if (meanvfflag) 297 { 298 // Compute unsigned flux 299 if (computeAbsFlux(bz_err, bz , dims, &absFlux, &mean_vf, &mean_vf_err, 300 &count_mask, mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 301 { 302 mean_vf = DRMS_MISSING_FLOAT; 303 mean_vf_err = DRMS_MISSING_FLOAT; 304 count_mask = DRMS_MISSING_INT; 305 } 306 307 drms_setkey_float(sharpoutrec, "USFLUX", mean_vf); 308 drms_setkey_float(sharpoutrec, "CMASK", count_mask); 309 drms_setkey_float(sharpoutrec, "ERRVF", mean_vf_err); 310 311 drms_setkey_float(sharpceaoutrec, "USFLUX", mean_vf); 312 drms_setkey_float(sharpceaoutrec, "CMASK", count_mask); 313 drms_setkey_float(sharpceaoutrec, "ERRVF", mean_vf_err); 314 }
315 arta 1.5
316 mbobra 1.1 /*** MEANPOT and TOTPOT, Example Function 13 (Requires Keiji's code) / 317 318 if (totpotflag \|\| meanpotflag) 319 { 320 // First compute potential field 321 for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; 322 greenpot(bpx, bpy, bpz, nx, ny); 323 324 // Then compute energy 325 if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, 326 &meanpot, &meanpot_err, &totpot, &totpot_err, 327 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 328 { 329 meanpot = DRMS_MISSING_FLOAT; 330 totpot = DRMS_MISSING_FLOAT; 331 meanpot_err = DRMS_MISSING_FLOAT; 332 totpot_err = DRMS_MISSING_FLOAT; 333 } 334 335 // Set sharp keys 336 drms_setkey_float(sharpoutrec, "MEANPOT", meanpot); 337 mbobra 1.1 drms_setkey_float(sharpoutrec, "TOTPOT", totpot); 338 drms_setkey_float(sharpoutrec, "ERRMPOT", meanpot_err); 339 drms_setkey_float(sharpoutrec, "ERRTPOT", totpot_err); 340 341 // Set sharp cea keys 342 drms_setkey_float(sharpceaoutrec, "MEANPOT", meanpot); 343 drms_setkey_float(sharpceaoutrec, "TOTPOT", totpot); 344 drms_setkey_float(sharpceaoutrec, "ERRMPOT", meanpot_err); 345 drms_setkey_float(sharpceaoutrec, "ERRTPOT", totpot_err); 346 } 347 348 /*** MEANGAM, Example Function 3 ********************************/ 349 350 if (meangamflag) 351 { 352 // First compute horizontal field 353 computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask); 354 355 // Then compute inclination angle, gamma 356 if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &mean_gamma, &mean_gamma_err,mask, bitmask)) 357 { 358 mbobra 1.1 mean_gamma = DRMS_MISSING_FLOAT; 359 mean_gamma_err = DRMS_MISSING_FLOAT; 360 } 361 362 // Set sharp keys 363 drms_setkey_float(sharpoutrec, "MEANGAM", mean_gamma); 364 drms_setkey_float(sharpoutrec, "ERRGAM", mean_gamma_err); 365 366 // Set sharp cea keys 367 drms_setkey_float(sharpceaoutrec, "MEANGAM", mean_gamma); 368 drms_setkey_float(sharpceaoutrec, "ERRGAM", mean_gamma_err); 369 } 370 371 /* MEANGBT, Example Function 5 (Requires Function 4) *********/ 372 373 if (meangbtflag) 374 { 375 // First, compute Bt 376 computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask); 377 378 // Then take the derivative of Bt 379 mbobra 1.1 if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err, 380 &mean_derivative_btotal_err)) 381 { 382 mean_derivative_btotal = DRMS_MISSING_FLOAT; 383 mean_derivative_btotal_err = DRMS_MISSING_FLOAT; 384 } 385 386 // Set sharp keys 387 drms_setkey_float(sharpoutrec, "MEANGBT", mean_derivative_btotal); 388 drms_setkey_float(sharpoutrec, "ERRBT", mean_derivative_btotal_err); 389 390 // Set sharp cea keys 391 drms_setkey_float(sharpceaoutrec, "MEANGBT", mean_derivative_btotal); 392 drms_setkey_float(sharpceaoutrec, "ERRBT", mean_derivative_btotal_err); 393 } 394 395 /* MEANGBH, Example Function 6 (Requires Function 2) *********/ 396 397 if (meangbhflag) 398 { 399 // First, compute Bh 400 mbobra 1.1 computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask); 401 402 // Then take the derivative of Bh 403 if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh)) 404 { 405 mean_derivative_bh = DRMS_MISSING_FLOAT; 406 mean_derivative_bh_err = DRMS_MISSING_FLOAT; 407 } 408 409 // Set sharp keys 410 drms_setkey_float(sharpoutrec, "MEANGBH", mean_derivative_bh); 411 drms_setkey_float(sharpoutrec, "ERRBH", mean_derivative_bh_err); 412 413 // Set sharp cea keys 414 drms_setkey_float(sharpceaoutrec, "MEANGBH", mean_derivative_bh); 415 drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err); 416 } 417 418 /* MEANGBZ, Example Function 7 ********************************/ 419 420 if (meangbzflag) 421 mbobra 1.1 { 422 // Compute Bz derivative 423 if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz)) 424 { 425 mean_derivative_bz = DRMS_MISSING_FLOAT; 426 mean_derivative_bz_err = DRMS_MISSING_FLOAT; 427 } 428 429 // Set sharp keys 430 drms_setkey_float(sharpoutrec, "MEANGBZ", mean_derivative_bz); 431 drms_setkey_float(sharpoutrec, "ERRBZ", mean_derivative_bz_err); 432 433 // Set sharp cea keys 434 drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz); 435 drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err); 436 } 437 438 /* MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) / 439 440 if (meanjzdflag \|\| totusjzflag) 441 { 442 mbobra 1.1 // First, compute Jz on all pixels 443 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 444 derx, dery); 445 446 // Then take the sums of the appropriate pixels 447 if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz, 448 &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1, 449 rsun_ref, rsun_obs, derx, dery)) 450 451 452 { 453 mean_jz = DRMS_MISSING_FLOAT; 454 us_i = DRMS_MISSING_FLOAT; 455 mean_jz_err = DRMS_MISSING_FLOAT; 456 us_i_err = DRMS_MISSING_FLOAT; 457 } 458 459 // Set sharp keys 460 drms_setkey_float(sharpoutrec, "MEANJZD", mean_jz); 461 drms_setkey_float(sharpoutrec, "TOTUSJZ", us_i); 462 drms_setkey_float(sharpoutrec, "ERRJZ", mean_jz_err); 463 mbobra 1.1 drms_setkey_float(sharpoutrec, "ERRUSI", us_i_err); 464 465 // Set sharp cea keys 466 drms_setkey_float(sharpceaoutrec, "MEANJZD", mean_jz); 467 drms_setkey_float(sharpceaoutrec, "TOTUSJZ", us_i); 468 drms_setkey_float(sharpceaoutrec, "ERRJZ", mean_jz_err); 469 drms_setkey_float(sharpceaoutrec, "ERRUSI", us_i_err); 470 } 471 472 /** MEANALP, Example Function 10 (Requires Function 8)*****/ 473 474 if (meanalpflag) 475 { 476 // First, compute Jz on all pixels 477 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 478 derx, dery); 479 480 // Then compute alpha quantities on select pixels 481 if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err, 482 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 483 { 484 mbobra 1.1 mean_alpha = DRMS_MISSING_FLOAT; 485 mean_alpha_err = DRMS_MISSING_FLOAT; 486 } 487 488 // Set sharp keys 489 drms_setkey_float(sharpoutrec, "MEANALP", mean_alpha); 490 drms_setkey_float(sharpoutrec, "ERRALP", mean_alpha_err); 491 492 // Set sharp cea keys 493 drms_setkey_float(sharpceaoutrec, "MEANALP", mean_alpha); 494 drms_setkey_float(sharpceaoutrec, "ERRALP", mean_alpha_err); 495 496 } 497 498 /* MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) / 499 500 if (meanjzhflag \|\| totusjhflag \|\| absnjzhflag) 501 { 502 // First, compute Jz on all pixels 503 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 504 derx, dery); 505 mbobra 1.1 506 // Then compute helicity quantities on select pixels 507 if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih, 508 &total_abs_ih, &total_us_ih_err, &total_abs_ih_err, mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 509 { 510 mean_ih = DRMS_MISSING_FLOAT; 511 total_us_ih = DRMS_MISSING_FLOAT; 512 total_abs_ih = DRMS_MISSING_FLOAT; 513 mean_ih_err = DRMS_MISSING_FLOAT; 514 total_us_ih_err = DRMS_MISSING_FLOAT; 515 total_abs_ih_err = DRMS_MISSING_FLOAT; 516 } 517 518 // Set sharp keys 519 drms_setkey_float(sharpoutrec, "MEANJZH", mean_ih); 520 drms_setkey_float(sharpoutrec, "TOTUSJH", total_us_ih); 521 drms_setkey_float(sharpoutrec, "ABSNJZH", total_abs_ih); 522 drms_setkey_float(sharpoutrec, "ERRMIH", mean_ih_err); 523 drms_setkey_float(sharpoutrec, "ERRTUI", total_us_ih_err); 524 drms_setkey_float(sharpoutrec, "ERRTAI", total_abs_ih_err); 525 526 mbobra 1.1 // Set sharp cea keys 527 drms_setkey_float(sharpceaoutrec, "MEANJZH", mean_ih); 528 drms_setkey_float(sharpceaoutrec, "TOTUSJH", total_us_ih); 529 drms_setkey_float(sharpceaoutrec, "ABSNJZH", total_abs_ih); 530 drms_setkey_float(sharpceaoutrec, "ERRMIH", mean_ih_err); 531 drms_setkey_float(sharpceaoutrec, "ERRTUI", total_us_ih_err); 532 drms_setkey_float(sharpceaoutrec, "ERRTAI", total_abs_ih_err); 533 } 534 535 /** SAVNCPP, Example Function 12 (Requires Function 8) ***************/ 536 537 if (savncppflag) 538 { 539 // First, compute Jz on all pixels 540 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 541 derx, dery); 542 543 // Then compute sums of Jz on select pixels 544 if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err, 545 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 546 { 547 mbobra 1.1 totaljz = DRMS_MISSING_FLOAT; 548 totaljz_err = DRMS_MISSING_FLOAT; 549 } 550 551 // Set sharp keys 552 drms_setkey_float(sharpoutrec, "SAVNCPP", totaljz); 553 drms_setkey_float(sharpoutrec, "ERRJHT", totaljz_err); 554 555 // Set sharp cea keys 556 drms_setkey_float(sharpceaoutrec, "SAVNCPP", totaljz); 557 drms_setkey_float(sharpceaoutrec, "ERRJHT", totaljz_err); 558 } 559 560 /* MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) ******/ 561 562 if (meanshrflag \|\| shrgt45flag) 563 { 564 // First compute potential field 565 for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; 566 greenpot(bpx, bpy, bpz, nx, ny); 567 568 mbobra 1.1 // Then compute shear angles 569 if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, 570 &meanshear_angle, &meanshear_angle_err, &area_w_shear_gt_45, 571 mask, bitmask)) 572 { 573 meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN 574 area_w_shear_gt_45 = DRMS_MISSING_FLOAT; 575 meanshear_angle_err= DRMS_MISSING_FLOAT; 576 } 577 578 // Set sharp keys 579 drms_setkey_float(sharpoutrec, "MEANSHR", meanshear_angle); 580 drms_setkey_float(sharpoutrec, "SHRGT45", area_w_shear_gt_45); 581 drms_setkey_float(sharpoutrec, "ERRMSHA", meanshear_angle_err); 582 583 // Set sharp cea keys 584 drms_setkey_float(sharpceaoutrec, "MEANSHR", meanshear_angle); 585 drms_setkey_float(sharpceaoutrec, "SHRGT45", area_w_shear_gt_45); 586 drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err); 587 } 588 589 mbobra 1.1 /*************************** END FLAGS ********************************/ 590 591 drms_free_array(bitmaskArray); 592 drms_free_array(maskArray); 593 drms_free_array(bxArray); 594 drms_free_array(byArray); 595 drms_free_array(bzArray); 596 drms_free_array(bx_errArray); 597 drms_free_array(by_errArray); 598 drms_free_array(bz_errArray); 599 } //endfor 600 601 drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD); 602 drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD); 603 604 //used for select calculations 605 free(bh); free(bt); free(jz); 606 free(bpx); free(bpy); free(bpz); 607 free(derx); free(dery); 608 free(derx_bt); free(dery_bt); 609 free(derx_bz); free(dery_bz); 610 mbobra 1.1 free(derx_bh); free(dery_bh); 611 free(bt_err); free(bh_err); free(jz_err); 612 free(jz_err_squared); free(jz_rms_err); 613 free(cvsinfoall);
614 arta 1.5
615 arta 1.3 return 0;
616 mbobra 1.1 } // DoIt 617 618

Karen Tian