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