ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/JSOC/proj/sharp/apps/update_sharp_keys.c
Revision: 1.10
Committed: Mon Jun 9 21:07:41 2014 UTC (9 years, 3 months ago) by mbobra
Content type: text/plain
Branch: MAIN
CVS Tags: Ver_8-6, Ver_8-5
Changes since 1.9: +26 -4 lines
Log Message: 
added check on CEA series array dimensions
if nx by ny are not constant for any given CEA harpnum, the program will exit

File Contents

# Content
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 * 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 * 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 * -- DEBUG flag (use like this debug=debug)
28 *
29 * AUTHOR : Monica Bobra
30 *
31 * Version : v0.0 Jun 14 2013
32 * v0.1 Feb 11 2014 Added different input and output series
33 *
34 * EXAMPLE :
35 * update_sharp_keys sharpseriesin=hmi.sharp_720s sharpceaseriesin=hmi.sharp_cea_720s //
36 * HARPNUM=1 sharpseriesout=hmi.sharp_720s sharpceaseriesout=hmi.sharp_cea_720s keylist=USFLUX,TOTPOT debug=debug
37 *
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 #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 char *version_id = "2014 Feb 12"; /* Version number */
72
73 ModuleArgs_t module_args[] =
74 {
75 {ARG_STRING, "sharpseriesin", NULL, "Input Sharp dataseries"},
76 {ARG_STRING, "sharpceaseriesin", NULL, "Input Sharp CEA dataseries"},
77 {ARG_STRING, "sharpseriesout", NULL, "Output Sharp dataseries"},
78 {ARG_STRING, "sharpceaseriesout", NULL, "Output Sharp CEA dataseries"},
79 {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 {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 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 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 float Rparam;
131 float meanshear_angle_err;
132 float totfx;
133 float totfy;
134 float totfz;
135 float totbsq;
136 float epsx;
137 float epsy;
138 float epsz;
139
140 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
150 int harpnum = params_get_int(&cmdparams, "HARPNUM");
151
152 char sharpquery[256];
153 char sharpceaquery[256];
154 sprintf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
155 printf(sharpceaquery, "%s[%d]\n", sharpceaseriesin,harpnum);
156
157 sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
158 sprintf(sharpquery, "%s[%d]\n", sharpseriesin,harpnum);
159
160 DRMS_RecordSet_t *sharpinrecset= drms_open_records(drms_env, sharpquery, &status);
161 if (status != DRMS_SUCCESS)
162 DIE("Problem opening sharp recordset.");
163 nrecs=sharpinrecset->n;
164 if (nrecs == 0)
165 DIE("Sharp recordset contains no records.");
166 printf("nrecs=%d\n",nrecs);
167
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 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 if (status != DRMS_SUCCESS)
182 DIE("Problem cloning sharp cea records.");
183 }
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
195 char *keylist = (char *) params_get_str(&cmdparams, "keylist");
196 char *debug = (char *) params_get_str(&cmdparams, "debug");
197
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 int r_valueflag = (strstr(keylist,"R_VALUE") != NULL);
216 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 int debugflag = (strstr(debug,"debug") != NULL);
224
225 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 //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
237 // 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 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 DRMS_Record_t *sharpoutrec, *sharpceaoutrec, *testrec;
266
267 // prepare to set CODEVER7 (CVS Version of the SHARP module)
268 char *cvsinfo0;
269 char *history0;
270 char *cvsinfo1 = strdup("$Id: update_sharp_keys.c,v 1.9 2014/06/06 19:04:51 mbobra Exp $");
271 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 // 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_ref*cdelt1_orig*RADSINDEG)/(dsun_obs)))*(1/RADSINDEG)*(3600.);
296 int scale = round(2.0/cdelt1);
297 int nx1 = nx/scale;
298 int ny1 = ny/scale;
299 int nxp = nx1+40;
300 int nyp = ny1+40;
301
302 // check to make sure that the dimensions passed into fsample() are adequate
303 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 DIE("Y-dimension of output array in fsample() is too large.");
307
308 // malloc some arrays for the R parameter calculation
309 float *rim = (float *)malloc(nx1*ny1*sizeof(float));
310 float *p1p0 = (float *)malloc(nx1*ny1*sizeof(float));
311 float *p1n0 = (float *)malloc(nx1*ny1*sizeof(float));
312 float *p1p = (float *)malloc(nx1*ny1*sizeof(float));
313 float *p1n = (float *)malloc(nx1*ny1*sizeof(float));
314 float *p1 = (float *)malloc(nx1*ny1*sizeof(float));
315 float *pmap = (float *)malloc(nxp*nyp*sizeof(float));
316 float *p1pad = (float *)malloc(nxp*nyp*sizeof(float));
317 float *pmapn = (float *)malloc(nx1*ny1*sizeof(float));
318
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 for (irec=0;irec<nrecs;irec++)
325 {
326
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 // Get emphemeris
337 sharpinrec = sharpinrecset->records[irec];
338 sharpceainrec = sharpceainrecset->records[irec];
339 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
348 sharpoutrec = sharpoutrecset->records[irec];
349 sharpceaoutrec = sharpceaoutrecset->records[irec];
350 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
360 if (debugflag)
361 {
362 printf("i=%d\n",irec);
363 }
364
365 // set CODEVER7 and HISTORY and DATE
366 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 drms_setkey_time(sharpoutrec,"DATE",tnow);
371 drms_setkey_time(sharpceaoutrec,"DATE",tnow);
372
373 // Get bx, by, bz, mask
374 // Use HARP (Turmon) bitmap as a threshold on spaceweather quantities
375 DRMS_Segment_t *bitmaskSeg = drms_segment_lookup(sharpceainrec, "bitmap");
376 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 DRMS_Segment_t *maskSeg = drms_segment_lookup(sharpceainrec, "conf_disambig");
381 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 //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 DRMS_Array_t *bxArray = drms_segment_read(bxSeg, DRMS_TYPE_FLOAT, &status);
393 float *bx = (float *) bxArray->data; // bx
394
395 DRMS_Segment_t *bySeg = drms_segment_lookup(sharpceainrec, "Bt");
396 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 DRMS_Segment_t *bzSeg = drms_segment_lookup(sharpceainrec, "Br");
401 DRMS_Array_t *bzArray = drms_segment_read(bzSeg, DRMS_TYPE_FLOAT, &status);
402 float *bz = (float *) bzArray->data; // bz
403
404 DRMS_Segment_t *bz_errSeg = drms_segment_lookup(sharpceainrec, "Br_err");
405 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 DRMS_Segment_t *by_errSeg = drms_segment_lookup(sharpceainrec, "Bt_err");
409 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 DRMS_Segment_t *bx_errSeg = drms_segment_lookup(sharpceainrec, "Bp_err");
413 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 /***** USFLUX, Example Function 1 *************************************/
417 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
437 /***** RPARAM, Example Function 14 *************************************/
438 if (r_valueflag)
439 {
440 if (computeR(bz_err, los , dims, &Rparam, cdelt1, rim, p1p0, p1n0,
441 p1p, p1n, p1, pmap, nx1, ny1, scale, p1pad, nxp, nyp, pmapn))
442
443 {
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 /***** 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 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 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 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 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 {
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 // 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 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 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
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 // 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 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 // 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 /***** 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 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 /******************************* END FLAGS **********************************/
764
765
766 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 drms_free_array(losArray);
775
776 } //endfor
777
778 drms_close_records(sharpinrecset, DRMS_FREE_RECORD);
779 drms_close_records(sharpceainrecset, DRMS_FREE_RECORD);
780
781 drms_close_records(sharpoutrecset, DRMS_INSERT_RECORD);
782 drms_close_records(sharpceaoutrecset, DRMS_INSERT_RECORD);
783
784 //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
795 free(rim);
796 free(p1p0);
797 free(p1n0);
798 free(p1p);
799 free(p1n);
800 free(p1);
801 free(pmap);
802 free(p1pad);
803 free(pmapn);
804
805 // free the arrays that are related to the lorentz calculation
806 free(fx); free(fy); free(fz);
807 printf("Success=%d\n",sameflag);
808 return 0;
809 } // DoIt