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