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