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

Karen Tian