1 mbobra 1.1 /*
2 * swharp_vectorB.
3 *
4 * Created by Xudong Sun on 8/22/11.
5 * Modified to -- include ALL spaceweather keywords by Monica Bobra 25 Aug 2011
6 * -- include potential field calculation
7 * -- run only on los data 5 March 2011
8 * Bz arrays
9 * Write out abs(B) as data segment and a few keywords as SW index
10 *
11 * Use:
12 * First use the bmap module to create the in= and mask= parameters.
13 *
14 * then run this module:
15 * swharp_vectorB "in=su_mbobra.test_mmap_me[][]" /
16 * "mask=su_mbobra.test_mmap_bitmap_me[][]" "out=su_mbobra.swharp_test_v2" "dzvalue=0.001"
17 */
18
19
20 #include <jsoc_main.h>
21 #include <stdio.h>
22 mbobra 1.1 #include <stdlib.h>
23 #include <math.h>
24 #include "sharp_functions.c"
25
26 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
27 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
28 #define IN_FILES 3 /* Number of input files */
29 #define PI (3.141592653589793)
30
31 /* declaring all the functions */
32
33 int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask);
34 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask);
35 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_gamma_ptr, int *mask);
36 int readFits(char *filename, float **image, int *dims);
37 int writeFits(char *filename, float *image, int *dims);
38 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask);
39 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask);
40 int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask);
41 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask);
42 int computeJz(float *by, float *bx, int *dims, float *jz, float *mean_jz_ptr, float *us_i_ptr, int *mask);
43 mbobra 1.1 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask);
44 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr, int *mask);
45 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int *mask);
46 void greenpot(float *bx, float *by, float *bz, int nnx, int nny);
47 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, float *meanpotptr, float *totpotptr, int *mask);
48 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, float *meanshear_angleptr, float *area_w_shear_gt_45ptr, float *meanshear_anglehptr, float *area_w_shear_gt_45hptr, int *mask);
49
50 char *module_name = "swharp_vectorB"; /* Module name */
51
52 ModuleArgs_t module_args[] =
53 {
54 {ARG_STRING, "in", NULL, "Input vec mag recordset."},
55 {ARG_STRING, "mask", NULL, "Input bitmap recordset."},
56 {ARG_STRING, "out", NULL, "Output series."},
57 {ARG_FLOAT, "dzvalue", NULL, "Monopole depth."},
58 {ARG_END}
59 };
60
61
62 int DoIt(void)
63 {
64 mbobra 1.1
65 int status = DRMS_SUCCESS;
66
67 char *inQuery, *outQuery; // input series query string
68 char *maskQuery; // mask series query string
69 DRMS_RecordSet_t *inRecSet, *outRecSet, *maskRecSet;
70 DRMS_Record_t *inRec, *outRec, *maskRec;
71 DRMS_Segment_t *inSegBx, *inSegBy, *inSegBz, *outSeg, *maskSeg;
72 DRMS_Array_t *inArrayBx, *inArrayBy, *inArrayBz, *outArray, *maskArray;
73 float *bx, *by, *bz, *outData, *bh, *bt, *jz, *bpx, *bpy, *bpz;
74 int *mask;
75 int dims[2], nxny, nx, ny; // dimensions; NAXIS1 = dims[0] which is the number of columns.
76 float mean_vf = 0.0;
77 float absFlux = 0.0;
78 float mean_hf = 0.0;
79 float mean_gamma = 0.0;
80 float mean_derivative_btotal = 0.0;
81 float mean_derivative_bh = 0.0;
82 float mean_derivative_bz = 0.0;
83 float mean_jz = 0.0;
84 float us_i = 0.0;
85 mbobra 1.1 float mean_alpha = 0.0;
86 float mean_ih = 0.0;
87 float total_us_ih = 0.0;
88 float total_abs_ih = 0.0;
89 float totaljz = 0.0;
90 float totpot =0.0;
91 float meanpot = 0.0;
92 float area_w_shear_gt_45 = 0.0;
93 float meanshear_angle = 0.0;
94 float area_w_shear_gt_45h = 0.0;
95 float meanshear_angleh = 0.0;
96 int nrecs, irec, i;
97
98 /* Input */
99
100 inQuery = (char *) params_get_str(&cmdparams, "in");
101 inRecSet = drms_open_records(drms_env, inQuery, &status);
102 if (status || inRecSet->n == 0) DIE("No input data found");
103 nrecs = inRecSet->n;
104
105 /* Mask */
106 mbobra 1.1
107 maskQuery = (char *) params_get_str(&cmdparams, "mask");
108 maskRecSet = drms_open_records(drms_env, maskQuery, &status);
109 if (status || maskRecSet->n == 0) DIE("No mask data found");
110 if (maskRecSet->n != nrecs) DIE("Mask and Input series do not have a 1:1 match");
111
112 /* Output */
113
114 outQuery = (char *) params_get_str(&cmdparams, "out");
115 outRecSet = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
116 if (status) DIE("Output recordset not created");
117
118 /* Do this for each record */
119
120 for (irec = 0; irec < nrecs; irec++)
121 {
122 /* Input record and data */
123
124 inRec = inRecSet->records[irec];
125 printf("Input Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
126
127 mbobra 1.1 maskRec = maskRecSet->records[irec];
128 printf("Mask Record #%d of #%d\n", irec+1, nrecs); fflush(stdout);
129
130 inSegBx = drms_segment_lookupnum(inRec, 0); /* Assume this is Bx equivalent */
131 inSegBy = drms_segment_lookupnum(inRec, 1); /* Assume this is By equivalent */
132 inSegBz = drms_segment_lookupnum(inRec, 2); /* Assume this is Bz equivalent */
133
134 maskSeg = drms_segment_lookupnum(maskRec, 0); /* This is the bitmap */
135
136 inArrayBx = drms_segment_read(inSegBx, DRMS_TYPE_FLOAT, &status);
137 if (status) DIE("No Bx data file found. \n");
138 inArrayBy = drms_segment_read(inSegBy, DRMS_TYPE_FLOAT, &status);
139 if (status) DIE("No By data file found. \n");
140 inArrayBz = drms_segment_read(inSegBz, DRMS_TYPE_FLOAT, &status);
141 if (status) DIE("No Bz data file found. \n");
142
143 maskArray = drms_segment_read(maskSeg, DRMS_TYPE_INT, &status);
144 if (status) DIE("No mask data file found. \n");
145
146 bx = (float *)inArrayBx->data;
147 by = (float *)inArrayBy->data;
148 mbobra 1.1 bz = (float *)inArrayBz->data;
149 mask = (int *)maskArray->data;
150
151 nx = dims[0] = inArrayBz->axis[0];
152 ny = dims[1] = inArrayBz->axis[1];
153 nxny = dims[0] * dims[1];
154 if (maskArray->axis[0] != nx || maskArray->axis[1] != ny) DIE("Mask and Input series are not of the same size");
155
156 /* This is to modify the data for each PROJECTION method */
157 int flag;
158 char* value1;
159
160 value1 = drms_getkey_string(inRec, "PROJECTION", &status);
161 flag = strcmp("LambertCylindrical",value1);
162 if (flag == 0)
163 {
164 int i, j;
165 for (j = 0; j < ny; j++)
166 {
167 for (i = 0; i < nx; i++)
168 {
169 mbobra 1.1 by[j * nx + i] = - by[j * nx + i];
170 }
171 }
172 }
173
174 /* Output data */
175
176 outRec = outRecSet->records[irec];
177 drms_setlink_static(outRec, "SRCLINK", inRec->recnum);
178
179 /*===========================================*/
180 /* Malloc some arrays */
181
182 bh = (float *)malloc(nx*ny*sizeof(float));
183 bt = (float *)malloc(nx*ny*sizeof(float));
184 jz = (float *)malloc(nx*ny*sizeof(float));
185 bpx = (float *)malloc(nx*ny*sizeof(float));
186 bpy = (float *)malloc(nx*ny*sizeof(float));
187 bpz = (float *)malloc(nx*ny*sizeof(float));
188
189 /*===========================================*/
190 mbobra 1.1 /* SW Keyword computation */
191
192 if (computeAbsFlux(bz, dims, &absFlux, &mean_vf, mask))
193 {
194 absFlux = 0.0 / 0.0; // If fail, fill in NaN
195 mean_vf = 0.0 / 0.0;
196 }
197 drms_setkey_float(outRec, "USFLUX", mean_vf);
198
199 for (i=0 ;i<nxny; i++){bpz[i]=bz[i];}
200 greenpot(bpx, bpy, bpz, nx, ny);
201
202 computeBh(bx, by, bz, bh, dims, &mean_hf, mask);
203
204 if (computeGamma(bx, by, bz, bh, dims, &mean_gamma, mask)) mean_gamma = 0.0 / 0.0;
205 drms_setkey_float(outRec, "MEANGAM", mean_gamma);
206
207 computeB_total(bx, by, bz, bt, dims, mask);
208
209 if (computeBtotalderivative(bt, dims, &mean_derivative_btotal, mask)) mean_derivative_btotal = 0.0 / 0.0;
210 drms_setkey_float(outRec, "MEANGBT", mean_derivative_btotal);
211 mbobra 1.1
212 if (computeBhderivative(bh, dims, &mean_derivative_bh, mask)) mean_derivative_bh = 0.0 / 0.0;
213 drms_setkey_float(outRec, "MEANGBH", mean_derivative_bh);
214
215 if (computeBzderivative(bz, dims, &mean_derivative_bz, mask)) mean_derivative_bz = 0.0 / 0.0; // If fail, fill in NaN
216 drms_setkey_float(outRec, "MEANGBZ", mean_derivative_bz);
217
218 if(computeJz(bx, by, dims, jz, &mean_jz, &us_i, mask))
219 {
220 mean_jz = 0.0 / 0.0;
221 us_i = 0.0 /0.0;
222 }
223 drms_setkey_float(outRec, "MEANJZD", mean_jz);
224 drms_setkey_float(outRec, "TOTUSJZ", us_i);
225
226 if (computeAlpha(bz, dims, jz, &mean_alpha, mask)) mean_alpha = 0.0 / 0.0;
227 drms_setkey_float(outRec, "MEANALP", mean_alpha);
228
229 if (computeHelicity(bz, dims, jz, &mean_ih, &total_us_ih, &total_abs_ih, mask))
230 {
231 mean_ih = 0.0/0.0;
232 mbobra 1.1 total_us_ih = 0.0/0.0;
233 total_abs_ih= 0.0/0.0;
234 }
235 drms_setkey_float(outRec, "MEANJZH", mean_ih);
236 drms_setkey_float(outRec, "TOTUSJH", total_us_ih);
237 drms_setkey_float(outRec, "ABSNJZH", total_abs_ih);
238
239 if (computeSumAbsPerPolarity(bz, jz, dims, &totaljz, mask)) totaljz = 0.0 / 0.0;
240 drms_setkey_float(outRec, "SAVNCPP", totaljz);
241
242 if (computeFreeEnergy(bx, by, bpx, bpy, dims, &meanpot, &totpot, mask))
243 {
244 meanpot = 0.0 / 0.0; // If fail, fill in NaN
245 totpot = 0.0 / 0.0;
246 }
247 drms_setkey_float(outRec, "MEANPOT", meanpot);
248 drms_setkey_float(outRec, "TOTPOT", totpot);
249
250 if (computeShearAngle(bx, by, bz, bpx, bpy, bpz, dims, &meanshear_angle, &area_w_shear_gt_45, &meanshear_angleh, &area_w_shear_gt_45h, mask))
251 {
252 meanshear_angle = 0.0 / 0.0; // If fail, fill in NaN
253 mbobra 1.1 area_w_shear_gt_45 = 0.0/0.0;
254 meanshear_angleh = 0.0 / 0.0; // If fail, fill in NaN
255 area_w_shear_gt_45h = 0.0/0.0;
256 }
257 printf("meanshear_angle=%f, area_w_shear_gt_45=%f, meanshear_angleh=%f, area_w_shear_gt_45h=%f\n",meanshear_angle,area_w_shear_gt_45,meanshear_angleh,area_w_shear_gt_45h);
258 drms_setkey_float(outRec, "MEANSHR", meanshear_angle);
259 drms_setkey_float(outRec, "SHRGT45", area_w_shear_gt_45);
260 //drms_setkey_float(outRec, "MEANSHRH", meanshear_angleh);
261 //drms_setkey_float(outRec, "SHRGT45H", area_w_shear_gt_45h);
262
263
264 /*===========================================*/
265 /* Set non-SW keywords */
266
267 drms_copykey(outRec, inRec, "T_REC");
268 drms_copykey(outRec, inRec, "HARPNUM");
269
270 /* Clean up */
271 drms_free_array(inArrayBz);
272 drms_free_array(maskArray);
273
274 mbobra 1.1 }
275
276 drms_close_records(inRecSet, DRMS_FREE_RECORD);
277 drms_close_records(outRecSet, DRMS_INSERT_RECORD);
278
279 return 0;
280
281 }
|
Return to
swharp_vectorB.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps