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