ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/JSOC/proj/sharp/apps/update_sharp_keys.c
Revision: 1.19
Committed: Wed May 26 19:25:17 2021 UTC (2 years, 3 months ago) by mbobra
Content type: text/plain
Branch: MAIN
CVS Tags: Ver_LATEST, Ver_9-5, HEAD
Changes since 1.18: +88 -96 lines
Log Message: 
Figured out the memory allocation problem.
It turns out that the arrays are not constant per T_REC for any given HARPNUM, e.g.:

> show_info -d key=T_REC ds=hmi.sharp_cea_720s"[3311][]" seg={Bt}
T_REC	Bt	Bt_info
2013.10.24_00:48:00_TAI	Bt.fits	"Gauss, fits, 547x411"
2013.10.24_01:00:00_TAI	Bt.fits	"Gauss, fits, 547x411"
2013.10.24_01:12:00_TAI	Bt.fits	"Gauss, fits, 546x411"

So I moved all of the code that created one fixed array size per HARPNUM into the loop.
Now it creates a new array per T_REC.

File Contents

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