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.11 //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.11 //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
567 //printf("jz[j * nx + i]=%f,i=%d,j=%d\n,",jz[j * nx + i],i,j); // tmp.txt
568 //printf("i=%d,j=%d,jz[j * nx + i]=%f\n,",i,j,jz[j * nx + i]); // tmp1.txt
569
570
|
571 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]) +
572 (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
573 jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);
574 count_mask++;
|
575 mbobra 1.5 }
|
576 mbobra 1.10 }
|
577 mbobra 1.5
|
578 mbobra 1.10 return 0;
579 }
|
580 mbobra 1.5
581 /*===========================================*/
582
|
583 mbobra 1.9
|
584 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
585 // Compute mean and total current on Jz array.
|
586 mbobra 1.6
|
587 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,
588 float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
589 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
|
590 mbobra 1.5
591 {
592
593 int nx = dims[0], ny = dims[1];
594 int i, j, count_mask=0;
595
596 if (nx <= 0 || ny <= 0) return 1;
597
|
598 mbobra 1.9 float curl,us_i,test_perimeter,mean_curl,err=0.0;
|
599 mbobra 1.5
600
601 /* 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*/
602 for (i = 0; i <= nx-1; i++)
603 {
604 for (j = 0; j <= ny-1; j++)
605 {
|
606 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
607 mbobra 1.4 if isnan(derx[j * nx + i]) continue;
608 if isnan(dery[j * nx + i]) continue;
|
609 mbobra 1.9 if isnan(jz[j * nx + i]) continue;
610 curl += (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
611 us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */
612 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
613 xudong 1.1 count_mask++;
|
614 mbobra 1.9 }
|
615 xudong 1.1 }
616
|
617 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 */
|
618 xudong 1.1 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
|
619 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
620
|
621 mbobra 1.4 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
|
622 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
623
624 printf("MEANJZD=%f\n",*mean_jz_ptr);
625 printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
626
627 printf("TOTUSJZ=%g\n",*us_i_ptr);
628 printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
629
|
630 xudong 1.1 return 0;
631
632 }
633
|
634 mbobra 1.5 /*===========================================*/
635
636 /* Example function 10: Twist Parameter, alpha */
|
637 xudong 1.1
|
638 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
639 // for alpha is calculated in the following way (different from Leka and Barnes' approach):
640
641 // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
642 // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
643 // avg alpha = avg Jz / avg Bz
|
644 xudong 1.1
|
645 mbobra 1.6 // The sign is assigned as follows:
646 // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
647 // If this value is > 0, then alpha is > 0.
648 // If this value is < 0, then alpha is <0.
649 //
650 // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
651 // If this value is > 0, then alpha is < 0.
652 // If this value is < 0, then alpha is > 0.
653
|
654 mbobra 1.5 // The units of alpha are in 1/Mm
|
655 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
656 //
657 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
658 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
659 // = 1/Mm
660
|
661 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)
|
662 mbobra 1.5
|
663 xudong 1.1 {
664 int nx = dims[0], ny = dims[1];
|
665 mbobra 1.9 int i, j, count_mask, a,b,c,d=0;
|
666 xudong 1.1
667 if (nx <= 0 || ny <= 0) return 1;
668
|
669 mbobra 1.9 float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6, sum_err=0.0;
|
670 xudong 1.1
671 for (i = 1; i < nx-1; i++)
672 {
673 for (j = 1; j < ny-1; j++)
674 {
|
675 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
676 mbobra 1.9 //if isnan(jz_smooth[j * nx + i]) continue;
677 if isnan(jz[j * nx + i]) continue;
|
678 xudong 1.1 if isnan(bz[j * nx + i]) continue;
|
679 mbobra 1.9 //if (jz_smooth[j * nx + i] == 0) continue;
680 if (jz[j * nx + i] == 0.0) continue;
681 if (bz_err[j * nx + i] == 0.0) continue;
682 if (bz[j * nx + i] == 0.0) continue;
683 if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i] ); a++;
684 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;
685 //if (bz[j * nx + i] > 0) sum3 += ( jz_smooth[j * nx + i]);
686 //if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
687 if (bz[j * nx + i] > 0) sum3 += ( jz[j * nx + i] ); c++;
688 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;
689 sum5 += bz[j * nx + i];
690 /* sum_err is a fractional uncertainty */
691 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.));
692 count_mask++;
|
693 xudong 1.1 }
694 }
|
695 mbobra 1.5
|
696 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 */
697
|
698 mbobra 1.5 /* Determine the sign of alpha */
699 if ((sum5 > 0) && (sum3 > 0)) sum=sum;
700 if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
701 if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
702 if ((sum5 < 0) && (sum4 > 0)) sum=-sum;
703
704 *mean_alpha_ptr = sum; /* Units are 1/Mm */
|
705 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
706
707 printf("a=%d\n",a);
708 printf("b=%d\n",b);
709 printf("d=%d\n",d);
710 printf("c=%d\n",c);
711
712 printf("MEANALP=%f\n",*mean_alpha_ptr);
713 printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);
714
|
715 xudong 1.1 return 0;
716 }
717
718 /*===========================================*/
|
719 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
|
720 xudong 1.1
721 // The current helicity is defined as Bz*Jz and the units are G^2 / m
722 // The units of Jz are in G/pix; the units of Bz are in G.
|
723 mbobra 1.9 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
|
724 xudong 1.1 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
|
725 mbobra 1.9 // = G^2 / m.
|
726 xudong 1.1
|
727 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,
728 float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
729 float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
730 xudong 1.1
731 {
732
733 int nx = dims[0], ny = dims[1];
734 int i, j, count_mask=0;
735
736 if (nx <= 0 || ny <= 0) return 1;
737
|
738 mbobra 1.9 float sum,sum2,sum_err=0.0;
|
739 xudong 1.1
|
740 mbobra 1.5 for (i = 0; i < nx; i++)
|
741 xudong 1.1 {
|
742 mbobra 1.5 for (j = 0; j < ny; j++)
|
743 xudong 1.1 {
|
744 mbobra 1.9 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
745 if isnan(jz[j * nx + i]) continue;
746 if isnan(bz[j * nx + i]) continue;
747 if (bz[j * nx + i] == 0.0) continue;
748 if (jz[j * nx + i] == 0.0) continue;
749 sum += (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
750 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
751 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));
752 count_mask++;
|
753 xudong 1.1 }
754 }
755
|
756 mbobra 1.9 *mean_ih_ptr = sum/count_mask ; /* Units are G^2 / m ; keyword is MEANJZH */
757 *total_us_ih_ptr = sum2 ; /* Units are G^2 / m ; keyword is TOTUSJH */
758 *total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
759
760 *mean_ih_err_ptr = (sqrt(sum_err*sum_err)) / (count_mask*100.) ; // error in the quantity MEANJZH
761 *total_us_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.) ; // error in the quantity TOTUSJH
762 *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.) ; // error in the quantity ABSNJZH
763
764 printf("MEANJZH=%f\n",*mean_ih_ptr);
765 printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
766
767 printf("TOTUSJH=%f\n",*total_us_ih_ptr);
768 printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
769
770 printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
771 printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
|
772 xudong 1.1
773 return 0;
774 }
775
776 /*===========================================*/
|
777 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
|
778 xudong 1.1
779 // The Sum of the Absolute Value per polarity is defined as the following:
780 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
781 // The units of jz are in G/pix. In this case, we would have the following:
782 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
783 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
|
784 mbobra 1.9 //
785 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
|
786 xudong 1.1
|
787 mbobra 1.9 int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
|
788 mbobra 1.3 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
789 xudong 1.1
790 {
791 int nx = dims[0], ny = dims[1];
792 int i, j, count_mask=0;
793
794 if (nx <= 0 || ny <= 0) return 1;
795
796 *totaljzptr=0.0;
|
797 mbobra 1.9 float sum1,sum2,err=0.0;
|
798 xudong 1.1
799 for (i = 0; i < nx; i++)
800 {
801 for (j = 0; j < ny; j++)
802 {
|
803 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
804 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
805 mbobra 1.9 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
806 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
807 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
808 count_mask++;
|
809 xudong 1.1 }
810 }
811
|
812 mbobra 1.9 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
813 *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
814 printf("SAVNCPP=%g\n",*totaljzptr);
815 printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
816
|
817 xudong 1.1 return 0;
818 }
819
820 /*===========================================*/
|
821 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
822 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
|
823 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,
824 // the integral is over B^2 dV and the 8*PI is divided at the end.
|
825 xudong 1.1 //
826 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
827 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
|
828 mbobra 1.9 // erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
829 // = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
830 // = erg/cm^3*(1.30501e15)
|
831 xudong 1.1 // = erg/cm(1/pix^2)
832
|
833 mbobra 1.9 int computeFreeEnergy(float *bx_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
834 float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
|
835 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
836
837 {
838 int nx = dims[0], ny = dims[1];
839 int i, j, count_mask=0;
840
841 if (nx <= 0 || ny <= 0) return 1;
842
843 *totpotptr=0.0;
844 *meanpotptr=0.0;
|
845 mbobra 1.9 float sum,err=0.0;
|
846 xudong 1.1
847 for (i = 0; i < nx; i++)
848 {
849 for (j = 0; j < ny; j++)
850 {
|
851 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
852 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
853 if isnan(by[j * nx + i]) continue;
|
854 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);
|
855 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]);
|
856 xudong 1.1 count_mask++;
857 }
858 }
859
|
860 mbobra 1.11 *meanpotptr = (sum/(8.*PI)) / (count_mask); /* Units are ergs per cubic centimeter */
|
861 mbobra 1.10 *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
|
862 mbobra 1.9
863 /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
|
864 mbobra 1.11 *totpotptr = (sum)/(8.*PI);
865 *totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
|
866 mbobra 1.9
867 printf("MEANPOT=%g\n",*meanpotptr);
868 printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
869
870 printf("TOTPOT=%g\n",*totpotptr);
871 printf("TOTPOT_err=%g\n",*totpot_err_ptr);
872
|
873 xudong 1.1 return 0;
874 }
875
876 /*===========================================*/
|
877 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 */
|
878 xudong 1.1
|
879 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,
880 float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
|
881 xudong 1.1 {
882 int nx = dims[0], ny = dims[1];
883 int i, j;
884
885 if (nx <= 0 || ny <= 0) return 1;
886
|
887 mbobra 1.9 //*area_w_shear_gt_45ptr=0.0;
888 //*meanshear_angleptr=0.0;
889 float dotproduct, magnitude_potential, magnitude_vector, shear_angle,err=0.0, sum = 0.0, count=0.0, count_mask=0.0;
|
890 xudong 1.1
891 for (i = 0; i < nx; i++)
892 {
893 for (j = 0; j < ny; j++)
894 {
|
895 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
896 xudong 1.1 if isnan(bpx[j * nx + i]) continue;
897 if isnan(bpy[j * nx + i]) continue;
898 if isnan(bpz[j * nx + i]) continue;
899 if isnan(bz[j * nx + i]) continue;
|
900 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
901 if isnan(by[j * nx + i]) continue;
|
902 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
903 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]);
|
904 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]));
905 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]) );
|
906 xudong 1.1 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
907 count ++;
908 sum += shear_angle ;
|
909 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])));
|
910 xudong 1.1 if (shear_angle > 45) count_mask ++;
911 }
912 }
913
914 /* For mean 3D shear angle, area with shear greater than 45*/
|
915 mbobra 1.9 *meanshear_angleptr = (sum)/(count); /* Units are degrees */
916 *meanshear_angle_err_ptr = (sqrt(err*err))/(count); // error in the quantity (sum)/(count_mask)
917 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */
918
919 printf("MEANSHR=%f\n",*meanshear_angleptr);
920 printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
|
921 xudong 1.1
922 return 0;
923 }
924
925
926 /*==================KEIJI'S CODE =========================*/
927
928 // #include <omp.h>
929 #include <math.h>
930
931 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
932 {
933 /* local workings */
934 int inx, iny, i, j, n;
935 /* local array */
936 float *pfpot, *rdist;
937 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
938 rdist=(float *)malloc(sizeof(float) *nnx*nny);
939 float *bztmp;
940 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
941 /* make nan */
942 xudong 1.1 // unsigned long long llnan = 0x7ff0000000000000;
943 // float NAN = (float)(llnan);
944
945 // #pragma omp parallel for private (inx)
946 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
947 // #pragma omp parallel for private (inx)
948 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
949 // #pragma omp parallel for private (inx)
950 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
951 // #pragma omp parallel for private (inx)
952 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
953 // #pragma omp parallel for private (inx)
954 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
955 {
956 float val0 = bz[nnx*iny + inx];
957 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
958 }}
959
960 // dz is the monopole depth
961 float dz = 0.001;
962
963 xudong 1.1 // #pragma omp parallel for private (inx)
964 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
965 {
966 float rdd, rdd1, rdd2;
967 float r;
968 rdd1 = (float)(inx);
969 rdd2 = (float)(iny);
970 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
971 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
972 }}
973
974 int iwindow;
975 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
976 float rwindow;
977 rwindow = (float)(iwindow);
978 rwindow = rwindow * rwindow + 0.01; // must be of square
979
980 rwindow = 1.0e2; // limit the window size to be 10.
981
982 rwindow = sqrt(rwindow);
983 iwindow = (int)(rwindow);
984 xudong 1.1
985 // #pragma omp parallel for private(inx)
986 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
987 {
988 float val0 = bz[nnx*iny + inx];
989 if (isnan(val0))
990 {
991 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
992 }
993 else
994 {
995 float sum;
996 sum = 0.0;
997 int j2, i2;
998 int j2s, j2e, i2s, i2e;
999 j2s = iny - iwindow;
1000 j2e = iny + iwindow;
1001 if (j2s < 0){j2s = 0;}
1002 if (j2e > nny){j2e = nny;}
1003 i2s = inx - iwindow;
1004 i2e = inx + iwindow;
1005 xudong 1.1 if (i2s < 0){i2s = 0;}
1006 if (i2e > nnx){i2e = nnx;}
1007
1008 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1009 {
1010 float val1 = bztmp[nnx*j2 + i2];
1011 float rr, r1, r2;
1012 // r1 = (float)(i2 - inx);
1013 // r2 = (float)(j2 - iny);
1014 // rr = r1*r1 + r2*r2;
1015 // if (rr < rwindow)
1016 // {
1017 int di, dj;
1018 di = abs(i2 - inx);
1019 dj = abs(j2 - iny);
1020 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1021 // }
1022 } }
1023 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1024 }
1025 } } // end of OpenMP parallelism
1026 xudong 1.1
1027 // #pragma omp parallel for private(inx)
1028 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1029 {
1030 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1031 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1032 } } // end of OpenMP parallelism
1033
1034 free(rdist);
1035 free(pfpot);
1036 free(bztmp);
1037 } // end of void func. greenpot
1038
1039
1040 /*===========END OF KEIJI'S CODE =========================*/
1041 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps