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