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