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.18 if isnan(jz[j * nx + i]) continue;
|
850 mbobra 1.9 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
851 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
852 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
853 count_mask++;
|
854 xudong 1.1 }
855 }
856
|
857 mbobra 1.9 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
858 *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
|
859 mbobra 1.16 //printf("SAVNCPP=%g\n",*totaljzptr);
860 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
|
861 mbobra 1.9
|
862 xudong 1.1 return 0;
863 }
864
865 /*===========================================*/
|
866 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
867 xudong 1.1 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
|
868 mbobra 1.11 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
869 // the integral is over B^2 dV and the 8*PI is divided at the end.
|
870 xudong 1.1 //
871 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
872 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
|
873 mbobra 1.9 // erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
874 // = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
875 // = erg/cm^3*(1.30501e15)
|
876 xudong 1.1 // = erg/cm(1/pix^2)
877
|
878 mbobra 1.9 int computeFreeEnergy(float *bx_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
879 float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
|
880 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
881
882 {
|
883 mbobra 1.14 int nx = dims[0];
884 int ny = dims[1];
|
885 mbobra 1.15 int i = 0;
886 int j = 0;
887 int count_mask = 0;
888 double sum = 0.0;
889 double sum1 = 0.0;
890 double err = 0.0;
891 *totpotptr = 0.0;
892 *meanpotptr = 0.0;
|
893 mbobra 1.14
894 if (nx <= 0 || ny <= 0) return 1;
|
895 xudong 1.1
896 for (i = 0; i < nx; i++)
897 {
898 for (j = 0; j < ny; j++)
899 {
|
900 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
901 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
902 if isnan(by[j * nx + i]) continue;
|
903 mbobra 1.13 sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
904 sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );
905 err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
|
906 xudong 1.1 count_mask++;
907 }
908 }
909
|
910 mbobra 1.13 *meanpotptr = (sum1/(8.*PI)) / (count_mask); /* Units are ergs per cubic centimeter */
|
911 mbobra 1.12 *meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
|
912 mbobra 1.9
913 /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
|
914 mbobra 1.11 *totpotptr = (sum)/(8.*PI);
915 *totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
|
916 mbobra 1.9
|
917 mbobra 1.16 //printf("MEANPOT=%g\n",*meanpotptr);
918 //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
|
919 mbobra 1.9
|
920 mbobra 1.16 //printf("TOTPOT=%g\n",*totpotptr);
921 //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
|
922 mbobra 1.9
|
923 xudong 1.1 return 0;
924 }
925
926 /*===========================================*/
|
927 mbobra 1.5 /* Example function 14: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
|
928 xudong 1.1
|
929 mbobra 1.9 int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
930 float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
|
931 xudong 1.1 {
|
932 mbobra 1.14 int nx = dims[0];
933 int ny = dims[1];
|
934 mbobra 1.15 int i = 0;
935 int j = 0;
936 int count_mask = 0;
937 double dotproduct = 0.0;
938 double magnitude_potential = 0.0;
939 double magnitude_vector = 0.0;
940 double shear_angle = 0.0;
941 double err = 0.0;
942 double sum = 0.0;
943 double count = 0.0;
944 *area_w_shear_gt_45ptr = 0.0;
945 *meanshear_angleptr = 0.0;
|
946 xudong 1.1
947 if (nx <= 0 || ny <= 0) return 1;
948
949 for (i = 0; i < nx; i++)
950 {
951 for (j = 0; j < ny; j++)
952 {
|
953 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
954 xudong 1.1 if isnan(bpx[j * nx + i]) continue;
955 if isnan(bpy[j * nx + i]) continue;
956 if isnan(bpz[j * nx + i]) continue;
957 if isnan(bz[j * nx + i]) continue;
|
958 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
959 if isnan(by[j * nx + i]) continue;
|
960 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
961 dotproduct = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
|
962 mbobra 1.9 magnitude_potential = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
963 magnitude_vector = sqrt( (bx[j * nx + i]*bx[j * nx + i]) + (by[j * nx + i]*by[j * nx + i]) + (bz[j * nx + i]*bz[j * nx + i]) );
|
964 xudong 1.1 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
965 count ++;
966 sum += shear_angle ;
|
967 mbobra 1.9 err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
|
968 xudong 1.1 if (shear_angle > 45) count_mask ++;
969 }
970 }
971
972 /* For mean 3D shear angle, area with shear greater than 45*/
|
973 mbobra 1.9 *meanshear_angleptr = (sum)/(count); /* Units are degrees */
974 *meanshear_angle_err_ptr = (sqrt(err*err))/(count); // error in the quantity (sum)/(count_mask)
|
975 mbobra 1.14 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */
|
976 mbobra 1.9
|
977 mbobra 1.16 //printf("MEANSHR=%f\n",*meanshear_angleptr);
978 //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
|
979 xudong 1.1
980 return 0;
981 }
982
983
984 /*==================KEIJI'S CODE =========================*/
985
986 // #include <omp.h>
987 #include <math.h>
988
989 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
990 {
991 /* local workings */
992 int inx, iny, i, j, n;
993 /* local array */
994 float *pfpot, *rdist;
995 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
996 rdist=(float *)malloc(sizeof(float) *nnx*nny);
997 float *bztmp;
998 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
999 /* make nan */
1000 xudong 1.1 // unsigned long long llnan = 0x7ff0000000000000;
1001 // float NAN = (float)(llnan);
1002
1003 // #pragma omp parallel for private (inx)
1004 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
1005 // #pragma omp parallel for private (inx)
1006 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
1007 // #pragma omp parallel for private (inx)
1008 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
1009 // #pragma omp parallel for private (inx)
1010 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
1011 // #pragma omp parallel for private (inx)
1012 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1013 {
1014 float val0 = bz[nnx*iny + inx];
1015 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
1016 }}
1017
1018 // dz is the monopole depth
1019 float dz = 0.001;
1020
1021 xudong 1.1 // #pragma omp parallel for private (inx)
1022 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1023 {
1024 float rdd, rdd1, rdd2;
1025 float r;
1026 rdd1 = (float)(inx);
1027 rdd2 = (float)(iny);
1028 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
1029 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
1030 }}
1031
1032 int iwindow;
1033 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
1034 float rwindow;
1035 rwindow = (float)(iwindow);
1036 rwindow = rwindow * rwindow + 0.01; // must be of square
1037
1038 rwindow = 1.0e2; // limit the window size to be 10.
1039
1040 rwindow = sqrt(rwindow);
1041 iwindow = (int)(rwindow);
1042 xudong 1.1
1043 // #pragma omp parallel for private(inx)
1044 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1045 {
1046 float val0 = bz[nnx*iny + inx];
1047 if (isnan(val0))
1048 {
1049 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
1050 }
1051 else
1052 {
1053 float sum;
1054 sum = 0.0;
1055 int j2, i2;
1056 int j2s, j2e, i2s, i2e;
1057 j2s = iny - iwindow;
1058 j2e = iny + iwindow;
1059 if (j2s < 0){j2s = 0;}
1060 if (j2e > nny){j2e = nny;}
1061 i2s = inx - iwindow;
1062 i2e = inx + iwindow;
1063 xudong 1.1 if (i2s < 0){i2s = 0;}
1064 if (i2e > nnx){i2e = nnx;}
1065
1066 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1067 {
1068 float val1 = bztmp[nnx*j2 + i2];
1069 float rr, r1, r2;
1070 // r1 = (float)(i2 - inx);
1071 // r2 = (float)(j2 - iny);
1072 // rr = r1*r1 + r2*r2;
1073 // if (rr < rwindow)
1074 // {
1075 int di, dj;
1076 di = abs(i2 - inx);
1077 dj = abs(j2 - iny);
1078 sum = sum + val1 * rdist[nnx * dj + di] * dz;
1079 // }
1080 } }
1081 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1082 }
1083 } } // end of OpenMP parallelism
1084 xudong 1.1
1085 // #pragma omp parallel for private(inx)
1086 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1087 {
1088 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1089 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1090 } } // end of OpenMP parallelism
1091
1092 free(rdist);
1093 free(pfpot);
1094 free(bztmp);
1095 } // end of void func. greenpot
1096
1097
1098 /*===========END OF KEIJI'S CODE =========================*/
|
1099 mbobra 1.14
1100 char *sw_functions_version() // Returns CVS version of sw_functions.c
1101 {
|
1102 mbobra 1.18 return strdup("$Id: sw_functions.c,v 1.17 2013/07/24 02:35:08 mbobra Exp $");
|
1103 mbobra 1.14 }
1104
|
1105 xudong 1.1 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps