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