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