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