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