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 mbobra 1.9 * -- DEBUG flag (use like this debug=debug)
22 mbobra 1.1 * 23 * AUTHOR : Monica Bobra 24 * 25 * Version : v0.0 Jun 14 2013
26 mbobra 1.6 * v0.1 Feb 11 2014 Added different input and output series
27 mbobra 1.1 * 28 * EXAMPLE :
29 mbobra 1.6 * update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
30 mbobra 1.9 * HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
31 mbobra 1.1 * 32 / 33 34 / Include files / 35 #include <stdio.h> 36 #include <stdlib.h> 37 #include <time.h> 38 #include <sys/time.h> 39 #include <math.h> 40 #include <string.h> 41 #include "jsoc_main.h" 42 #include "astro.h" 43 #include "fstats.h" 44 #include "cartography.c" 45 #include "fresize.h" 46 #include "finterpolate.h" 47 #include "img2helioVector.c" 48 #include "copy_me_keys.c" 49 #include "errorprop.c" 50 #include "sw_functions.c" 51 #include <mkl_blas.h> 52 mbobra 1.1 #include <mkl_service.h> 53 #include <mkl_lapack.h> 54 #include <mkl_vml_functions.h> 55 #include <omp.h> 56 57 / Define values / 58 #define PI (M_PI) 59 #define RADSINDEG (PI/180.) 60 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0])) 61 #define DIE(msg) {fflush(stdout); fprintf(stderr,"%s, status=%d\n", msg, status); return(status);} 62 #define SHOW(msg) {printf("%s", msg); fflush(stdout);} 63 64 char module_name = "update_sharp_keys"; /* Module name */
65 mbobra 1.6 char version_id = "2014 Feb 12"; / Version number */
66 mbobra 1.1 67 ModuleArgs_t module_args[] = 68 {
69 mbobra 1.9 {ARG_STRING, "sharpseriesin", NULL, "Input Sharp dataseries"}, 70 {ARG_STRING, "sharpceaseriesin", NULL, "Input Sharp CEA dataseries"},
71 mbobra 1.6 {ARG_STRING, "sharpseriesout", NULL, "Output Sharp dataseries"}, 72 {ARG_STRING, "sharpceaseriesout", NULL, "Output Sharp CEA dataseries"},
73 mbobra 1.9 {ARG_INT, "HARPNUM", NULL, "HARP number"}, 74 {ARG_STRING, "keylist", NULL, "comma separated list of keywords to update"}, 75 {ARG_STRING, "debug", NULL, "debug mode functionality. use like this: debug=debug"},
76 mbobra 1.1 {ARG_END} 77 }; 78 79 80 int DoIt(void) 81 { 82 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ); 83 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ); 84 85 int status = DRMS_SUCCESS; 86 int nrecs, irec; 87 88 /* Keywords */ 89 float mean_vf; 90 float absFlux; 91 float count_mask; 92 float mean_hf; 93 float mean_gamma; 94 float mean_derivative_btotal; 95 float mean_derivative_bh; 96 float mean_derivative_bz; 97 mbobra 1.1 float mean_jz; 98 float us_i; 99 float mean_alpha; 100 float mean_ih; 101 float total_us_ih; 102 float total_abs_ih; 103 float totaljz; 104 float totpot; 105 float meanpot; 106 float area_w_shear_gt_45; 107 float meanshear_angle; 108 float area_w_shear_gt_45h; 109 float meanshear_angleh; 110 float mean_derivative_btotal_err; 111 float mean_vf_err; 112 float mean_gamma_err; 113 float mean_derivative_bh_err; 114 float mean_derivative_bz_err; 115 float mean_jz_err; 116 float us_i_err; 117 float mean_alpha_err; 118 mbobra 1.1 float mean_ih_err; 119 float total_us_ih_err; 120 float total_abs_ih_err; 121 float totaljz_err; 122 float meanpot_err; 123 float totpot_err;
124 mbobra 1.6 float Rparam;
125 mbobra 1.1 float meanshear_angle_err;
126 mbobra 1.8 float totfx; 127 float totfy; 128 float totfz; 129 float totbsq; 130 float epsx; 131 float epsy; 132 float epsz; 133
134 mbobra 1.6 char sharpseriesin = (char ) params_get_str(&cmdparams, "sharpseriesin"); 135 char sharpceaseriesin = (char ) params_get_str(&cmdparams, "sharpceaseriesin"); 136 char sharpseriesout = (char ) params_get_str(&cmdparams, "sharpseriesout"); 137 char sharpceaseriesout = (char ) params_get_str(&cmdparams, "sharpceaseriesout"); 138 139 int sameflag = strcmp(sharpseriesin, sharpseriesout) == 0; 140 int testflag = strcmp(sharpceaseriesin, sharpceaseriesout) == 0; 141 if (testflag != sameflag) 142 DIE("Either both outputs must be the same as their inputs, or both be different.");
143 mbobra 1.1 144 int harpnum = params_get_int(&cmdparams, "HARPNUM");
145 mbobra 1.6
146 mbobra 1.1 char sharpquery[256]; 147 char sharpceaquery[256];
148 mbobra 1.6 sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum); 149 printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
150 mbobra 1.1
151 mbobra 1.6 sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum); 152 sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
153 mbobra 1.1
154 mbobra 1.6 DRMS_RecordSet_t *sharpinrecset= drms_open_records(drms_env, sharpquery, &status);
155 mbobra 1.1 if (status != DRMS_SUCCESS) 156 DIE("Problem opening sharp recordset.");
157 mbobra 1.6 nrecs=sharpinrecset->n; 158 if (nrecs == 0)
159 mbobra 1.1 DIE("Sharp recordset contains no records.");
160 mbobra 1.6 printf("nrecs=%d\n",nrecs);
161 mbobra 1.1 162 DRMS_RecordSet_t *sharpceainrecset = drms_open_records(drms_env, sharpceaquery, &status); 163 if (status != DRMS_SUCCESS) 164 DIE("Problem opening sharp cea recordset.");
165 mbobra 1.6 if (sharpceainrecset->n != nrecs) 166 DIE("Sharp cea recordset differs from sharp recordset."); 167 168 DRMS_RecordSet_t sharpoutrecset, sharpceaoutrecset; 169 if (sameflag) 170 { 171 sharpoutrecset = drms_clone_records(sharpinrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status); 172 if (status != DRMS_SUCCESS) 173 DIE("Problem cloning sharp records."); 174 sharpceaoutrecset = drms_clone_records(sharpceainrecset, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
175 mbobra 1.1 if (status != DRMS_SUCCESS) 176 DIE("Problem cloning sharp cea records.");
177 mbobra 1.6 } 178 else 179 { 180 sharpoutrecset = drms_create_records(drms_env, nrecs, sharpseriesout, DRMS_PERMANENT, &status); 181 if (status != DRMS_SUCCESS) 182 DIE("Problem creating sharp records."); 183 sharpceaoutrecset = drms_create_records(drms_env, nrecs, sharpceaseriesout, DRMS_PERMANENT, &status); 184 if (status != DRMS_SUCCESS) 185 DIE("Problem creating sharp cea records."); 186 } 187
188 mbobra 1.1 189 char keylist = (char ) params_get_str(&cmdparams, "keylist");
190 mbobra 1.9 char debug = (char ) params_get_str(&cmdparams, "debug");
191 mbobra 1.1 192 // Flags to indicate which keyword will be recalculated 193 int meanvfflag = (strstr(keylist,"USFLUX") != NULL); // generalize so that lowercase is acceptable 194 int totpotflag = (strstr(keylist,"TOTPOT") != NULL); 195 int meangamflag = (strstr(keylist,"MEANGAM") != NULL); 196 int meangbtflag = (strstr(keylist,"MEANGBT") != NULL); 197 int meangbzflag = (strstr(keylist,"MEANGBZ") != NULL); 198 int meangbhflag = (strstr(keylist,"MEANGBH") != NULL); 199 int meanjzdflag = (strstr(keylist,"MEANJZD") != NULL); 200 int totusjzflag = (strstr(keylist,"TOTUSJZ") != NULL); 201 int meanalpflag = (strstr(keylist,"MEANALP") != NULL); 202 int meanjzhflag = (strstr(keylist,"MEANJZH") != NULL); 203 int totusjhflag = (strstr(keylist,"TOTUSJH") != NULL); 204 int absnjzhflag = (strstr(keylist,"ABSNJZH") != NULL); 205 int savncppflag = (strstr(keylist,"SAVNCPP") != NULL); 206 int meanpotflag = (strstr(keylist,"MEANPOT") != NULL); 207 int meanshrflag = (strstr(keylist,"MEANSHR") != NULL); 208 int shrgt45flag = (strstr(keylist,"SHRGT45") != NULL);
209 mbobra 1.6 int r_valueflag = (strstr(keylist,"R_VALUE") != NULL);
210 mbobra 1.8 int totbsqflag = (strstr(keylist,"TOTBSQ") != NULL); 211 int totfxflag = (strstr(keylist,"TOTFX") != NULL); 212 int totfyflag = (strstr(keylist,"TOTFY") != NULL); 213 int totfzflag = (strstr(keylist,"TOTFZ") != NULL); 214 int epsxflag = (strstr(keylist,"EPSX") != NULL); 215 int epsyflag = (strstr(keylist,"EPSY") != NULL); 216 int epszflag = (strstr(keylist,"EPSZ") != NULL);
217 mbobra 1.9 int debugflag = (strstr(debug,"debug") != NULL);
218 mbobra 1.8
219 mbobra 1.1 DRMS_Record_t sharpinrec = sharpinrecset->records[0]; 220 DRMS_Record_t sharpceainrec = sharpceainrecset->records[0]; 221 DRMS_Segment_t inseg = drms_segment_lookup(sharpceainrec, "Br"); 222 int nx = inseg->axis[0]; 223 int ny = inseg->axis[1]; 224 int nxny = nx ny; 225 int dims[2] = {nx, ny};
226 mbobra 1.6
227 mbobra 1.1 // Temp arrays 228 float bh = (float ) (malloc(nxny * sizeof(float))); 229 float bt = (float ) (malloc(nxny * sizeof(float))); 230 float jz = (float ) (malloc(nxny * sizeof(float))); 231 float bpx = (float ) (malloc(nxny * sizeof(float))); 232 float bpy = (float ) (malloc(nxny * sizeof(float))); 233 float bpz = (float ) (malloc(nxny * sizeof(float))); 234 float derx = (float ) (malloc(nxny * sizeof(float))); 235 float dery = (float ) (malloc(nxny * sizeof(float))); 236 float derx_bt = (float ) (malloc(nxny * sizeof(float))); 237 float dery_bt = (float ) (malloc(nxny * sizeof(float))); 238 float derx_bh = (float ) (malloc(nxny * sizeof(float))); 239 float dery_bh = (float ) (malloc(nxny * sizeof(float))); 240 float derx_bz = (float ) (malloc(nxny * sizeof(float))); 241 float dery_bz = (float ) (malloc(nxny * sizeof(float))); 242 float bt_err = (float ) (malloc(nxny * sizeof(float))); 243 float bh_err = (float ) (malloc(nxny * sizeof(float))); 244 float jz_err = (float ) (malloc(nxny * sizeof(float))); 245 float jz_err_squared = (float ) (malloc(nxny * sizeof(float))); 246 float jz_rms_err = (float ) (malloc(nxny * sizeof(float))); 247 float jz_err_squared_smooth = (float ) (malloc(nxny * sizeof(float))); 248 mbobra 1.1 float jz_smooth = (float ) (malloc(nxny * sizeof(float))); 249 250 // ephemeris variables 251 float cdelt1_orig, cdelt1, dsun_obs, imcrpix1, imcrpix2, crpix1, crpix2; 252 double rsun_ref, rsun_obs; 253 254 // outrecords
255 mbobra 1.6 DRMS_Record_t sharpoutrec, sharpceaoutrec, *testrec;
256 mbobra 1.1 257 // prepare to set CODEVER7 (CVS Version of the SHARP module) 258 char cvsinfo0; 259 char history0;
260 mbobra 1.9 char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.8 2014/06/05 21:27:32 mbobra Exp $");
261 mbobra 1.1 char cvsinfo2 = sw_functions_version(); 262 char cvsinfoall = (char *)malloc(2048); 263 char historyofthemodule[2048]; 264 cvsinfo0 = drms_getkey_string(sharpinrec, "CODEVER7", &status); 265 strcpy(cvsinfoall,cvsinfo0); 266 strcat(cvsinfoall,"\n"); 267 strcat(cvsinfoall,cvsinfo2); 268 269 // prepare to set HISTORY (History of the data) 270 char timebuf[1024]; 271 float UNIX_epoch = -220924792.000; 272 sprint_time(timebuf, (double)time(NULL) + UNIX_epoch, "ISO", 0); 273 sprintf(historyofthemodule,"The following keywords were re-computed on %s: %s.",timebuf,keylist); 274 printf("historyofthemodule=%s\n",historyofthemodule); 275 history0 = drms_getkey_string(sharpinrec, "HISTORY", &status); 276 strcat(historyofthemodule,"\n"); 277 strcat(historyofthemodule,history0); 278
279 mbobra 1.6 // grab one record such that we can malloc the R-parameter calculation arrays outside of the loop 280 // assume scale = round(2./cdelt1) does not change for any given HARP 281 testrec = sharpceainrecset->records[nrecs/2]; 282 dsun_obs = drms_getkey_float(testrec, "DSUN_OBS", &status); 283 rsun_ref = drms_getkey_double(testrec, "RSUN_REF", &status); 284 cdelt1_orig = drms_getkey_float(testrec, "CDELT1", &status); 285 cdelt1 = (atan((rsun_refcdelt1_origRADSINDEG)/(dsun_obs)))(1/RADSINDEG)(3600.); 286 int scale = round(2.0/cdelt1); 287 int nx1 = nx/scale; 288 int ny1 = ny/scale;
289 mbobra 1.7 int nxp = nx1+40; 290 int nyp = ny1+40;
291 mbobra 1.6 292 if (nx1 > floor((nx-1)/scale + 1) ) 293 DIE("X-dimension of output array in fsample() is too large."); 294 if (ny1 > floor((ny-1)/scale + 1) ) 295 DIE("Y-dimension of output array in fsample() is too large."); 296 297 // malloc some arrays for the R parameter calculation
298 mbobra 1.7 float rim = (float )malloc(nx1ny1sizeof(float)); 299 float p1p0 = (float )malloc(nx1ny1sizeof(float)); 300 float p1n0 = (float )malloc(nx1ny1sizeof(float)); 301 float p1p = (float )malloc(nx1ny1sizeof(float)); 302 float p1n = (float )malloc(nx1ny1sizeof(float)); 303 float p1 = (float )malloc(nx1ny1sizeof(float)); 304 float pmap = (float )malloc(nxpnypsizeof(float)); 305 float p1pad = (float )malloc(nxpnypsizeof(float)); 306 float pmapn = (float )malloc(nx1ny1sizeof(float));
307 mbobra 1.8 308 // define some arrays for the lorentz force calculation 309 float fx = (float ) (malloc(nxny * sizeof(float))); 310 float fy = (float ) (malloc(nxny * sizeof(float))); 311 float fz = (float ) (malloc(nxny * sizeof(float))); 312
313 mbobra 1.7
314 mbobra 1.1 for (irec=0;irec<nrecs;irec++) 315 { 316 // Get emphemeris
317 mbobra 1.7 sharpinrec = sharpinrecset->records[irec];
318 mbobra 1.6 sharpceainrec = sharpceainrecset->records[irec];
319 mbobra 1.7 cdelt1_orig = drms_getkey_float(sharpceainrec, "CDELT1", &status); 320 dsun_obs = drms_getkey_float(sharpinrec, "DSUN_OBS", &status); 321 rsun_ref = drms_getkey_double(sharpinrec, "RSUN_REF", &status); 322 rsun_obs = drms_getkey_double(sharpinrec, "RSUN_OBS", &status); 323 imcrpix1 = drms_getkey_float(sharpinrec, "IMCRPIX1", &status); 324 imcrpix2 = drms_getkey_float(sharpinrec, "IMCRPIX2", &status); 325 crpix1 = drms_getkey_float(sharpinrec, "CRPIX1", &status); 326 crpix2 = drms_getkey_float(sharpinrec, "CRPIX2", &status);
327 mbobra 1.6
328 mbobra 1.1 sharpoutrec = sharpoutrecset->records[irec]; 329 sharpceaoutrec = sharpceaoutrecset->records[irec];
330 mbobra 1.6 if (!sameflag) 331 { 332 drms_copykeys(sharpoutrec, sharpinrec, 0, kDRMS_KeyClass_Explicit); 333 drms_copykeys(sharpceaoutrec, sharpceainrec, 0, kDRMS_KeyClass_Explicit); 334 } 335 336 TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */ 337 tnow = (double)time(NULL); 338 tnow += UNIX_epoch;
339 mbobra 1.1
340 mbobra 1.9 if (debugflag) 341 { 342 printf("i=%d\n",irec); 343 } 344
345 mbobra 1.6 // set CODEVER7 and HISTORY and DATE
346 mbobra 1.1 drms_setkey_string(sharpoutrec, "CODEVER7", cvsinfoall); 347 drms_setkey_string(sharpoutrec,"HISTORY",historyofthemodule); 348 drms_setkey_string(sharpceaoutrec, "CODEVER7", cvsinfoall); 349 drms_setkey_string(sharpceaoutrec,"HISTORY",historyofthemodule);
350 mbobra 1.6 drms_setkey_time(sharpoutrec,"DATE",tnow); 351 drms_setkey_time(sharpceaoutrec,"DATE",tnow);
352 mbobra 1.1 353 // Get bx, by, bz, mask 354 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
355 mbobra 1.6 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
356 mbobra 1.1 DRMS_Array_t bitmaskArray = drms_segment_read(bitmaskSeg, DRMS_TYPE_INT, &status); 357 int bitmask = (int *) bitmaskArray->data; // get the previously made mask array 358 359 //Use conf_disambig map as a threshold on spaceweather quantities
360 mbobra 1.6 DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");
361 mbobra 1.1 DRMS_Array_t maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status); 362 if (status != DRMS_SUCCESS) 363 DIE("No CONF_DISAMBIG segment."); 364 int mask = (int *) maskArray->data; // get the previously made mask array 365
366 mbobra 1.6 //Use magnetogram map to compute R 367 DRMS_Segment_t losSeg = drms_segment_lookup(sharpceainrec, "magnetogram"); 368 DRMS_Array_t losArray = drms_segment_read(losSeg, DRMS_TYPE_FLOAT, &status); 369 float los = (float ) losArray->data; // los 370 371 DRMS_Segment_t *bxSeg = drms_segment_lookup(sharpceainrec, "Bp");
372 mbobra 1.1 DRMS_Array_t bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status); 373 float bx = (float *) bxArray->data; // bx 374
375 mbobra 1.6 DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
376 mbobra 1.1 DRMS_Array_t byArray = drms_segment_read(bySeg, DRMS_TYPE_FLOAT, &status); 377 float by = (float ) byArray->data; // by 378 for (int i = 0; i < nxny; i++) by[i] = -1; 379
380 mbobra 1.6 DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
381 mbobra 1.1 DRMS_Array_t bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status); 382 float bz = (float *) bzArray->data; // bz 383
384 mbobra 1.6 DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
385 mbobra 1.1 DRMS_Array_t bz_errArray = drms_segment_read(bz_errSeg, DRMS_TYPE_FLOAT, &status); 386 float bz_err = (float *) bz_errArray->data; // bz_err 387
388 mbobra 1.6 DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
389 mbobra 1.1 DRMS_Array_t by_errArray = drms_segment_read(by_errSeg, DRMS_TYPE_FLOAT, &status); 390 float by_err = (float *) by_errArray->data; // by_err 391
392 mbobra 1.6 DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
393 mbobra 1.1 DRMS_Array_t bx_errArray = drms_segment_read(bx_errSeg, DRMS_TYPE_FLOAT, &status); 394 float bx_err = (float *) bx_errArray->data; // bx_err 395
396 mbobra 1.6 /*** USFLUX, Example Function 1 ***********************************/
397 mbobra 1.1 if (meanvfflag) 398 { 399 // Compute unsigned flux 400 if (computeAbsFlux(bz_err, bz , dims, &absFlux, &mean_vf, &mean_vf_err, 401 &count_mask, mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 402 { 403 mean_vf = DRMS_MISSING_FLOAT; 404 mean_vf_err = DRMS_MISSING_FLOAT; 405 count_mask = DRMS_MISSING_INT; 406 } 407 408 drms_setkey_float(sharpoutrec, "USFLUX", mean_vf); 409 drms_setkey_float(sharpoutrec, "CMASK", count_mask); 410 drms_setkey_float(sharpoutrec, "ERRVF", mean_vf_err); 411 412 drms_setkey_float(sharpceaoutrec, "USFLUX", mean_vf); 413 drms_setkey_float(sharpceaoutrec, "CMASK", count_mask); 414 drms_setkey_float(sharpceaoutrec, "ERRVF", mean_vf_err); 415 }
416 mbobra 1.6 417 /*** RPARAM, Example Function 14 ***********************************/ 418 if (r_valueflag) 419 { 420 if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
421 mbobra 1.7 p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn)) 422
423 mbobra 1.6 { 424 Rparam = DRMS_MISSING_FLOAT; 425 } 426 427 drms_setkey_float(sharpoutrec, "R_VALUE", Rparam); 428 drms_setkey_float(sharpceaoutrec, "R_VALUE", Rparam); 429 430 } 431
432 mbobra 1.1 /*** MEANPOT and TOTPOT, Example Function 13 (Requires Keiji's code) / 433 434 if (totpotflag \|\| meanpotflag) 435 { 436 // First compute potential field 437 for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; 438 greenpot(bpx, bpy, bpz, nx, ny); 439 440 // Then compute energy 441 if (computeFreeEnergy(bx_err, by_err, bx, by, bpx, bpy, dims, 442 &meanpot, &meanpot_err, &totpot, &totpot_err, 443 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 444 { 445 meanpot = DRMS_MISSING_FLOAT; 446 totpot = DRMS_MISSING_FLOAT; 447 meanpot_err = DRMS_MISSING_FLOAT; 448 totpot_err = DRMS_MISSING_FLOAT; 449 } 450 451 // Set sharp keys 452 drms_setkey_float(sharpoutrec, "MEANPOT", meanpot); 453 mbobra 1.1 drms_setkey_float(sharpoutrec, "TOTPOT", totpot); 454 drms_setkey_float(sharpoutrec, "ERRMPOT", meanpot_err); 455 drms_setkey_float(sharpoutrec, "ERRTPOT", totpot_err); 456 457 // Set sharp cea keys 458 drms_setkey_float(sharpceaoutrec, "MEANPOT", meanpot); 459 drms_setkey_float(sharpceaoutrec, "TOTPOT", totpot); 460 drms_setkey_float(sharpceaoutrec, "ERRMPOT", meanpot_err); 461 drms_setkey_float(sharpceaoutrec, "ERRTPOT", totpot_err); 462 } 463 464 /*** MEANGAM, Example Function 3 ********************************/ 465 466 if (meangamflag) 467 { 468 // First compute horizontal field 469 computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask); 470 471 // Then compute inclination angle, gamma 472 if (computeGamma(bz_err, bh_err, bx, by, bz, bh, dims, &mean_gamma, &mean_gamma_err,mask, bitmask)) 473 { 474 mbobra 1.1 mean_gamma = DRMS_MISSING_FLOAT; 475 mean_gamma_err = DRMS_MISSING_FLOAT; 476 } 477 478 // Set sharp keys 479 drms_setkey_float(sharpoutrec, "MEANGAM", mean_gamma); 480 drms_setkey_float(sharpoutrec, "ERRGAM", mean_gamma_err); 481 482 // Set sharp cea keys 483 drms_setkey_float(sharpceaoutrec, "MEANGAM", mean_gamma); 484 drms_setkey_float(sharpceaoutrec, "ERRGAM", mean_gamma_err); 485 } 486 487 /* MEANGBT, Example Function 5 (Requires Function 4) *********/ 488 489 if (meangbtflag) 490 { 491 // First, compute Bt 492 computeB_total(bx_err, by_err, bz_err, bt_err, bx, by, bz, bt, dims, mask, bitmask); 493 494 // Then take the derivative of Bt 495 mbobra 1.1 if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask, bitmask, derx_bt, dery_bt, bt_err, 496 &mean_derivative_btotal_err)) 497 { 498 mean_derivative_btotal = DRMS_MISSING_FLOAT; 499 mean_derivative_btotal_err = DRMS_MISSING_FLOAT; 500 } 501 502 // Set sharp keys 503 drms_setkey_float(sharpoutrec, "MEANGBT", mean_derivative_btotal); 504 drms_setkey_float(sharpoutrec, "ERRBT", mean_derivative_btotal_err); 505 506 // Set sharp cea keys 507 drms_setkey_float(sharpceaoutrec, "MEANGBT", mean_derivative_btotal); 508 drms_setkey_float(sharpceaoutrec, "ERRBT", mean_derivative_btotal_err); 509 } 510 511 /* MEANGBH, Example Function 6 (Requires Function 2) *********/ 512 513 if (meangbhflag) 514 { 515 // First, compute Bh 516 mbobra 1.1 computeBh(bx_err, by_err, bh_err, bx, by, bz, bh, dims, &mean_hf, mask, bitmask); 517 518 // Then take the derivative of Bh 519 if (computeBhderivative(bh, bh_err, dims, &mean_derivative_bh, &mean_derivative_bh_err, mask, bitmask, derx_bh, dery_bh)) 520 { 521 mean_derivative_bh = DRMS_MISSING_FLOAT; 522 mean_derivative_bh_err = DRMS_MISSING_FLOAT; 523 } 524 525 // Set sharp keys 526 drms_setkey_float(sharpoutrec, "MEANGBH", mean_derivative_bh); 527 drms_setkey_float(sharpoutrec, "ERRBH", mean_derivative_bh_err); 528 529 // Set sharp cea keys 530 drms_setkey_float(sharpceaoutrec, "MEANGBH", mean_derivative_bh); 531 drms_setkey_float(sharpceaoutrec, "ERRBH", mean_derivative_bh_err); 532 } 533 534 /* MEANGBZ, Example Function 7 ********************************/ 535 536 if (meangbzflag) 537 mbobra 1.1 { 538 // Compute Bz derivative 539 if (computeBzderivative(bz, bz_err, dims, &mean_derivative_bz, &mean_derivative_bz_err, mask, bitmask, derx_bz, dery_bz)) 540 { 541 mean_derivative_bz = DRMS_MISSING_FLOAT; 542 mean_derivative_bz_err = DRMS_MISSING_FLOAT; 543 } 544 545 // Set sharp keys 546 drms_setkey_float(sharpoutrec, "MEANGBZ", mean_derivative_bz); 547 drms_setkey_float(sharpoutrec, "ERRBZ", mean_derivative_bz_err); 548 549 // Set sharp cea keys 550 drms_setkey_float(sharpceaoutrec, "MEANGBZ", mean_derivative_bz); 551 drms_setkey_float(sharpceaoutrec, "ERRBZ", mean_derivative_bz_err); 552 } 553 554 /* MEANJZD and TOTUSJZ, Example Function 9 (Requires Function 8) / 555 556 if (meanjzdflag \|\| totusjzflag) 557 { 558 mbobra 1.1 // First, compute Jz on all pixels 559 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 560 derx, dery); 561 562 // Then take the sums of the appropriate pixels 563 if(computeJzsmooth(bx, by, dims, jz, jz_smooth, jz_err, jz_rms_err, jz_err_squared_smooth, &mean_jz, 564 &mean_jz_err, &us_i, &us_i_err, mask, bitmask, cdelt1, 565 rsun_ref, rsun_obs, derx, dery)) 566 567 568 { 569 mean_jz = DRMS_MISSING_FLOAT; 570 us_i = DRMS_MISSING_FLOAT; 571 mean_jz_err = DRMS_MISSING_FLOAT; 572 us_i_err = DRMS_MISSING_FLOAT; 573 } 574 575 // Set sharp keys 576 drms_setkey_float(sharpoutrec, "MEANJZD", mean_jz); 577 drms_setkey_float(sharpoutrec, "TOTUSJZ", us_i); 578 drms_setkey_float(sharpoutrec, "ERRJZ", mean_jz_err); 579 mbobra 1.1 drms_setkey_float(sharpoutrec, "ERRUSI", us_i_err); 580 581 // Set sharp cea keys 582 drms_setkey_float(sharpceaoutrec, "MEANJZD", mean_jz); 583 drms_setkey_float(sharpceaoutrec, "TOTUSJZ", us_i); 584 drms_setkey_float(sharpceaoutrec, "ERRJZ", mean_jz_err); 585 drms_setkey_float(sharpceaoutrec, "ERRUSI", us_i_err); 586 } 587 588 /** MEANALP, Example Function 10 (Requires Function 8)*****/ 589 590 if (meanalpflag) 591 { 592 // First, compute Jz on all pixels 593 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 594 derx, dery); 595 596 // Then compute alpha quantities on select pixels 597 if (computeAlpha(jz_err, bz_err, bz, dims, jz, jz_smooth, &mean_alpha, &mean_alpha_err, 598 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 599 { 600 mbobra 1.1 mean_alpha = DRMS_MISSING_FLOAT; 601 mean_alpha_err = DRMS_MISSING_FLOAT; 602 } 603 604 // Set sharp keys 605 drms_setkey_float(sharpoutrec, "MEANALP", mean_alpha); 606 drms_setkey_float(sharpoutrec, "ERRALP", mean_alpha_err); 607 608 // Set sharp cea keys 609 drms_setkey_float(sharpceaoutrec, "MEANALP", mean_alpha); 610 drms_setkey_float(sharpceaoutrec, "ERRALP", mean_alpha_err); 611 612 } 613 614 /* MEANJZH, TOTUSJH, ABSNJZH, Example Function 11 (Requires Function 8) / 615 616 if (meanjzhflag \|\| totusjhflag \|\| absnjzhflag) 617 { 618 // First, compute Jz on all pixels 619 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 620 derx, dery); 621 mbobra 1.1 622 // Then compute helicity quantities on select pixels 623 if (computeHelicity(jz_err, jz_rms_err, bz_err, bz, dims, jz, &mean_ih, &mean_ih_err, &total_us_ih, 624 &total_abs_ih, &total_us_ih_err, &total_abs_ih_err, mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 625 { 626 mean_ih = DRMS_MISSING_FLOAT; 627 total_us_ih = DRMS_MISSING_FLOAT; 628 total_abs_ih = DRMS_MISSING_FLOAT; 629 mean_ih_err = DRMS_MISSING_FLOAT; 630 total_us_ih_err = DRMS_MISSING_FLOAT; 631 total_abs_ih_err = DRMS_MISSING_FLOAT; 632 } 633 634 // Set sharp keys 635 drms_setkey_float(sharpoutrec, "MEANJZH", mean_ih); 636 drms_setkey_float(sharpoutrec, "TOTUSJH", total_us_ih); 637 drms_setkey_float(sharpoutrec, "ABSNJZH", total_abs_ih); 638 drms_setkey_float(sharpoutrec, "ERRMIH", mean_ih_err); 639 drms_setkey_float(sharpoutrec, "ERRTUI", total_us_ih_err); 640 drms_setkey_float(sharpoutrec, "ERRTAI", total_abs_ih_err); 641 642 mbobra 1.1 // Set sharp cea keys 643 drms_setkey_float(sharpceaoutrec, "MEANJZH", mean_ih); 644 drms_setkey_float(sharpceaoutrec, "TOTUSJH", total_us_ih); 645 drms_setkey_float(sharpceaoutrec, "ABSNJZH", total_abs_ih); 646 drms_setkey_float(sharpceaoutrec, "ERRMIH", mean_ih_err); 647 drms_setkey_float(sharpceaoutrec, "ERRTUI", total_us_ih_err); 648 drms_setkey_float(sharpceaoutrec, "ERRTAI", total_abs_ih_err); 649 } 650 651 /** SAVNCPP, Example Function 12 (Requires Function 8) ***************/ 652 653 if (savncppflag) 654 { 655 // First, compute Jz on all pixels 656 computeJz(bx_err, by_err, bx, by, dims, jz, jz_err, jz_err_squared, mask, bitmask, cdelt1, rsun_ref, rsun_obs, 657 derx, dery); 658 659 // Then compute sums of Jz on select pixels 660 if (computeSumAbsPerPolarity(jz_err, bz_err, bz, jz, dims, &totaljz, &totaljz_err, 661 mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 662 { 663 mbobra 1.1 totaljz = DRMS_MISSING_FLOAT; 664 totaljz_err = DRMS_MISSING_FLOAT; 665 } 666 667 // Set sharp keys 668 drms_setkey_float(sharpoutrec, "SAVNCPP", totaljz); 669 drms_setkey_float(sharpoutrec, "ERRJHT", totaljz_err); 670 671 // Set sharp cea keys 672 drms_setkey_float(sharpceaoutrec, "SAVNCPP", totaljz); 673 drms_setkey_float(sharpceaoutrec, "ERRJHT", totaljz_err); 674 } 675 676 /* MEANSHR and SHRGT45, Example Function 14 (Requires Keiji's code) ********/ 677 678 if (meanshrflag \|\| shrgt45flag) 679 { 680 // First compute potential field 681 for (int i = 0; i < nxny; i++) bpz[i] = bz[i]; 682 greenpot(bpx, bpy, bpz, nx, ny); 683 684 mbobra 1.1 // Then compute shear angles 685 if (computeShearAngle(bx_err, by_err, bh_err, bx, by, bz, bpx, bpy, bpz, dims, 686 &meanshear_angle, &meanshear_angle_err, &area_w_shear_gt_45, 687 mask, bitmask)) 688 { 689 meanshear_angle = DRMS_MISSING_FLOAT; // If fail, fill in NaN 690 area_w_shear_gt_45 = DRMS_MISSING_FLOAT; 691 meanshear_angle_err= DRMS_MISSING_FLOAT; 692 } 693 694 // Set sharp keys 695 drms_setkey_float(sharpoutrec, "MEANSHR", meanshear_angle); 696 drms_setkey_float(sharpoutrec, "SHRGT45", area_w_shear_gt_45); 697 drms_setkey_float(sharpoutrec, "ERRMSHA", meanshear_angle_err); 698 699 // Set sharp cea keys 700 drms_setkey_float(sharpceaoutrec, "MEANSHR", meanshear_angle); 701 drms_setkey_float(sharpceaoutrec, "SHRGT45", area_w_shear_gt_45); 702 drms_setkey_float(sharpceaoutrec, "ERRMSHA", meanshear_angle_err); 703 } 704
705 mbobra 1.8 /*** TOTFX, TOTFY, TOTFX, TOTBSQ, EPSX, EPSY, EPSZ , Example Function 16 (Lorentz forces) ********/ 706 707 if (totfxflag \|\| totfyflag \|\| totfzflag \|\| totbsqflag \|\| epsxflag \|\| epsyflag \|\| epszflag) 708 { 709 710 // Compute Lorentz forces 711 if (computeLorentz(bx, by, bz, fx, fy, fz, dims, &totfx, &totfy, &totfz, &totbsq, 712 &epsx, &epsy, &epsz, mask, bitmask, cdelt1, rsun_ref, rsun_obs)) 713 { 714 totfx = DRMS_MISSING_FLOAT; 715 totfy = DRMS_MISSING_FLOAT; 716 totfz = DRMS_MISSING_FLOAT; 717 totbsq = DRMS_MISSING_FLOAT; 718 epsx = DRMS_MISSING_FLOAT; 719 epsy = DRMS_MISSING_FLOAT; 720 epsz = DRMS_MISSING_FLOAT; 721 } 722 723 // Set sharp keys 724 drms_setkey_float(sharpoutrec, "TOTFX", totfx); 725 drms_setkey_float(sharpoutrec, "TOTFY", totfy); 726 mbobra 1.8 drms_setkey_float(sharpoutrec, "TOTFZ", totfz); 727 drms_setkey_float(sharpoutrec, "TOTBSQ",totbsq); 728 drms_setkey_float(sharpoutrec, "EPSX", epsx); 729 drms_setkey_float(sharpoutrec, "EPSY", epsy); 730 drms_setkey_float(sharpoutrec, "EPSZ", epsz); 731 732 // Set sharp cea keys 733 drms_setkey_float(sharpceaoutrec, "TOTFX", totfx); 734 drms_setkey_float(sharpceaoutrec, "TOTFY", totfy); 735 drms_setkey_float(sharpceaoutrec, "TOTFZ", totfz); 736 drms_setkey_float(sharpceaoutrec, "TOTBSQ",totbsq); 737 drms_setkey_float(sharpceaoutrec, "EPSX", epsx); 738 drms_setkey_float(sharpceaoutrec, "EPSY", epsy); 739 drms_setkey_float(sharpceaoutrec, "EPSZ", epsz); 740 741 } 742
743 mbobra 1.1 /***************************** END FLAGS ********************************/ 744 745 drms_free_array(bitmaskArray); 746 drms_free_array(maskArray); 747 drms_free_array(bxArray); 748 drms_free_array(byArray); 749 drms_free_array(bzArray); 750 drms_free_array(bx_errArray); 751 drms_free_array(by_errArray); 752 drms_free_array(bz_errArray);
753 mbobra 1.6 drms_free_array(losArray);
754 mbobra 1.1 } //endfor 755
756 mbobra 1.6 drms_close_records(sharpinrecset, DRMS_FREE_RECORD); 757 drms_close_records(sharpceainrecset, DRMS_FREE_RECORD); 758
759 mbobra 1.1 drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD); 760 drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD); 761 762 //used for select calculations 763 free(bh); free(bt); free(jz); 764 free(bpx); free(bpy); free(bpz); 765 free(derx); free(dery); 766 free(derx_bt); free(dery_bt); 767 free(derx_bz); free(dery_bz); 768 free(derx_bh); free(dery_bh); 769 free(bt_err); free(bh_err); free(jz_err); 770 free(jz_err_squared); free(jz_rms_err); 771 free(cvsinfoall);
772 mbobra 1.7
773 mbobra 1.6 free(rim); 774 free(p1p0); 775 free(p1n0); 776 free(p1p); 777 free(p1n); 778 free(p1); 779 free(pmap);
780 mbobra 1.7 free(p1pad); 781 free(pmapn);
782 mbobra 1.6
783 mbobra 1.8 // free the arrays that are related to the lorentz calculation 784 free(fx); free(fy); free(fz);
785 mbobra 1.9 printf("Success=%d\n",sameflag);
786 arta 1.3 return 0;
787 mbobra 1.1 } // DoIt

Karen Tian