1 mbobra 1.37
|
2 xudong 1.1 /*===========================================
|
3 xudong 1.27
|
4 mbobra 1.40 The following functions calculate these spaceweather indices from the vector magnetic field data:
|
5 xudong 1.27
6 USFLUX Total unsigned flux in Maxwells
7 MEANGAM Mean inclination angle, gamma, in degrees
8 MEANGBT Mean value of the total field gradient, in Gauss/Mm
9 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
10 MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
11 MEANJZD Mean vertical current density, in mA/m2
12 TOTUSJZ Total unsigned vertical current, in Amperes
13 MEANALP Mean twist parameter, alpha, in 1/Mm
14 MEANJZH Mean current helicity in G2/m
15 TOTUSJH Total unsigned current helicity in G2/m
16 ABSNJZH Absolute value of the net current helicity in G2/m
17 SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
18 MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
19 TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
20 MEANSHR Mean shear angle (measured using Btotal) in degrees
|
21 mbobra 1.40 CMASK The total number of pixels that contributed to the calculation of all the indices listed above
22
23 And these spaceweather indices from the line-of-sight magnetic field data:
24 USFLUXL Total unsigned flux in Maxwells
25 MEANGBL Mean value of the line-of-sight field gradient, in Gauss/Mm
26 CMASKL The total number of pixels that contributed to the calculation of USFLUXL and MEANGBL
|
27 mbobra 1.29 R_VALUE Karel Schrijver's R parameter
|
28 xudong 1.27
29 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
30 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
31 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
32 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
33 contain nearest-neighbor bitmaps).
34
35 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
36 and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:
37
38 For conf_disambig:
39 50 : not all solutions agree (weak field method applied)
40 60 : not all solutions agree (weak field + annealed)
41 90 : all solutions agree (strong field + annealed)
42 0 : not disambiguated
43
44 For bitmap:
45 1 : weak field outside smooth bounding curve
46 2 : strong field outside smooth bounding curve
47 33 : weak field inside smooth bounding curve
48 34 : strong field inside smooth bounding curve
49 xudong 1.27
50 Written by Monica Bobra 15 August 2012
51 Potential Field code (appended) written by Keiji Hayashi
52 Error analysis modification 21 October 2013
53
54 ===========================================*/
|
55 xudong 1.1 #include <math.h>
|
56 mbobra 1.9 #include <mkl.h>
|
57 xudong 1.1
|
58 mbobra 1.9 #define PI (M_PI)
|
59 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
60
61 /*===========================================*/
62
63 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
64
|
65 xudong 1.27 // To compute the unsigned flux, we simply calculate
|
66 xudong 1.1 // flux = surface integral [(vector Bz) dot (normal vector)],
67 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
68 // However, since the field is radial, we will assume cos theta = 1.
69 // Therefore the pixels only need to be corrected for the projection.
70
71 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
72 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
73 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
|
74 mbobra 1.20 // =Gauss*cm^2
|
75 xudong 1.1
|
76 xudong 1.27 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
77 float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
|
78 mbobra 1.9 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
79 xudong 1.1
80 {
|
81 xudong 1.27
|
82 mbobra 1.14 int nx = dims[0];
83 int ny = dims[1];
|
84 mbobra 1.15 int i = 0;
85 int j = 0;
86 int count_mask = 0;
87 double sum = 0.0;
88 double err = 0.0;
|
89 mbobra 1.14 *absFlux = 0.0;
90 *mean_vf_ptr = 0.0;
|
91 xudong 1.27
92
|
93 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
94 xudong 1.27
95 for (i = 0; i < nx; i++)
|
96 xudong 1.1 {
|
97 mbobra 1.34 for (j = 0; j < ny; j++)
98 {
|
99 xudong 1.27 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
100 if isnan(bz[j * nx + i]) continue;
101 sum += (fabs(bz[j * nx + i]));
102 err += bz_err[j * nx + i]*bz_err[j * nx + i];
103 count_mask++;
|
104 mbobra 1.34 }
|
105 xudong 1.1 }
|
106 xudong 1.27
107 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
108 *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
109 *count_mask_ptr = count_mask;
110 return 0;
|
111 xudong 1.1 }
112
113 /*===========================================*/
|
114 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
|
115 xudong 1.1 // Native units of Bh are Gauss
116
|
117 xudong 1.27 int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
|
118 mbobra 1.3 float *mean_hf_ptr, int *mask, int *bitmask)
|
119 xudong 1.1
120 {
|
121 xudong 1.27
|
122 mbobra 1.14 int nx = dims[0];
123 int ny = dims[1];
|
124 mbobra 1.15 int i = 0;
|
125 xudong 1.27 int j = 0;
|
126 mbobra 1.15 int count_mask = 0;
|
127 xudong 1.27 double sum = 0.0;
|
128 mbobra 1.9 *mean_hf_ptr = 0.0;
|
129 xudong 1.27
|
130 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
131 xudong 1.27
132 for (i = 0; i < nx; i++)
133 {
|
134 mbobra 1.5 for (j = 0; j < ny; j++)
|
135 xudong 1.27 {
136 if isnan(bx[j * nx + i])
137 {
138 bh[j * nx + i] = NAN;
139 bh_err[j * nx + i] = NAN;
140 continue;
141 }
142 if isnan(by[j * nx + i])
143 {
144 bh[j * nx + i] = NAN;
145 bh_err[j * nx + i] = NAN;
146 continue;
147 }
148 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
149 sum += bh[j * nx + i];
150 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];
151 count_mask++;
152 }
153 }
154
|
155 xudong 1.1 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
|
156 xudong 1.27
|
157 xudong 1.1 return 0;
158 }
159
160 /*===========================================*/
161 /* Example function 3: Calculate Gamma in units of degrees */
162 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
|
163 xudong 1.27 //
|
164 mbobra 1.20 // Error analysis calculations are done in radians (since derivatives are only true in units of radians),
165 // and multiplied by (180./PI) at the end for consistency in units.
|
166 xudong 1.1
|
167 mbobra 1.9 int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
168 float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
|
169 xudong 1.1 {
|
170 mbobra 1.14 int nx = dims[0];
171 int ny = dims[1];
|
172 mbobra 1.15 int i = 0;
173 int j = 0;
174 int count_mask = 0;
175 double sum = 0.0;
176 double err = 0.0;
177 *mean_gamma_ptr = 0.0;
|
178 xudong 1.27
|
179 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
180 xudong 1.27
181 for (i = 0; i < nx; i++)
182 {
183 for (j = 0; j < ny; j++)
184 {
185 if (bh[j * nx + i] > 100)
186 {
187 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
188 if isnan(bz[j * nx + i]) continue;
189 if isnan(bz_err[j * nx + i]) continue;
190 if isnan(bh_err[j * nx + i]) continue;
191 if isnan(bh[j * nx + i]) continue;
192 if (bz[j * nx + i] == 0) continue;
193 sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
194 err += (1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))))*(1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])))) *
195 ( ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) +
196 ((bh[j * nx + i]*bh[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i])) );
197 count_mask++;
198 }
199 }
200 }
201 xudong 1.27
202 *mean_gamma_ptr = sum/count_mask;
203 *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
204 //printf("MEANGAM=%f\n",*mean_gamma_ptr);
205 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
206 return 0;
|
207 xudong 1.1 }
208
209 /*===========================================*/
210 /* Example function 4: Calculate B_Total*/
211 // Native units of B_Total are in gauss
212
|
213 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)
|
214 xudong 1.1 {
|
215 xudong 1.27
|
216 mbobra 1.14 int nx = dims[0];
217 int ny = dims[1];
|
218 mbobra 1.15 int i = 0;
219 int j = 0;
220 int count_mask = 0;
|
221 xudong 1.1
222 if (nx <= 0 || ny <= 0) return 1;
|
223 xudong 1.27
224 for (i = 0; i < nx; i++)
225 {
226 for (j = 0; j < ny; j++)
227 {
228 if isnan(bx[j * nx + i])
229 {
230 bt[j * nx + i] = NAN;
231 bt_err[j * nx + i] = NAN;
232 continue;
233 }
234 if isnan(by[j * nx + i])
235 {
236 bt[j * nx + i] = NAN;
237 bt_err[j * nx + i] = NAN;
238 continue;
239 }
240 if isnan(bz[j * nx + i])
241 {
242 bt[j * nx + i] = NAN;
243 bt_err[j * nx + i] = NAN;
244 xudong 1.27 continue;
245 }
246 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]);
247 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];
248 }
249 }
250 return 0;
|
251 xudong 1.1 }
252
253 /*===========================================*/
254 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
255
|
256 mbobra 1.37 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, float *err_termAt, float *err_termBt)
|
257 xudong 1.1 {
|
258 xudong 1.27
|
259 mbobra 1.14 int nx = dims[0];
260 int ny = dims[1];
|
261 mbobra 1.15 int i = 0;
262 int j = 0;
263 int count_mask = 0;
|
264 xudong 1.27 double sum = 0.0;
|
265 mbobra 1.15 double err = 0.0;
|
266 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
|
267 xudong 1.27
|
268 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
269 xudong 1.27
270 /* brute force method of calculating the derivative (no consideration for edges) */
271 for (i = 1; i <= nx-2; i++)
272 {
|
273 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
274 xudong 1.27 {
|
275 mbobra 1.34 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
|
276 mbobra 1.37 err_termAt[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
|
277 xudong 1.27 }
278 }
279
280 /* brute force method of calculating the derivative (no consideration for edges) */
281 for (i = 0; i <= nx-1; i++)
282 {
|
283 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
284 xudong 1.27 {
|
285 mbobra 1.34 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
|
286 mbobra 1.37 err_termBt[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
|
287 xudong 1.27 }
288 }
289
|
290 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
291 ignore the edges for the error terms as those arrays have been initialized to zero.
292 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
|
293 xudong 1.27 i=0;
294 for (j = 0; j <= ny-1; j++)
295 {
296 derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
297 }
298
299 i=nx-1;
300 for (j = 0; j <= ny-1; j++)
301 {
302 derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
303 }
304
305 j=0;
306 for (i = 0; i <= nx-1; i++)
307 {
308 dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
309 }
310
311 j=ny-1;
312 for (i = 0; i <= nx-1; i++)
313 {
314 xudong 1.27 dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
315 }
|
316 mbobra 1.34
317 // Calculate the sum only
|
318 xudong 1.27 for (i = 1; i <= nx-2; i++)
319 {
320 for (j = 1; j <= ny-2; j++)
321 {
322 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
323 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
324 if isnan(bt[j * nx + i]) continue;
325 if isnan(bt[(j+1) * nx + i]) continue;
326 if isnan(bt[(j-1) * nx + i]) continue;
327 if isnan(bt[j * nx + i-1]) continue;
328 if isnan(bt[j * nx + i+1]) continue;
329 if isnan(bt_err[j * nx + i]) continue;
330 if isnan(derx_bt[j * nx + i]) continue;
331 if isnan(dery_bt[j * nx + i]) continue;
332 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 */
|
333 mbobra 1.37 err += err_termBt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ))+
334 err_termAt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] )) ;
|
335 xudong 1.27 count_mask++;
336 }
337 }
338
339 *mean_derivative_btotal_ptr = (sum)/(count_mask);
340 *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask);
341 //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
342 //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
343
344 return 0;
|
345 xudong 1.1 }
346
347
348 /*===========================================*/
349 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
350
|
351 mbobra 1.37 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, float *err_termAh, float *err_termBh)
|
352 xudong 1.1 {
|
353 xudong 1.27
354 int nx = dims[0];
355 int ny = dims[1];
356 int i = 0;
357 int j = 0;
358 int count_mask = 0;
359 double sum= 0.0;
360 double err =0.0;
361 *mean_derivative_bh_ptr = 0.0;
362
363 if (nx <= 0 || ny <= 0) return 1;
364
365 /* brute force method of calculating the derivative (no consideration for edges) */
366 for (i = 1; i <= nx-2; i++)
367 {
|
368 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
369 xudong 1.27 {
|
370 mbobra 1.34 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
|
371 mbobra 1.37 err_termAh[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
|
372 xudong 1.27 }
373 }
374
375 /* brute force method of calculating the derivative (no consideration for edges) */
376 for (i = 0; i <= nx-1; i++)
377 {
|
378 mbobra 1.34 for (j = 1; j <= ny-2; j++)
379 {
380 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
|
381 mbobra 1.37 err_termBh[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
|
382 mbobra 1.34 }
|
383 xudong 1.27 }
384
|
385 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
386 ignore the edges for the error terms as those arrays have been initialized to zero.
387 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
|
388 xudong 1.27 i=0;
389 for (j = 0; j <= ny-1; j++)
390 {
391 derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
392 }
393
394 i=nx-1;
395 for (j = 0; j <= ny-1; j++)
396 {
397 derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
398 }
399
400 j=0;
401 for (i = 0; i <= nx-1; i++)
402 {
403 dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
404 }
405
406 j=ny-1;
407 for (i = 0; i <= nx-1; i++)
408 {
409 xudong 1.27 dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
410 }
411
412
413 for (i = 0; i <= nx-1; i++)
414 {
415 for (j = 0; j <= ny-1; j++)
416 {
417 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
418 if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
419 if isnan(bh[j * nx + i]) continue;
420 if isnan(bh[(j+1) * nx + i]) continue;
421 if isnan(bh[(j-1) * nx + i]) continue;
422 if isnan(bh[j * nx + i-1]) continue;
423 if isnan(bh[j * nx + i+1]) continue;
424 if isnan(bh_err[j * nx + i]) continue;
425 if isnan(derx_bh[j * nx + i]) continue;
426 if isnan(dery_bh[j * nx + i]) continue;
427 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 */
|
428 mbobra 1.37 err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ))+
429 err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] )) ;
|
430 xudong 1.27 count_mask++;
431 }
432 }
433
434 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
435 *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
436 //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
437 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
438
439 return 0;
|
440 xudong 1.1 }
441
442 /*===========================================*/
443 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
444
|
445 mbobra 1.37 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, float *err_termA, float *err_termB)
|
446 xudong 1.1 {
|
447 xudong 1.27
448 int nx = dims[0];
449 int ny = dims[1];
450 int i = 0;
451 int j = 0;
452 int count_mask = 0;
|
453 mbobra 1.34 double sum = 0.0;
|
454 xudong 1.27 double err = 0.0;
|
455 mbobra 1.34 *mean_derivative_bz_ptr = 0.0;
|
456 xudong 1.27
|
457 mbobra 1.34 if (nx <= 0 || ny <= 0) return 1;
|
458 xudong 1.27
459 /* brute force method of calculating the derivative (no consideration for edges) */
460 for (i = 1; i <= nx-2; i++)
461 {
|
462 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
463 xudong 1.27 {
|
464 mbobra 1.34 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
|
465 mbobra 1.37 err_termA[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
|
466 xudong 1.27 }
467 }
468
469 /* brute force method of calculating the derivative (no consideration for edges) */
470 for (i = 0; i <= nx-1; i++)
471 {
|
472 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
473 xudong 1.27 {
|
474 mbobra 1.34 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
|
475 mbobra 1.37 err_termB[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
|
476 xudong 1.27 }
477 }
478
|
479 mbobra 1.34 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
|
480 mbobra 1.35 ignore the edges for the error terms as those arrays have been initialized to zero.
|
481 mbobra 1.34 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
|
482 xudong 1.27 i=0;
483 for (j = 0; j <= ny-1; j++)
484 {
485 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
486 }
487
488 i=nx-1;
489 for (j = 0; j <= ny-1; j++)
490 {
491 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
492 }
493
494 j=0;
495 for (i = 0; i <= nx-1; i++)
496 {
497 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
498 }
499
500 j=ny-1;
501 for (i = 0; i <= nx-1; i++)
502 {
503 xudong 1.27 dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
504 }
505
506
507 for (i = 0; i <= nx-1; i++)
508 {
509 for (j = 0; j <= ny-1; j++)
510 {
511 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
512 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
513 if isnan(bz[j * nx + i]) continue;
514 if isnan(bz[(j+1) * nx + i]) continue;
515 if isnan(bz[(j-1) * nx + i]) continue;
516 if isnan(bz[j * nx + i-1]) continue;
517 if isnan(bz[j * nx + i+1]) continue;
518 if isnan(bz_err[j * nx + i]) continue;
519 if isnan(derx_bz[j * nx + i]) continue;
520 if isnan(dery_bz[j * nx + i]) continue;
|
521 mbobra 1.34 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 */
|
522 mbobra 1.37 err += err_termB[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) +
523 err_termA[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) ;
|
524 xudong 1.27 count_mask++;
525 }
526 }
527
|
528 mbobra 1.34 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
|
529 xudong 1.27 *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
530 //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
531 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
532
|
533 xudong 1.1 return 0;
534 }
535
536 /*===========================================*/
537 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
538
539 // In discretized space like data pixels,
540 // the current (or curl of B) is calculated as
541 // the integration of the field Bx and By along
542 // the circumference of the data pixel divided by the area of the pixel.
543 // One form of differencing (a word for the differential operator
|
544 xudong 1.27 // in the discretized space) of the curl is expressed as
|
545 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
546 // +dy * (By(i+1,j)+By(i,j)) / 2
547 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
|
548 xudong 1.27 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
549 //
|
550 xudong 1.1 //
551 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
552 // one must perform the following unit conversions:
|
553 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
554 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
|
555 xudong 1.27 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
|
556 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
|
557 mbobra 1.4 //
|
558 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
559 // In that case, we would have the following:
560 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
561 // jz * (35.0)
562 //
563 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
|
564 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)
565 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
|
566 xudong 1.1
567
|
568 xudong 1.27 // Comment out random number generator, which can only run on solar3
|
569 mbobra 1.33 // int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
|
570 xudong 1.27 // int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
|
571 mbobra 1.9 // float *noiseby, float *noisebz)
572
573 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
|
574 mbobra 1.33 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)
|
575 mbobra 1.9
|
576 xudong 1.1
|
577 xudong 1.27 {
578 int nx = dims[0];
579 int ny = dims[1];
580 int i = 0;
581 int j = 0;
582 int count_mask = 0;
583
584 if (nx <= 0 || ny <= 0) return 1;
585
586 /* Calculate the derivative*/
587 /* brute force method of calculating the derivative (no consideration for edges) */
588
589 for (i = 1; i <= nx-2; i++)
590 {
|
591 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
592 xudong 1.27 {
|
593 mbobra 1.34 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
594 err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);
|
595 xudong 1.27 }
596 }
597
598 for (i = 0; i <= nx-1; i++)
599 {
|
600 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
601 xudong 1.27 {
|
602 mbobra 1.34 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
603 err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
|
604 xudong 1.27 }
605 }
|
606 mbobra 1.33
|
607 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
608 ignore the edges for the error terms as those arrays have been initialized to zero.
609 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
610
|
611 xudong 1.27 i=0;
612 for (j = 0; j <= ny-1; j++)
613 {
|
614 mbobra 1.33 derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
|
615 xudong 1.27 }
616
617 i=nx-1;
618 for (j = 0; j <= ny-1; j++)
619 {
|
620 mbobra 1.33 derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
|
621 xudong 1.27 }
622
623 j=0;
624 for (i = 0; i <= nx-1; i++)
625 {
|
626 mbobra 1.33 dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
|
627 xudong 1.27 }
628
629 j=ny-1;
630 for (i = 0; i <= nx-1; i++)
631 {
|
632 mbobra 1.33 dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
|
633 xudong 1.27 }
|
634 mbobra 1.33
|
635 xudong 1.27
|
636 mbobra 1.33 for (i = 0; i <= nx-1; i++)
|
637 xudong 1.27 {
|
638 mbobra 1.33 for (j = 0; j <= ny-1; j++)
|
639 xudong 1.27 {
|
640 mbobra 1.33 // calculate jz at all points
|
641 xudong 1.27 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
|
642 mbobra 1.33 jz_err[j * nx + i] = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;
|
643 xudong 1.27 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
644 mbobra 1.33 count_mask++;
|
645 xudong 1.27 }
646 }
647 return 0;
648 }
|
649 mbobra 1.33
|
650 mbobra 1.5 /*===========================================*/
651
|
652 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
|
653 xudong 1.27 // Compute mean and total current on Jz array.
|
654 mbobra 1.6
|
655 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,
|
656 xudong 1.27 float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
|
657 mbobra 1.9 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
|
658 mbobra 1.5
659 {
|
660 xudong 1.27
661 int nx = dims[0];
662 int ny = dims[1];
663 int i = 0;
664 int j = 0;
665 int count_mask = 0;
|
666 mbobra 1.15 double curl = 0.0;
|
667 xudong 1.27 double us_i = 0.0;
668 double err = 0.0;
669
|
670 mbobra 1.5 if (nx <= 0 || ny <= 0) return 1;
|
671 xudong 1.27
672 /* 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*/
673 for (i = 0; i <= nx-1; i++)
674 {
675 for (j = 0; j <= ny-1; j++)
676 {
677 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
678 if isnan(derx[j * nx + i]) continue;
679 if isnan(dery[j * nx + i]) continue;
680 if isnan(jz[j * nx + i]) continue;
681 curl += (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
682 us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */
683 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
684 count_mask++;
685 }
686 }
687
688 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
689 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
690 *mean_jz_err_ptr = (sqrt(err)/count_mask)*((1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.)); // error in the quantity MEANJZD
691
692 xudong 1.27 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
693 *us_i_err_ptr = (sqrt(err))*((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
694
695 //printf("MEANJZD=%f\n",*mean_jz_ptr);
696 //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
697
698 //printf("TOTUSJZ=%g\n",*us_i_ptr);
699 //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
700
|
701 xudong 1.1 return 0;
|
702 xudong 1.27
|
703 xudong 1.1 }
704
|
705 mbobra 1.5 /*===========================================*/
706
707 /* Example function 10: Twist Parameter, alpha */
|
708 xudong 1.1
|
709 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
|
710 mbobra 1.19 // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
|
711 xudong 1.27
712 // numerator = sum of all Jz*Bz
713 // denominator = sum of Bz*Bz
714 // alpha = numerator/denominator
|
715 mbobra 1.6
|
716 mbobra 1.5 // The units of alpha are in 1/Mm
|
717 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
718 //
|
719 xudong 1.27 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
|
720 xudong 1.1 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
721 // = 1/Mm
722
|
723 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)
|
724 mbobra 1.5
|
725 xudong 1.1 {
|
726 xudong 1.27 int nx = dims[0];
727 int ny = dims[1];
728 int i = 0;
729 int j = 0;
|
730 mbobra 1.19 double alpha_total = 0.0;
|
731 xudong 1.27 double C = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
732 double total = 0.0;
733 double A = 0.0;
734 double B = 0.0;
735
|
736 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
737 xudong 1.27 for (i = 1; i < nx-1; i++)
738 {
739 for (j = 1; j < ny-1; j++)
740 {
741 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
742 if isnan(jz[j * nx + i]) continue;
743 if isnan(bz[j * nx + i]) continue;
744 if (jz[j * nx + i] == 0.0) continue;
745 if (bz[j * nx + i] == 0.0) continue;
746 A += jz[j*nx+i]*bz[j*nx+i];
747 B += bz[j*nx+i]*bz[j*nx+i];
748 }
749 }
750
751 for (i = 1; i < nx-1; i++)
752 {
753 for (j = 1; j < ny-1; j++)
754 {
755 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
756 if isnan(jz[j * nx + i]) continue;
757 if isnan(bz[j * nx + i]) continue;
758 xudong 1.27 if (jz[j * nx + i] == 0.0) continue;
759 if (bz[j * nx + i] == 0.0) continue;
760 total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
761 }
762 }
763
764 /* Determine the absolute value of alpha. The units for alpha are 1/Mm */
765 alpha_total = ((A/B)*C);
766 *mean_alpha_ptr = alpha_total;
767 *mean_alpha_err_ptr = (C/B)*(sqrt(total));
768
|
769 xudong 1.1 return 0;
770 }
771
772 /*===========================================*/
|
773 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
|
774 xudong 1.1
775 // The current helicity is defined as Bz*Jz and the units are G^2 / m
776 // The units of Jz are in G/pix; the units of Bz are in G.
|
777 xudong 1.27 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
778 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
|
779 mbobra 1.9 // = G^2 / m.
|
780 xudong 1.1
|
781 xudong 1.27 int computeHelicity(float *jz_err, float *jz_rms_err, float *bz_err, float *bz, int *dims, float *jz, float *mean_ih_ptr,
782 float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
|
783 mbobra 1.9 float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
784 xudong 1.1
785 {
|
786 xudong 1.27
787 int nx = dims[0];
788 int ny = dims[1];
789 int i = 0;
790 int j = 0;
791 int count_mask = 0;
|
792 mbobra 1.20 double sum = 0.0;
793 double sum2 = 0.0;
794 double err = 0.0;
|
795 xudong 1.1
796 if (nx <= 0 || ny <= 0) return 1;
|
797 xudong 1.27
798 for (i = 0; i < nx; i++)
|
799 xudong 1.1 {
|
800 xudong 1.27 for (j = 0; j < ny; j++)
|
801 xudong 1.1 {
|
802 xudong 1.27 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
803 if isnan(jz[j * nx + i]) continue;
804 if isnan(bz[j * nx + i]) continue;
805 if isnan(jz_err[j * nx + i]) continue;
806 if isnan(bz_err[j * nx + i]) continue;
807 if (bz[j * nx + i] == 0.0) continue;
808 if (jz[j * nx + i] == 0.0) continue;
809 sum += (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
810 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
811 err += (jz_err[j * nx + i]*jz_err[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]) + (bz_err[j * nx + i]*bz_err[j * nx + i]*jz[j * nx + i]*jz[j * nx + i]);
812 count_mask++;
813 }
814 }
815
816 *mean_ih_ptr = sum/count_mask ; /* Units are G^2 / m ; keyword is MEANJZH */
|
817 mbobra 1.9 *total_us_ih_ptr = sum2 ; /* Units are G^2 / m ; keyword is TOTUSJH */
818 *total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
|
819 xudong 1.27
820 *mean_ih_err_ptr = (sqrt(err)/count_mask)*(1/cdelt1)*(rsun_obs/rsun_ref) ; // error in the quantity MEANJZH
821 *total_us_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ; // error in the quantity TOTUSJH
822 *total_abs_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ; // error in the quantity ABSNJZH
823
824 //printf("MEANJZH=%f\n",*mean_ih_ptr);
825 //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
826
827 //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
828 //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
829
830 //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
831 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
832
|
833 xudong 1.1 return 0;
834 }
835
836 /*===========================================*/
|
837 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
|
838 xudong 1.1
839 // The Sum of the Absolute Value per polarity is defined as the following:
|
840 mbobra 1.36 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
|
841 xudong 1.1 // The units of jz are in G/pix. In this case, we would have the following:
842 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
843 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
|
844 mbobra 1.9 //
845 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
|
846 xudong 1.1
|
847 xudong 1.27 int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
|
848 mbobra 1.3 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
849 xudong 1.1
|
850 xudong 1.27 {
851 int nx = dims[0];
852 int ny = dims[1];
853 int i=0;
854 int j=0;
855 int count_mask=0;
|
856 mbobra 1.32 double sum1=0.0;
|
857 xudong 1.27 double sum2=0.0;
858 double err=0.0;
|
859 mbobra 1.33 *totaljzptr=0.0;
|
860 xudong 1.27
|
861 mbobra 1.33 if (nx <= 0 || ny <= 0) return 1;
|
862 xudong 1.27
|
863 mbobra 1.33 for (i = 0; i < nx; i++)
|
864 xudong 1.27 {
|
865 mbobra 1.33 for (j = 0; j < ny; j++)
|
866 xudong 1.27 {
867 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
868 if isnan(bz[j * nx + i]) continue;
869 if isnan(jz[j * nx + i]) continue;
870 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
871 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
872 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
873 count_mask++;
874 }
875 }
|
876 xudong 1.1
|
877 mbobra 1.32 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are Amperes per arcsecond */
|
878 xudong 1.27 *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
879 //printf("SAVNCPP=%g\n",*totaljzptr);
880 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
881
|
882 mbobra 1.33 return 0;
|
883 xudong 1.1 }
884
885 /*===========================================*/
|
886 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
887 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
|
888 xudong 1.27 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
889 // the integral is over B^2 dV and the 8*PI is divided at the end.
|
890 xudong 1.1 //
891 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
892 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
|
893 mbobra 1.9 // erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
894 // = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
895 // = erg/cm^3*(1.30501e15)
|
896 xudong 1.1 // = erg/cm(1/pix^2)
897
|
898 xudong 1.27 int computeFreeEnergy(float *bx_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
899 float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
|
900 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
901
902 {
|
903 xudong 1.27 int nx = dims[0];
904 int ny = dims[1];
905 int i = 0;
906 int j = 0;
907 int count_mask = 0;
|
908 mbobra 1.31 double sum = 0.0;
|
909 xudong 1.27 double sum1 = 0.0;
910 double err = 0.0;
911 *totpotptr = 0.0;
|
912 mbobra 1.31 *meanpotptr = 0.0;
|
913 xudong 1.27
|
914 mbobra 1.31 if (nx <= 0 || ny <= 0) return 1;
|
915 xudong 1.27
|
916 mbobra 1.31 for (i = 0; i < nx; i++)
|
917 xudong 1.27 {
|
918 mbobra 1.31 for (j = 0; j < ny; j++)
|
919 xudong 1.27 {
920 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
921 if isnan(bx[j * nx + i]) continue;
922 if isnan(by[j * nx + i]) continue;
923 sum += ( ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
924 sum1 += ( ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) );
925 err += 4.0*(bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])*(bx_err[j * nx + i]*bx_err[j * nx + i]) +
926 4.0*(by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])*(by_err[j * nx + i]*by_err[j * nx + i]);
927 count_mask++;
928 }
929 }
930
931 /* Units of meanpotptr are ergs per centimeter */
|
932 mbobra 1.20 *meanpotptr = (sum1) / (count_mask*8.*PI) ; /* Units are ergs per cubic centimeter */
|
933 xudong 1.27 *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
934
935 /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
936 *totpotptr = (sum)/(8.*PI);
937 *totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
938
939 //printf("MEANPOT=%g\n",*meanpotptr);
940 //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
941
942 //printf("TOTPOT=%g\n",*totpotptr);
943 //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
944
|
945 mbobra 1.31 return 0;
|
946 xudong 1.1 }
947
948 /*===========================================*/
|
949 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 */
|
950 xudong 1.1
|
951 mbobra 1.20 int computeShearAngle(float *bx_err, float *by_err, float *bz_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
|
952 mbobra 1.9 float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
|
953 mbobra 1.21
954
|
955 xudong 1.27 {
956 int nx = dims[0];
957 int ny = dims[1];
958 int i = 0;
959 int j = 0;
960 float count_mask = 0;
961 float count = 0;
962 double dotproduct = 0.0;
963 double magnitude_potential = 0.0;
964 double magnitude_vector = 0.0;
965 double shear_angle = 0.0;
966 double denominator = 0.0;
967 double term1 = 0.0;
968 double term2 = 0.0;
969 double term3 = 0.0;
970 double sumsum = 0.0;
971 double err = 0.0;
972 double part1 = 0.0;
973 double part2 = 0.0;
974 double part3 = 0.0;
975 *area_w_shear_gt_45ptr = 0.0;
|
976 mbobra 1.31 *meanshear_angleptr = 0.0;
|
977 xudong 1.1
|
978 mbobra 1.31 if (nx <= 0 || ny <= 0) return 1;
|
979 xudong 1.27
|
980 mbobra 1.31 for (i = 0; i < nx; i++)
|
981 xudong 1.27 {
|
982 mbobra 1.31 for (j = 0; j < ny; j++)
|
983 xudong 1.27 {
984 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
985 if isnan(bpx[j * nx + i]) continue;
986 if isnan(bpy[j * nx + i]) continue;
987 if isnan(bpz[j * nx + i]) continue;
988 if isnan(bz[j * nx + i]) continue;
989 if isnan(bx[j * nx + i]) continue;
990 if isnan(by[j * nx + i]) continue;
991 if isnan(bx_err[j * nx + i]) continue;
992 if isnan(by_err[j * nx + i]) continue;
993 if isnan(bz_err[j * nx + i]) continue;
994
995 /* For mean 3D shear angle, percentage with shear greater than 45*/
996 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]);
997 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]));
998 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]) );
999 //printf("dotproduct=%f\n",dotproduct);
1000 //printf("magnitude_potential=%f\n",magnitude_potential);
1001 //printf("magnitude_vector=%f\n",magnitude_vector);
1002
1003 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
1004 xudong 1.27 sumsum += shear_angle;
1005 //printf("shear_angle=%f\n",shear_angle);
1006 count ++;
1007
1008 if (shear_angle > 45) count_mask ++;
1009
1010 // For the error analysis
1011
1012 term1 = bx[j * nx + i]*by[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bz[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i];
1013 term2 = bx[j * nx + i]*bx[j * nx + i]*bpy[j * nx + i] - bx[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bx[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i] - bz[j * nx + i]*by[j * nx + i]*bpz[j * nx + i];
1014 term3 = bx[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bx[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*by[j * nx + i]*bpz[j * nx + i] - by[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i];
1015
1016 part1 = bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i];
1017 part2 = bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i];
1018 part3 = bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i];
1019
1020 // denominator is squared
1021 denominator = part1*part1*part1*part2*(1.0-((part3*part3)/(part1*part2)));
1022
1023 err = (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
1024 (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
1025 xudong 1.27 (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) ;
1026
1027 }
1028 }
1029 /* For mean 3D shear angle, area with shear greater than 45*/
1030 *meanshear_angleptr = (sumsum)/(count); /* Units are degrees */
|
1031 mbobra 1.39
1032 // For the error in the mean 3D shear angle:
1033 // If count_mask is 0, then we run into a divide by zero error. In this case, set *meanshear_angle_err_ptr to NAN
1034 // If count_mask is greater than zero, then compute the error.
1035 if (count_mask == 0)
1036 *meanshear_angle_err_ptr = NAN;
1037 else
1038 *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
|
1039 xudong 1.27
1040 /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */
1041 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.0);
1042
1043 //printf("MEANSHR=%f\n",*meanshear_angleptr);
|
1044 mbobra 1.39 //printf("ERRMSHA=%f\n",*meanshear_angle_err_ptr);
|
1045 xudong 1.27 //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
|
1046 mbobra 1.39 return 0;
|
1047 xudong 1.27 }
|
1048 xudong 1.1
|
1049 xudong 1.27 /*===========================================*/
|
1050 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 */
1051 //
1052 // Note that there is a restriction on the function fsample()
1053 // If the following occurs:
1054 // nx_out > floor((ny_in-1)/scale + 1)
1055 // ny_out > floor((ny_in-1)/scale + 1),
1056 // where n*_out are the dimensions of the output array and n*_in
1057 // are the dimensions of the input array, fsample() will usually result
1058 // in a segfault (though not always, depending on how the segfault was accessed.)
1059
|
1060 xudong 1.27 int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
1061 float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
|
1062 mbobra 1.30 float *pmap, int nx1, int ny1,
1063 int scale, float *p1pad, int nxp, int nyp, float *pmapn)
1064
|
1065 mbobra 1.29 {
|
1066 xudong 1.27 int nx = dims[0];
1067 int ny = dims[1];
1068 int i = 0;
1069 int j = 0;
|
1070 mbobra 1.30 int index, index1;
|
1071 xudong 1.27 double sum = 0.0;
1072 double err = 0.0;
1073 *Rparam = 0.0;
1074 struct fresize_struct fresboxcar, fresgauss;
|
1075 mbobra 1.30 struct fint_struct fints;
|
1076 xudong 1.27 float sigma = 10.0/2.3548;
1077
|
1078 mbobra 1.30 // set up convolution kernels
|
1079 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1);
1080 init_fresize_gaussian(&fresgauss,sigma,20,1);
|
1081 mbobra 1.29
|
1082 mbobra 1.30 // =============== [STEP 1] ===============
1083 // bin the line-of-sight magnetogram down by a factor of scale
|
1084 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
|
1085 mbobra 1.29
|
1086 mbobra 1.30 // =============== [STEP 2] ===============
1087 // identify positive and negative pixels greater than +/- 150 gauss
1088 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
|
1089 xudong 1.27 for (i = 0; i < nx1; i++)
1090 {
1091 for (j = 0; j < ny1; j++)
1092 {
1093 index = j * nx1 + i;
1094 if (rim[index] > 150)
1095 p1p0[index]=1.0;
1096 else
1097 p1p0[index]=0.0;
1098 if (rim[index] < -150)
1099 p1n0[index]=1.0;
1100 else
1101 p1n0[index]=0.0;
1102 }
1103 }
|
1104 mbobra 1.30
1105 // =============== [STEP 3] ===============
1106 // smooth each of the negative and positive pixel bitmaps
|
1107 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1108 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
|
1109 mbobra 1.30
1110 // =============== [STEP 4] ===============
1111 // find the pixels for which p1p and p1n are both equal to 1.
1112 // this defines the polarity inversion line
|
1113 xudong 1.27 for (i = 0; i < nx1; i++)
1114 {
1115 for (j = 0; j < ny1; j++)
1116 {
1117 index = j * nx1 + i;
|
1118 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
|
1119 xudong 1.27 p1[index]=1.0;
1120 else
1121 p1[index]=0.0;
1122 }
1123 }
|
1124 mbobra 1.30
1125 // pad p1 with zeroes so that the gaussian colvolution in step 5
1126 // does not cut off data within hwidth of the edge
1127
1128 // step i: zero p1pad
1129 for (i = 0; i < nxp; i++)
1130 {
1131 for (j = 0; j < nyp; j++)
1132 {
1133 index = j * nxp + i;
1134 p1pad[index]=0.0;
1135 }
1136 }
1137
1138 // step ii: place p1 at the center of p1pad
1139 for (i = 0; i < nx1; i++)
1140 {
1141 for (j = 0; j < ny1; j++)
1142 {
1143 index = j * nx1 + i;
1144 index1 = (j+20) * nxp + (i+20);
1145 mbobra 1.30 p1pad[index1]=p1[index];
1146 }
1147 }
1148
1149 // =============== [STEP 5] ===============
1150 // convolve the polarity inversion line map with a gaussian
1151 // to identify the region near the plarity inversion line
1152 // the resultant array is called pmap
1153 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);
1154
1155
1156 // select out the nx1 x ny1 non-padded array within pmap
1157 for (i = 0; i < nx1; i++)
1158 {
1159 for (j = 0; j < ny1; j++)
1160 {
1161 index = j * nx1 + i;
1162 index1 = (j+20) * nxp + (i+20);
1163 pmapn[index]=pmap[index1];
1164 }
1165 }
1166 mbobra 1.30
1167 // =============== [STEP 6] ===============
1168 // the R parameter is calculated
|
1169 xudong 1.27 for (i = 0; i < nx1; i++)
1170 {
1171 for (j = 0; j < ny1; j++)
1172 {
1173 index = j * nx1 + i;
|
1174 mbobra 1.32 if isnan(pmapn[index]) continue;
1175 if isnan(rim[index]) continue;
|
1176 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
|
1177 xudong 1.27 }
1178 }
1179
1180 if (sum < 1.0)
1181 *Rparam = 0.0;
1182 else
1183 *Rparam = log10(sum);
|
1184 mbobra 1.30
|
1185 mbobra 1.32 //printf("R_VALUE=%f\n",*Rparam);
1186
|
1187 xudong 1.27 free_fresize(&fresboxcar);
1188 free_fresize(&fresgauss);
|
1189 mbobra 1.30
|
1190 xudong 1.27 return 0;
|
1191 mbobra 1.30
|
1192 xudong 1.1 }
1193
|
1194 mbobra 1.31 /*===========================================*/
1195 /* Example function 16: Lorentz force as defined in Fisher, 2012 */
1196 //
1197 // This calculation is adapted from Xudong's code
1198 // at /proj/cgem/lorentz/apps/lorentz.c
1199
1200 int computeLorentz(float *bx, float *by, float *bz, float *fx, float *fy, float *fz, int *dims,
1201 float *totfx_ptr, float *totfy_ptr, float *totfz_ptr, float *totbsq_ptr,
1202 float *epsx_ptr, float *epsy_ptr, float *epsz_ptr, int *mask, int *bitmask,
1203 float cdelt1, double rsun_ref, double rsun_obs)
1204
1205 {
1206
1207 int nx = dims[0];
1208 int ny = dims[1];
1209 int nxny = nx*ny;
1210 int j = 0;
1211 int index;
1212 double totfx = 0, totfy = 0, totfz = 0;
1213 double bsq = 0, totbsq = 0;
1214 double epsx = 0, epsy = 0, epsz = 0;
1215 mbobra 1.31 double area = cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
1216 double k_h = -1.0 * area / (4. * PI) / 1.0e20;
1217 double k_z = area / (8. * PI) / 1.0e20;
1218
1219 if (nx <= 0 || ny <= 0) return 1;
1220
1221 for (int i = 0; i < nxny; i++)
1222 {
1223 if ( mask[i] < 70 || bitmask[i] < 30 ) continue;
1224 if isnan(bx[i]) continue;
1225 if isnan(by[i]) continue;
1226 if isnan(bz[i]) continue;
|
1227 mbobra 1.32 fx[i] = bx[i] * bz[i] * k_h;
1228 fy[i] = by[i] * bz[i] * k_h;
|
1229 mbobra 1.31 fz[i] = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
1230 bsq = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
1231 totfx += fx[i]; totfy += fy[i]; totfz += fz[i];
1232 totbsq += bsq;
1233 }
1234
1235 *totfx_ptr = totfx;
1236 *totfy_ptr = totfy;
1237 *totfz_ptr = totfz;
1238 *totbsq_ptr = totbsq;
1239 *epsx_ptr = (totfx / k_h) / totbsq;
1240 *epsy_ptr = (totfy / k_h) / totbsq;
1241 *epsz_ptr = (totfz / k_z) / totbsq;
1242
|
1243 mbobra 1.32 //printf("TOTBSQ=%f\n",*totbsq_ptr);
1244
|
1245 mbobra 1.31 return 0;
1246
1247 }
1248
|
1249 mbobra 1.38 /*===========================================*/
1250
1251 /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field */
1252
1253 // To compute the unsigned flux, we simply calculate
1254 // flux = surface integral [(vector LOS) dot (normal vector)],
1255 // = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
1256 // However, since the field is radial, we will assume cos theta = 1.
1257 // Therefore the pixels only need to be corrected for the projection.
1258
1259 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
1260 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
1261 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
1262 // =Gauss*cm^2
1263
|
1264 mbobra 1.40 int computeAbsFlux_los(float *los, int *dims, float *absFlux_los,
1265 float *mean_vf_los_ptr, float *count_mask_los_ptr,
|
1266 mbobra 1.38 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
1267
1268 {
1269
1270 int nx = dims[0];
1271 int ny = dims[1];
1272 int i = 0;
1273 int j = 0;
|
1274 mbobra 1.40 int count_mask_los = 0;
|
1275 mbobra 1.38 double sum = 0.0;
|
1276 mbobra 1.40 *absFlux_los = 0.0;
1277 *mean_vf_los_ptr = 0.0;
|
1278 mbobra 1.38
1279
1280 if (nx <= 0 || ny <= 0) return 1;
1281
1282 for (i = 0; i < nx; i++)
1283 {
1284 for (j = 0; j < ny; j++)
1285 {
1286 if ( bitmask[j * nx + i] < 30 ) continue;
|
1287 mbobra 1.43 if isnan(los[j * nx + i]) continue;
|
1288 mbobra 1.38 sum += (fabs(los[j * nx + i]));
|
1289 mbobra 1.40 count_mask_los++;
|
1290 mbobra 1.38 }
1291 }
1292
|
1293 mbobra 1.40 *mean_vf_los_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
1294 *count_mask_los_ptr = count_mask_los;
1295
|
1296 mbobra 1.42 //printf("USFLUXL=%f\n",*mean_vf_los_ptr);
1297 //printf("CMASKL=%f\n",*count_mask_los_ptr);
|
1298 mbobra 1.38
1299 return 0;
1300 }
1301
1302 /*===========================================*/
1303 /* Example function 18: Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
1304
1305 int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
1306 {
1307
1308 int nx = dims[0];
1309 int ny = dims[1];
1310 int i = 0;
1311 int j = 0;
1312 int count_mask = 0;
1313 double sum = 0.0;
1314 *mean_derivative_los_ptr = 0.0;
1315
1316 if (nx <= 0 || ny <= 0) return 1;
1317
1318 /* brute force method of calculating the derivative (no consideration for edges) */
1319 mbobra 1.38 for (i = 1; i <= nx-2; i++)
1320 {
1321 for (j = 0; j <= ny-1; j++)
1322 {
1323 derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
1324 }
1325 }
1326
1327 /* brute force method of calculating the derivative (no consideration for edges) */
1328 for (i = 0; i <= nx-1; i++)
1329 {
1330 for (j = 1; j <= ny-2; j++)
1331 {
1332 dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
1333 }
1334 }
1335
1336 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
1337 ignore the edges for the error terms as those arrays have been initialized to zero.
1338 this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
1339 i=0;
1340 mbobra 1.38 for (j = 0; j <= ny-1; j++)
1341 {
1342 derx_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
1343 }
1344
1345 i=nx-1;
1346 for (j = 0; j <= ny-1; j++)
1347 {
1348 derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
1349 }
1350
1351 j=0;
1352 for (i = 0; i <= nx-1; i++)
1353 {
1354 dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5;
1355 }
1356
1357 j=ny-1;
1358 for (i = 0; i <= nx-1; i++)
1359 {
1360 dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
1361 mbobra 1.38 }
1362
1363
1364 for (i = 0; i <= nx-1; i++)
1365 {
1366 for (j = 0; j <= ny-1; j++)
1367 {
1368 if ( bitmask[j * nx + i] < 30 ) continue;
1369 if isnan(los[j * nx + i]) continue;
1370 if isnan(los[(j+1) * nx + i]) continue;
1371 if isnan(los[(j-1) * nx + i]) continue;
1372 if isnan(los[j * nx + i-1]) continue;
1373 if isnan(los[j * nx + i+1]) continue;
1374 if isnan(derx_los[j * nx + i]) continue;
1375 if isnan(dery_los[j * nx + i]) continue;
1376 sum += sqrt( derx_los[j * nx + i]*derx_los[j * nx + i] + dery_los[j * nx + i]*dery_los[j * nx + i] ); /* Units of Gauss */
1377 count_mask++;
1378 }
1379 }
1380
1381 *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
|
1382 mbobra 1.40
|
1383 mbobra 1.42 //printf("MEANGBL=%f\n",*mean_derivative_los_ptr);
|
1384 mbobra 1.38
1385 return 0;
1386 }
1387
|
1388 xudong 1.1 /*==================KEIJI'S CODE =========================*/
1389
1390 // #include <omp.h>
1391 #include <math.h>
1392
1393 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
1394 {
|
1395 xudong 1.27 /* local workings */
1396 int inx, iny, i, j, n;
1397 /* local array */
1398 float *pfpot, *rdist;
1399 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
1400 rdist=(float *)malloc(sizeof(float) *nnx*nny);
1401 float *bztmp;
1402 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
1403 /* make nan */
1404 // unsigned long long llnan = 0x7ff0000000000000;
1405 // float NAN = (float)(llnan);
1406
1407 // #pragma omp parallel for private (inx)
1408 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
1409 // #pragma omp parallel for private (inx)
1410 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
1411 // #pragma omp parallel for private (inx)
1412 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
1413 // #pragma omp parallel for private (inx)
1414 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
1415 // #pragma omp parallel for private (inx)
1416 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1417 {
1418 float val0 = bz[nnx*iny + inx];
1419 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
1420 }}
1421
1422 // dz is the monopole depth
1423 float dz = 0.001;
1424
1425 // #pragma omp parallel for private (inx)
1426 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1427 {
1428 float rdd, rdd1, rdd2;
1429 float r;
1430 rdd1 = (float)(inx);
1431 rdd2 = (float)(iny);
1432 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
1433 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
1434 }}
1435
1436 int iwindow;
1437 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
1438 float rwindow;
1439 rwindow = (float)(iwindow);
1440 rwindow = rwindow * rwindow + 0.01; // must be of square
1441
1442 rwindow = 1.0e2; // limit the window size to be 10.
1443
1444 rwindow = sqrt(rwindow);
1445 iwindow = (int)(rwindow);
1446
1447 // #pragma omp parallel for private(inx)
1448 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
|
1449 xudong 1.1 {
|
1450 xudong 1.27 float val0 = bz[nnx*iny + inx];
1451 if (isnan(val0))
1452 {
1453 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
1454 }
1455 else
1456 {
1457 float sum;
1458 sum = 0.0;
1459 int j2, i2;
1460 int j2s, j2e, i2s, i2e;
1461 j2s = iny - iwindow;
1462 j2e = iny + iwindow;
1463 if (j2s < 0){j2s = 0;}
1464 if (j2e > nny){j2e = nny;}
1465 i2s = inx - iwindow;
1466 i2e = inx + iwindow;
1467 if (i2s < 0){i2s = 0;}
|
1468 mbobra 1.37 if (i2e > nnx){i2e = nnx;}
1469
1470 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1471 {
|
1472 xudong 1.27 float val1 = bztmp[nnx*j2 + i2];
|
1473 mbobra 1.37 float rr, r1, r2;
1474 // r1 = (float)(i2 - inx);
1475 // r2 = (float)(j2 - iny);
1476 // rr = r1*r1 + r2*r2;
1477 // if (rr < rwindow)
1478 // {
|
1479 xudong 1.27 int di, dj;
1480 di = abs(i2 - inx);
1481 dj = abs(j2 - iny);
1482 sum = sum + val1 * rdist[nnx * dj + di] * dz;
|
1483 mbobra 1.37 // }
|
1484 xudong 1.27 } }
1485 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1486 }
1487 } } // end of OpenMP parallelism
1488
1489 // #pragma omp parallel for private(inx)
1490 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
|
1491 xudong 1.1 {
|
1492 xudong 1.27 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1493 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1494 } } // end of OpenMP parallelism
|
1495 mbobra 1.37
|
1496 xudong 1.27 free(rdist);
1497 free(pfpot);
1498 free(bztmp);
|
1499 xudong 1.1 } // end of void func. greenpot
1500
1501
1502 /*===========END OF KEIJI'S CODE =========================*/
|
1503 mbobra 1.14
1504 char *sw_functions_version() // Returns CVS version of sw_functions.c
1505 {
|
1506 xudong 1.27 return strdup("$Id");
|
1507 mbobra 1.14 }
1508
|
1509 xudong 1.1 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps