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