1 xudong 1.1 /*===========================================
2
|
3 mbobra 1.5 The following 14 functions calculate the following spaceweather indices:
|
4 xudong 1.1
5 USFLUX Total unsigned flux in Maxwells
6 MEANGAM Mean inclination angle, gamma, in degrees
7 MEANGBT Mean value of the total field gradient, in Gauss/Mm
8 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
9 MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
10 MEANJZD Mean vertical current density, in mA/m2
11 TOTUSJZ Total unsigned vertical current, in Amperes
12 MEANALP Mean twist parameter, alpha, in 1/Mm
13 MEANJZH Mean current helicity in G2/m
14 TOTUSJH Total unsigned current helicity in G2/m
15 ABSNJZH Absolute value of the net current helicity in G2/m
16 SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
17 MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
18 TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
19 MEANSHR Mean shear angle (measured using Btotal) in degrees
20
|
21 mbobra 1.3 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
22 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
23 coordinate bitmaps are interpolated.
24
25 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
26 and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:
27
28 For conf_disambig:
29 50 : not all solutions agree (weak field method applied)
30 60 : not all solutions agree (weak field + annealed)
31 90 : all solutions agree (strong field + annealed)
32 0 : not disambiguated
33
34 For bitmap:
35 1 : weak field outside smooth bounding curve
36 2 : strong field outside smooth bounding curve
37 33 : weak field inside smooth bounding curve
38 34 : strong field inside smooth bounding curve
|
39 xudong 1.1
40 Written by Monica Bobra 15 August 2012
41 Potential Field code (appended) written by Keiji Hayashi
42
43 ===========================================*/
44 #include <math.h>
|
45 mbobra 1.9 #include <mkl.h>
|
46 xudong 1.1
|
47 mbobra 1.9 #define PI (M_PI)
|
48 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
49
50 /*===========================================*/
51
52 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
53
54 // To compute the unsigned flux, we simply calculate
55 // flux = surface integral [(vector Bz) dot (normal vector)],
56 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
57 // However, since the field is radial, we will assume cos theta = 1.
58 // Therefore the pixels only need to be corrected for the projection.
59
60 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
61 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
62 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
63 // =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
64 // =(1.30501e15)Gauss*cm^2
65
66 // The disambig mask value selects only the pixels with values of 5 or 7 -- that is,
67 // 5: pixels for which the radial acute disambiguation solution was chosen
68 // 7: pixels for which the radial acute and NRWA disambiguation agree
69 xudong 1.1
|
70 mbobra 1.9 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
71 float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
72 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
73 xudong 1.1
74 {
75
|
76 mbobra 1.14 int nx = dims[0];
77 int ny = dims[1];
|
78 mbobra 1.15 int i = 0;
79 int j = 0;
80 int count_mask = 0;
81 double sum = 0.0;
82 double err = 0.0;
|
83 mbobra 1.14 *absFlux = 0.0;
84 *mean_vf_ptr = 0.0;
85
86
|
87 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
88
|
89 mbobra 1.5 for (i = 0; i < nx; i++)
|
90 xudong 1.1 {
|
91 mbobra 1.5 for (j = 0; j < ny; j++)
|
92 xudong 1.1 {
|
93 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
94 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
95 xudong 1.1 sum += (fabs(bz[j * nx + i]));
|
96 mbobra 1.14 //printf("i,j,bz[j * nx + i]=%d,%d,%f\n",i,j,bz[j * nx + i]);
|
97 mbobra 1.9 err += bz_err[j * nx + i]*bz_err[j * nx + i];
|
98 xudong 1.1 count_mask++;
99 }
100 }
101
|
102 mbobra 1.9 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
103 *mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux
104 *count_mask_ptr = count_mask;
|
105 mbobra 1.16 //printf("cdelt1=%f\n",cdelt1);
106 //printf("rsun_obs=%f\n",rsun_obs);
107 //printf("rsun_ref=%f\n",rsun_ref);
108 //printf("CMASK=%g\n",*count_mask_ptr);
109 //printf("USFLUX=%g\n",*mean_vf_ptr);
110 //printf("sum=%f\n",sum);
111 //printf("USFLUX_err=%g\n",*mean_vf_err_ptr);
|
112 xudong 1.1 return 0;
113 }
114
115 /*===========================================*/
|
116 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
|
117 xudong 1.1 // Native units of Bh are Gauss
118
|
119 mbobra 1.9 int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
|
120 mbobra 1.3 float *mean_hf_ptr, int *mask, int *bitmask)
|
121 xudong 1.1
122 {
123
|
124 mbobra 1.14 int nx = dims[0];
125 int ny = dims[1];
|
126 mbobra 1.15 int i = 0;
127 int j = 0;
128 int count_mask = 0;
129 double sum = 0.0;
|
130 mbobra 1.9 *mean_hf_ptr = 0.0;
|
131 xudong 1.1
132 if (nx <= 0 || ny <= 0) return 1;
133
|
134 mbobra 1.5 for (i = 0; i < nx; i++)
|
135 xudong 1.1 {
|
136 mbobra 1.5 for (j = 0; j < ny; j++)
|
137 xudong 1.1 {
|
138 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
139 if isnan(by[j * nx + i]) continue;
|
140 xudong 1.1 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
141 sum += bh[j * nx + i];
|
142 mbobra 1.9 bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];
|
143 xudong 1.1 count_mask++;
144 }
145 }
146
147 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
|
148 mbobra 1.9
|
149 xudong 1.1 return 0;
150 }
151
152 /*===========================================*/
153 /* Example function 3: Calculate Gamma in units of degrees */
154 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
|
155 mbobra 1.9 // Redo calculation in radians for error analysis (since derivatives are only true in units of radians).
|
156 xudong 1.1
|
157 mbobra 1.9 int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
158 float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
|
159 xudong 1.1 {
|
160 mbobra 1.14 int nx = dims[0];
161 int ny = dims[1];
|
162 mbobra 1.15 int i = 0;
163 int j = 0;
164 int count_mask = 0;
165 double sum = 0.0;
166 double err = 0.0;
167 *mean_gamma_ptr = 0.0;
|
168 xudong 1.1
169 if (nx <= 0 || ny <= 0) return 1;
170
171 for (i = 0; i < nx; i++)
172 {
173 for (j = 0; j < ny; j++)
174 {
175 if (bh[j * nx + i] > 100)
176 {
|
177 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
178 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
179 mbobra 1.9 if isnan(bz_err[j * nx + i]) continue;
180 if isnan(bh_err[j * nx + i]) continue;
181 if (bz[j * nx + i] == 0) continue;
182 sum += (atan(fabs(bz[j * nx + i]/bh[j * nx + i] )))*(180./PI);
183 err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i]))) * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);
184 count_mask++;
|
185 xudong 1.1 }
186 }
187 }
188
189 *mean_gamma_ptr = sum/count_mask;
|
190 mbobra 1.14 *mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.0); // error in the quantity (sum)/(count_mask)
|
191 mbobra 1.16 //printf("MEANGAM=%f\n",*mean_gamma_ptr);
192 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
|
193 xudong 1.1 return 0;
194 }
195
196 /*===========================================*/
197 /* Example function 4: Calculate B_Total*/
198 // Native units of B_Total are in gauss
199
|
200 mbobra 1.9 int computeB_total(float *bx_err, float *by_err, float *bz_err, float *bt_err, float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask)
|
201 xudong 1.1 {
202
|
203 mbobra 1.14 int nx = dims[0];
204 int ny = dims[1];
|
205 mbobra 1.15 int i = 0;
206 int j = 0;
207 int count_mask = 0;
|
208 xudong 1.1
209 if (nx <= 0 || ny <= 0) return 1;
210
211 for (i = 0; i < nx; i++)
212 {
213 for (j = 0; j < ny; j++)
214 {
|
215 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
216 if isnan(by[j * nx + i]) continue;
217 if isnan(bz[j * nx + i]) continue;
|
218 xudong 1.1 bt[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]);
|
219 mbobra 1.9 bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];
|
220 xudong 1.1 }
221 }
222 return 0;
223 }
224
225 /*===========================================*/
226 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
227
|
228 mbobra 1.9 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr)
|
229 xudong 1.1 {
230
|
231 mbobra 1.14 int nx = dims[0];
232 int ny = dims[1];
|
233 mbobra 1.15 int i = 0;
234 int j = 0;
235 int count_mask = 0;
236 double sum = 0.0;
237 double err = 0.0;
|
238 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
|
239 xudong 1.1
240 if (nx <= 0 || ny <= 0) return 1;
241
242 /* brute force method of calculating the derivative (no consideration for edges) */
|
243 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
244 xudong 1.1 {
245 for (j = 0; j <= ny-1; j++)
246 {
|
247 mbobra 1.10 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
|
248 xudong 1.1 }
249 }
250
251 /* brute force method of calculating the derivative (no consideration for edges) */
252 for (i = 0; i <= nx-1; i++)
253 {
|
254 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
255 xudong 1.1 {
|
256 mbobra 1.10 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
|
257 xudong 1.1 }
258 }
259
260
|
261 mbobra 1.10 /* consider the edges */
|
262 xudong 1.1 i=0;
263 for (j = 0; j <= ny-1; j++)
264 {
|
265 mbobra 1.10 derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
|
266 xudong 1.1 }
267
268 i=nx-1;
269 for (j = 0; j <= ny-1; j++)
270 {
|
271 mbobra 1.10 derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
|
272 mbobra 1.9 }
273
|
274 xudong 1.1 j=0;
275 for (i = 0; i <= nx-1; i++)
276 {
|
277 mbobra 1.10 dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
|
278 xudong 1.1 }
279
280 j=ny-1;
281 for (i = 0; i <= nx-1; i++)
282 {
|
283 mbobra 1.10 dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
|
284 xudong 1.1 }
285
286
287 for (i = 0; i <= nx-1; i++)
288 {
289 for (j = 0; j <= ny-1; j++)
290 {
|
291 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
292 mbobra 1.5 if isnan(derx_bt[j * nx + i]) continue;
293 if isnan(dery_bt[j * nx + i]) continue;
|
294 xudong 1.1 sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ); /* Units of Gauss */
|
295 mbobra 1.9 err += (2.0)*bt_err[j * nx + i]*bt_err[j * nx + i];
|
296 xudong 1.1 count_mask++;
297 }
298 }
299
|
300 mbobra 1.9 *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
301 *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
|
302 mbobra 1.16 //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
303 //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
|
304 xudong 1.1 return 0;
305 }
306
307
308 /*===========================================*/
309 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
310
|
311 mbobra 1.9 int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)
|
312 xudong 1.1 {
313
|
314 mbobra 1.14 int nx = dims[0];
315 int ny = dims[1];
|
316 mbobra 1.15 int i = 0;
317 int j = 0;
318 int count_mask = 0;
319 double sum= 0.0;
320 double err =0.0;
|
321 mbobra 1.14 *mean_derivative_bh_ptr = 0.0;
|
322 xudong 1.1
323 if (nx <= 0 || ny <= 0) return 1;
324
325 /* brute force method of calculating the derivative (no consideration for edges) */
|
326 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
327 xudong 1.1 {
328 for (j = 0; j <= ny-1; j++)
329 {
|
330 mbobra 1.10 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
|
331 xudong 1.1 }
332 }
333
334 /* brute force method of calculating the derivative (no consideration for edges) */
335 for (i = 0; i <= nx-1; i++)
336 {
|
337 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
338 xudong 1.1 {
|
339 mbobra 1.10 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
|
340 xudong 1.1 }
341 }
342
343
|
344 mbobra 1.10 /* consider the edges */
|
345 xudong 1.1 i=0;
346 for (j = 0; j <= ny-1; j++)
347 {
|
348 mbobra 1.10 derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
|
349 xudong 1.1 }
350
351 i=nx-1;
352 for (j = 0; j <= ny-1; j++)
353 {
|
354 mbobra 1.10 derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
|
355 mbobra 1.9 }
356
|
357 xudong 1.1 j=0;
358 for (i = 0; i <= nx-1; i++)
359 {
|
360 mbobra 1.10 dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
|
361 xudong 1.1 }
362
363 j=ny-1;
364 for (i = 0; i <= nx-1; i++)
365 {
|
366 mbobra 1.10 dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
|
367 xudong 1.1 }
368
369
370 for (i = 0; i <= nx-1; i++)
371 {
372 for (j = 0; j <= ny-1; j++)
373 {
|
374 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
375 mbobra 1.5 if isnan(derx_bh[j * nx + i]) continue;
376 if isnan(dery_bh[j * nx + i]) continue;
|
377 xudong 1.1 sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ); /* Units of Gauss */
|
378 mbobra 1.9 err += (2.0)*bh_err[j * nx + i]*bh_err[j * nx + i];
|
379 xudong 1.1 count_mask++;
380 }
381 }
382
|
383 mbobra 1.9 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
384 *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
|
385 mbobra 1.16 //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
386 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
|
387 mbobra 1.9
|
388 xudong 1.1 return 0;
389 }
390
391 /*===========================================*/
392 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
393
|
394 mbobra 1.9 int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)
|
395 xudong 1.1 {
396
|
397 mbobra 1.14 int nx = dims[0];
398 int ny = dims[1];
|
399 mbobra 1.15 int i = 0;
400 int j = 0;
401 int count_mask = 0;
402 double sum = 0.0;
403 double err = 0.0;
|
404 mbobra 1.14 *mean_derivative_bz_ptr = 0.0;
|
405 xudong 1.1
406 if (nx <= 0 || ny <= 0) return 1;
407
408 /* brute force method of calculating the derivative (no consideration for edges) */
|
409 mbobra 1.10 for (i = 1; i <= nx-2; i++)
|
410 xudong 1.1 {
411 for (j = 0; j <= ny-1; j++)
412 {
|
413 mbobra 1.10 if isnan(bz[j * nx + i]) continue;
414 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
|
415 xudong 1.1 }
416 }
417
418 /* brute force method of calculating the derivative (no consideration for edges) */
419 for (i = 0; i <= nx-1; i++)
420 {
|
421 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
422 xudong 1.1 {
|
423 mbobra 1.10 if isnan(bz[j * nx + i]) continue;
424 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
|
425 xudong 1.1 }
426 }
427
428
|
429 mbobra 1.10 /* consider the edges */
|
430 xudong 1.1 i=0;
431 for (j = 0; j <= ny-1; j++)
432 {
|
433 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
434 mbobra 1.10 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
|
435 xudong 1.1 }
436
437 i=nx-1;
438 for (j = 0; j <= ny-1; j++)
439 {
|
440 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
441 mbobra 1.10 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
|
442 xudong 1.1 }
443
444 j=0;
445 for (i = 0; i <= nx-1; i++)
446 {
|
447 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
448 mbobra 1.10 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
|
449 xudong 1.1 }
450
451 j=ny-1;
452 for (i = 0; i <= nx-1; i++)
453 {
|
454 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
455 mbobra 1.10 dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
|
456 xudong 1.1 }
457
458
459 for (i = 0; i <= nx-1; i++)
460 {
461 for (j = 0; j <= ny-1; j++)
462 {
463 // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
|
464 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
465 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
466 mbobra 1.9 //if isnan(bz_err[j * nx + i]) continue;
|
467 mbobra 1.4 if isnan(derx_bz[j * nx + i]) continue;
468 if isnan(dery_bz[j * nx + i]) continue;
|
469 xudong 1.1 sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] ); /* Units of Gauss */
|
470 mbobra 1.9 err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];
|
471 xudong 1.1 count_mask++;
472 }
473 }
474
475 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
|
476 mbobra 1.9 *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
|
477 mbobra 1.16 //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
478 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
|
479 mbobra 1.9
|
480 xudong 1.1 return 0;
481 }
482
483 /*===========================================*/
484 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
485
486 // In discretized space like data pixels,
487 // the current (or curl of B) is calculated as
488 // the integration of the field Bx and By along
489 // the circumference of the data pixel divided by the area of the pixel.
490 // One form of differencing (a word for the differential operator
|
491 mbobra 1.10 // in the discretized space) of the curl is expressed as
|
492 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
493 // +dy * (By(i+1,j)+By(i,j)) / 2
494 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
495 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
496 //
497 //
498 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
499 // one must perform the following unit conversions:
|
500 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
501 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
502 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
|
503 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
|
504 mbobra 1.4 //
|
505 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
506 // In that case, we would have the following:
507 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
508 // jz * (35.0)
509 //
510 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
|
511 mbobra 1.8 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)
512 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
|
513 xudong 1.1
514
|
515 mbobra 1.9 // Comment out random number generator, which can only run on solar3
516 //int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
517 // int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
518 // float *noiseby, float *noisebz)
519
520 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
521 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
522
|
523 xudong 1.1
|
524 mbobra 1.10 {
|
525 mbobra 1.14 int nx = dims[0];
526 int ny = dims[1];
|
527 mbobra 1.15 int i = 0;
528 int j = 0;
529 int count_mask = 0;
|
530 mbobra 1.10
531 if (nx <= 0 || ny <= 0) return 1;
|
532 xudong 1.1
|
533 mbobra 1.9 /* Calculate the derivative*/
|
534 xudong 1.1 /* brute force method of calculating the derivative (no consideration for edges) */
|
535 mbobra 1.10
536 for (i = 1; i <= nx-2; i++)
|
537 xudong 1.1 {
538 for (j = 0; j <= ny-1; j++)
539 {
|
540 mbobra 1.12 if isnan(by[j * nx + i]) continue;
|
541 mbobra 1.10 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
|
542 xudong 1.1 }
543 }
544
545 for (i = 0; i <= nx-1; i++)
546 {
|
547 mbobra 1.10 for (j = 1; j <= ny-2; j++)
|
548 xudong 1.1 {
|
549 mbobra 1.12 if isnan(bx[j * nx + i]) continue;
|
550 mbobra 1.10 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
|
551 xudong 1.1 }
552 }
553
|
554 mbobra 1.10 // consider the edges
|
555 xudong 1.1 i=0;
556 for (j = 0; j <= ny-1; j++)
557 {
|
558 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
559 mbobra 1.10 derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
|
560 xudong 1.1 }
561
562 i=nx-1;
563 for (j = 0; j <= ny-1; j++)
564 {
|
565 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
566 mbobra 1.10 derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
567 }
|
568 mbobra 1.9
|
569 xudong 1.1 j=0;
570 for (i = 0; i <= nx-1; i++)
571 {
|
572 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
573 mbobra 1.10 dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
|
574 xudong 1.1 }
575
576 j=ny-1;
|
577 mbobra 1.11 for (i = 0; i <= nx-1; i++)
|
578 xudong 1.1 {
|
579 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
580 mbobra 1.10 dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
|
581 mbobra 1.9 }
582
|
583 xudong 1.1
584 for (i = 0; i <= nx-1; i++)
585 {
586 for (j = 0; j <= ny-1; j++)
587 {
|
588 mbobra 1.10 // calculate jz at all points
|
589 mbobra 1.15 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
590 jz_err[j * nx + i] = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +
|
591 mbobra 1.10 (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
|
592 mbobra 1.15 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
593 mbobra 1.10 count_mask++;
|
594 mbobra 1.5 }
|
595 mbobra 1.10 }
|
596 mbobra 1.5
|
597 mbobra 1.10 return 0;
598 }
|
599 mbobra 1.5
600 /*===========================================*/
601
|
602 mbobra 1.9
|
603 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
604 // Compute mean and total current on Jz array.
|
605 mbobra 1.6
|
606 mbobra 1.9 int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth,
607 float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
608 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
|
609 mbobra 1.5
610 {
611
|
612 mbobra 1.14 int nx = dims[0];
613 int ny = dims[1];
|
614 mbobra 1.15 int i = 0;
615 int j = 0;
616 int count_mask = 0;
617 double curl = 0.0;
618 double us_i = 0.0;
619 double err = 0.0;
|
620 mbobra 1.5
621 if (nx <= 0 || ny <= 0) return 1;
622
623 /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions*/
624 for (i = 0; i <= nx-1; i++)
625 {
626 for (j = 0; j <= ny-1; j++)
627 {
|
628 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
629 mbobra 1.4 if isnan(derx[j * nx + i]) continue;
630 if isnan(dery[j * nx + i]) continue;
|
631 mbobra 1.9 if isnan(jz[j * nx + i]) continue;
632 curl += (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
633 us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */
634 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
635 xudong 1.1 count_mask++;
|
636 mbobra 1.9 }
|
637 xudong 1.1 }
638
|
639 mbobra 1.15 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
|
640 xudong 1.1 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
|
641 mbobra 1.9 *mean_jz_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // error in the quantity MEANJZD
642
|
643 mbobra 1.4 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
|
644 mbobra 1.9 *us_i_err_ptr = (sqrt(err))*fabs((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
645
|
646 mbobra 1.16 //printf("MEANJZD=%f\n",*mean_jz_ptr);
647 //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
|
648 mbobra 1.9
|
649 mbobra 1.16 //printf("TOTUSJZ=%g\n",*us_i_ptr);
650 //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
|
651 mbobra 1.9
|
652 xudong 1.1 return 0;
653
654 }
655
|
656 mbobra 1.5 /*===========================================*/
657
658 /* Example function 10: Twist Parameter, alpha */
|
659 xudong 1.1
|
660 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
661 // for alpha is calculated in the following way (different from Leka and Barnes' approach):
662
663 // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
664 // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
665 // avg alpha = avg Jz / avg Bz
|
666 xudong 1.1
|
667 mbobra 1.6 // The sign is assigned as follows:
668 // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
669 // If this value is > 0, then alpha is > 0.
670 // If this value is < 0, then alpha is <0.
671 //
672 // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
673 // If this value is > 0, then alpha is < 0.
674 // If this value is < 0, then alpha is > 0.
675
|
676 mbobra 1.5 // The units of alpha are in 1/Mm
|
677 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
678 //
679 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
680 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
681 // = 1/Mm
682
|
683 mbobra 1.9 int computeAlpha(float *jz_err, float *bz_err, float *bz, int *dims, float *jz, float *jz_smooth, float *mean_alpha_ptr, float *mean_alpha_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
684 mbobra 1.5
|
685 xudong 1.1 {
|
686 mbobra 1.14 int nx = dims[0];
687 int ny = dims[1];
|
688 mbobra 1.15 int i = 0;
689 int j = 0;
690 int count_mask = 0;
691 double a = 0.0;
692 double b = 0.0;
693 double c = 0.0;
694 double d = 0.0;
695 double sum1 = 0.0;
696 double sum2 = 0.0;
697 double sum3 = 0.0;
698 double sum4 = 0.0;
699 double sum = 0.0;
700 double sum5 = 0.0;
701 double sum6 = 0.0;
702 double sum_err = 0.0;
|
703 xudong 1.1
704 if (nx <= 0 || ny <= 0) return 1;
705
706 for (i = 1; i < nx-1; i++)
707 {
708 for (j = 1; j < ny-1; j++)
709 {
|
710 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
711 mbobra 1.9 if isnan(jz[j * nx + i]) continue;
|
712 xudong 1.1 if isnan(bz[j * nx + i]) continue;
|
713 mbobra 1.9 if (jz[j * nx + i] == 0.0) continue;
714 if (bz_err[j * nx + i] == 0.0) continue;
715 if (bz[j * nx + i] == 0.0) continue;
716 if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i] ); a++;
717 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;
718 if (bz[j * nx + i] > 0) sum3 += ( jz[j * nx + i] ); c++;
719 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;
720 sum5 += bz[j * nx + i];
721 /* sum_err is a fractional uncertainty */
722 sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs( ( (jz[j * nx + i]) / (bz[j * nx + i]) ) *(1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
723 count_mask++;
|
724 xudong 1.1 }
725 }
|
726 mbobra 1.5
|
727 mbobra 1.9 sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */
728
|
729 mbobra 1.5 /* Determine the sign of alpha */
730 if ((sum5 > 0) && (sum3 > 0)) sum=sum;
731 if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
732 if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
733 if ((sum5 < 0) && (sum4 > 0)) sum=-sum;
734
735 *mean_alpha_ptr = sum; /* Units are 1/Mm */
|
736 mbobra 1.14 *mean_alpha_err_ptr = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.0); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent
737
|
738 mbobra 1.16 //printf("MEANALP=%f\n",*mean_alpha_ptr);
739 //printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);
|
740 mbobra 1.9
|
741 xudong 1.1 return 0;
742 }
743
744 /*===========================================*/
|
745 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
|
746 xudong 1.1
747 // The current helicity is defined as Bz*Jz and the units are G^2 / m
748 // The units of Jz are in G/pix; the units of Bz are in G.
|
749 mbobra 1.9 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
|
750 xudong 1.1 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
|
751 mbobra 1.9 // = G^2 / m.
|
752 xudong 1.1
|
753 mbobra 1.9 int computeHelicity(float *jz_err, float *jz_rms_err, float *bz_err, float *bz, int *dims, float *jz, float *mean_ih_ptr,
754 float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
755 float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
756 xudong 1.1
757 {
758
|
759 mbobra 1.14 int nx = dims[0];
760 int ny = dims[1];
|
761 mbobra 1.15 int i = 0;
762 int j = 0;
763 int count_mask = 0;
764 double sum = 0.0;
765 double sum2 = 0.0;
766 double sum_err = 0.0;
|
767 xudong 1.1
768 if (nx <= 0 || ny <= 0) return 1;
769
|
770 mbobra 1.5 for (i = 0; i < nx; i++)
|
771 xudong 1.1 {
|
772 mbobra 1.5 for (j = 0; j < ny; j++)
|
773 xudong 1.1 {
|
774 mbobra 1.9 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
775 if isnan(jz[j * nx + i]) continue;
776 if isnan(bz[j * nx + i]) continue;
777 if (bz[j * nx + i] == 0.0) continue;
778 if (jz[j * nx + i] == 0.0) continue;
779 sum += (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
|
780 mbobra 1.14 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
|
781 mbobra 1.9 sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));
782 count_mask++;
|
783 xudong 1.1 }
784 }
785
|
786 mbobra 1.9 *mean_ih_ptr = sum/count_mask ; /* Units are G^2 / m ; keyword is MEANJZH */
787 *total_us_ih_ptr = sum2 ; /* Units are G^2 / m ; keyword is TOTUSJH */
788 *total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
789
|
790 mbobra 1.14 *mean_ih_err_ptr = (sqrt(sum_err*sum_err)) / (count_mask*100.0) ; // error in the quantity MEANJZH
791 *total_us_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0) ; // error in the quantity TOTUSJH
792 *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0) ; // error in the quantity ABSNJZH
|
793 mbobra 1.9
|
794 mbobra 1.16 //printf("MEANJZH=%f\n",*mean_ih_ptr);
795 //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
|
796 mbobra 1.9
|
797 mbobra 1.16 //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
798 //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
|
799 mbobra 1.9
|
800 mbobra 1.16 //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
801 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
|
802 xudong 1.1
803 return 0;
804 }
805
806 /*===========================================*/
|
807 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
|
808 xudong 1.1
809 // The Sum of the Absolute Value per polarity is defined as the following:
810 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
811 // The units of jz are in G/pix. In this case, we would have the following:
812 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
813 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
|
814 mbobra 1.9 //
815 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
|
816 xudong 1.1
|
817 mbobra 1.9 int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
|
818 mbobra 1.3 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
819 xudong 1.1
820 {
|
821 mbobra 1.14 int nx = dims[0];
822 int ny = dims[1];
823 int i=0;
824 int j=0;
825 int count_mask=0;
|
826 mbobra 1.15 double sum1=0.0;
827 double sum2=0.0;
828 double err=0.0;
|
829 mbobra 1.14 *totaljzptr=0.0;
|
830 xudong 1.1
831 if (nx <= 0 || ny <= 0) return 1;
832
833 for (i = 0; i < nx; i++)
834 {
835 for (j = 0; j < ny; j++)
836 {
|
837 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
838 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
839 mbobra 1.9 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
840 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
841 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
842 count_mask++;
|
843 xudong 1.1 }
844 }
845
|
846 mbobra 1.9 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
847 *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
|
848 mbobra 1.16 //printf("SAVNCPP=%g\n",*totaljzptr);
849 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
|
850 mbobra 1.9
|
851 xudong 1.1 return 0;
852 }
853
854 /*===========================================*/
|
855 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
856 xudong 1.1 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
|
857 mbobra 1.11 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
858 // the integral is over B^2 dV and the 8*PI is divided at the end.
|
859 xudong 1.1 //
860 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
861 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
|
862 mbobra 1.9 // erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
863 // = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
864 // = erg/cm^3*(1.30501e15)
|
865 xudong 1.1 // = erg/cm(1/pix^2)
866
|
867 mbobra 1.9 int computeFreeEnergy(float *bx_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
868 float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
|
869 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
870
871 {
|
872 mbobra 1.14 int nx = dims[0];
873 int ny = dims[1];
|
874 mbobra 1.15 int i = 0;
875 int j = 0;
876 int count_mask = 0;
877 double sum = 0.0;
878 double sum1 = 0.0;
879 double err = 0.0;
880 *totpotptr = 0.0;
881 *meanpotptr = 0.0;
|
882 mbobra 1.14
883 if (nx <= 0 || ny <= 0) return 1;
|
884 xudong 1.1
885 for (i = 0; i < nx; i++)
886 {
887 for (j = 0; j < ny; j++)
888 {
|
889 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
890 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
891 if isnan(by[j * nx + i]) continue;
|
892 mbobra 1.13 sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
893 sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );
894 err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
|
895 xudong 1.1 count_mask++;
896 }
897 }
898
|
899 mbobra 1.13 *meanpotptr = (sum1/(8.*PI)) / (count_mask); /* Units are ergs per cubic centimeter */
|
900 mbobra 1.12 *meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
|
901 mbobra 1.9
902 /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
|
903 mbobra 1.11 *totpotptr = (sum)/(8.*PI);
904 *totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
|
905 mbobra 1.9
|
906 mbobra 1.16 //printf("MEANPOT=%g\n",*meanpotptr);
907 //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
|
908 mbobra 1.9
|
909 mbobra 1.16 //printf("TOTPOT=%g\n",*totpotptr);
910 //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
|
911 mbobra 1.9
|
912 xudong 1.1 return 0;
913 }
914
915 /*===========================================*/
|
916 mbobra 1.5 /* Example function 14: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
|
917 xudong 1.1
|
918 mbobra 1.9 int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
919 float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
|
920 xudong 1.1 {
|
921 mbobra 1.14 int nx = dims[0];
922 int ny = dims[1];
|
923 mbobra 1.15 int i = 0;
924 int j = 0;
925 int count_mask = 0;
926 double dotproduct = 0.0;
927 double magnitude_potential = 0.0;
928 double magnitude_vector = 0.0;
929 double shear_angle = 0.0;
930 double err = 0.0;
931 double sum = 0.0;
932 double count = 0.0;
933 *area_w_shear_gt_45ptr = 0.0;
934 *meanshear_angleptr = 0.0;
|
935 xudong 1.1
936 if (nx <= 0 || ny <= 0) return 1;
937
938 for (i = 0; i < nx; i++)
939 {
940 for (j = 0; j < ny; j++)
941 {
|
942 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
943 xudong 1.1 if isnan(bpx[j * nx + i]) continue;
944 if isnan(bpy[j * nx + i]) continue;
945 if isnan(bpz[j * nx + i]) continue;
946 if isnan(bz[j * nx + i]) continue;
|
947 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
948 if isnan(by[j * nx + i]) continue;
|
949 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
950 dotproduct = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
|
951 mbobra 1.9 magnitude_potential = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
952 magnitude_vector = sqrt( (bx[j * nx + i]*bx[j * nx + i]) + (by[j * nx + i]*by[j * nx + i]) + (bz[j * nx + i]*bz[j * nx + i]) );
|
953 xudong 1.1 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
954 count ++;
955 sum += shear_angle ;
|
956 mbobra 1.9 err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
|
957 xudong 1.1 if (shear_angle > 45) count_mask ++;
958 }
959 }
960
961 /* For mean 3D shear angle, area with shear greater than 45*/
|
962 mbobra 1.9 *meanshear_angleptr = (sum)/(count); /* Units are degrees */
963 *meanshear_angle_err_ptr = (sqrt(err*err))/(count); // error in the quantity (sum)/(count_mask)
|
964 mbobra 1.14 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */
|
965 mbobra 1.9
|
966 mbobra 1.16 //printf("MEANSHR=%f\n",*meanshear_angleptr);
967 //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
|
968 xudong 1.1
969 return 0;
970 }
971
972
973 /*==================KEIJI'S CODE =========================*/
974
975 // #include <omp.h>
976 #include <math.h>
977
978 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
979 {
980 /* local workings */
981 int inx, iny, i, j, n;
982 /* local array */
983 float *pfpot, *rdist;
984 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
985 rdist=(float *)malloc(sizeof(float) *nnx*nny);
986 float *bztmp;
987 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
988 /* make nan */
989 xudong 1.1 // unsigned long long llnan = 0x7ff0000000000000;
990 // float NAN = (float)(llnan);
991
992 // #pragma omp parallel for private (inx)
993 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
994 // #pragma omp parallel for private (inx)
995 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
996 // #pragma omp parallel for private (inx)
997 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
998 // #pragma omp parallel for private (inx)
999 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
1000 // #pragma omp parallel for private (inx)
1001 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1002 {
1003 float val0 = bz[nnx*iny + inx];
1004 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
1005 }}
1006
1007 // dz is the monopole depth
1008 float dz = 0.001;
1009
1010 xudong 1.1 // #pragma omp parallel for private (inx)
1011 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1012 {
1013 float rdd, rdd1, rdd2;
1014 float r;
1015 rdd1 = (float)(inx);
1016 rdd2 = (float)(iny);
1017 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
1018 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
1019 }}
1020
1021 int iwindow;
1022 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
1023 float rwindow;
1024 rwindow = (float)(iwindow);
1025 rwindow = rwindow * rwindow + 0.01; // must be of square
1026
1027 rwindow = 1.0e2; // limit the window size to be 10.
1028
1029 rwindow = sqrt(rwindow);
1030 iwindow = (int)(rwindow);
1031 xudong 1.1
1032 // #pragma omp parallel for private(inx)
1033 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1034 {
1035 float val0 = bz[nnx*iny + inx];
1036 if (isnan(val0))
1037 {
1038 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
1039 }
1040 else
1041 {
1042 float sum;
1043 sum = 0.0;
1044 int j2, i2;
1045 int j2s, j2e, i2s, i2e;
1046 j2s = iny - iwindow;
1047 j2e = iny + iwindow;
1048 if (j2s < 0){j2s = 0;}
1049 if (j2e > nny){j2e = nny;}
1050 i2s = inx - iwindow;
1051 i2e = inx + iwindow;
1052 xudong 1.1 if (i2s < 0){i2s = 0;}
1053 if (i2e > nnx){i2e = nnx;}
1054
1055 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1056 {
1057 float val1 = bztmp[nnx*j2 + i2];
1058 float rr, r1, r2;
1059 // r1 = (float)(i2 - inx);
1060 // r2 = (float)(j2 - iny);
1061 // rr = r1*r1 + r2*r2;
1062 // if (rr < rwindow)
1063 // {
1064 int di, dj;
1065 di = abs(i2 - inx);
1066 dj = abs(j2 - iny);
1067 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1068 // }
1069 } }
1070 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1071 }
1072 } } // end of OpenMP parallelism
1073 xudong 1.1
1074 // #pragma omp parallel for private(inx)
1075 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1076 {
1077 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1078 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1079 } } // end of OpenMP parallelism
1080
1081 free(rdist);
1082 free(pfpot);
1083 free(bztmp);
1084 } // end of void func. greenpot
1085
1086
1087 /*===========END OF KEIJI'S CODE =========================*/
|
1088 mbobra 1.14
1089 char *sw_functions_version() // Returns CVS version of sw_functions.c
1090 {
|
1091 mbobra 1.16 return strdup("$Id: sw_functions.c,v 1.15 2013/07/08 23:02:22 mbobra Exp $");
|
1092 mbobra 1.14 }
1093
|
1094 xudong 1.1 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps