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