1 mbobra 1.37
|
2 xudong 1.1 /*===========================================
|
3 xudong 1.27
4 The following 14 functions calculate the following spaceweather indices:
5
6 USFLUX Total unsigned flux in Maxwells
7 MEANGAM Mean inclination angle, gamma, in degrees
8 MEANGBT Mean value of the total field gradient, in Gauss/Mm
9 MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
10 MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
11 MEANJZD Mean vertical current density, in mA/m2
12 TOTUSJZ Total unsigned vertical current, in Amperes
13 MEANALP Mean twist parameter, alpha, in 1/Mm
14 MEANJZH Mean current helicity in G2/m
15 TOTUSJH Total unsigned current helicity in G2/m
16 ABSNJZH Absolute value of the net current helicity in G2/m
17 SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
18 MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
19 TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
20 MEANSHR Mean shear angle (measured using Btotal) in degrees
|
21 mbobra 1.29 R_VALUE Karel Schrijver's R parameter
|
22 xudong 1.27
23 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
24 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
25 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
26 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
27 contain nearest-neighbor bitmaps).
28
29 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
30 and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:
31
32 For conf_disambig:
33 50 : not all solutions agree (weak field method applied)
34 60 : not all solutions agree (weak field + annealed)
35 90 : all solutions agree (strong field + annealed)
36 0 : not disambiguated
37
38 For bitmap:
39 1 : weak field outside smooth bounding curve
40 2 : strong field outside smooth bounding curve
41 33 : weak field inside smooth bounding curve
42 34 : strong field inside smooth bounding curve
43 xudong 1.27
44 Written by Monica Bobra 15 August 2012
45 Potential Field code (appended) written by Keiji Hayashi
46 Error analysis modification 21 October 2013
47
48 ===========================================*/
|
49 xudong 1.1 #include <math.h>
|
50 mbobra 1.9 #include <mkl.h>
|
51 xudong 1.1
|
52 mbobra 1.9 #define PI (M_PI)
|
53 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
54
55 /*===========================================*/
56
57 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
58
|
59 xudong 1.27 // To compute the unsigned flux, we simply calculate
|
60 xudong 1.1 // flux = surface integral [(vector Bz) dot (normal vector)],
61 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
62 // However, since the field is radial, we will assume cos theta = 1.
63 // Therefore the pixels only need to be corrected for the projection.
64
65 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
66 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
67 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
|
68 mbobra 1.20 // =Gauss*cm^2
|
69 xudong 1.1
|
70 xudong 1.27 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
71 float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
|
72 mbobra 1.9 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
73 xudong 1.1
74 {
|
75 xudong 1.27
|
76 mbobra 1.14 int nx = dims[0];
77 int ny = dims[1];
|
78 mbobra 1.15 int i = 0;
79 int j = 0;
80 int count_mask = 0;
81 double sum = 0.0;
82 double err = 0.0;
|
83 mbobra 1.14 *absFlux = 0.0;
84 *mean_vf_ptr = 0.0;
|
85 xudong 1.27
86
|
87 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
88 xudong 1.27
89 for (i = 0; i < nx; i++)
|
90 xudong 1.1 {
|
91 mbobra 1.34 for (j = 0; j < ny; j++)
92 {
|
93 xudong 1.27 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
94 if isnan(bz[j * nx + i]) continue;
95 sum += (fabs(bz[j * nx + i]));
96 err += bz_err[j * nx + i]*bz_err[j * nx + i];
97 count_mask++;
|
98 mbobra 1.34 }
|
99 xudong 1.1 }
|
100 xudong 1.27
101 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
102 *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
103 *count_mask_ptr = count_mask;
104 return 0;
|
105 xudong 1.1 }
106
107 /*===========================================*/
|
108 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
|
109 xudong 1.1 // Native units of Bh are Gauss
110
|
111 xudong 1.27 int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
|
112 mbobra 1.3 float *mean_hf_ptr, int *mask, int *bitmask)
|
113 xudong 1.1
114 {
|
115 xudong 1.27
|
116 mbobra 1.14 int nx = dims[0];
117 int ny = dims[1];
|
118 mbobra 1.15 int i = 0;
|
119 xudong 1.27 int j = 0;
|
120 mbobra 1.15 int count_mask = 0;
|
121 xudong 1.27 double sum = 0.0;
|
122 mbobra 1.9 *mean_hf_ptr = 0.0;
|
123 xudong 1.27
|
124 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
125 xudong 1.27
126 for (i = 0; i < nx; i++)
127 {
|
128 mbobra 1.5 for (j = 0; j < ny; j++)
|
129 xudong 1.27 {
130 if isnan(bx[j * nx + i])
131 {
132 bh[j * nx + i] = NAN;
133 bh_err[j * nx + i] = NAN;
134 continue;
135 }
136 if isnan(by[j * nx + i])
137 {
138 bh[j * nx + i] = NAN;
139 bh_err[j * nx + i] = NAN;
140 continue;
141 }
142 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
143 sum += bh[j * nx + i];
144 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];
145 count_mask++;
146 }
147 }
148
|
149 xudong 1.1 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
|
150 xudong 1.27
|
151 xudong 1.1 return 0;
152 }
153
154 /*===========================================*/
155 /* Example function 3: Calculate Gamma in units of degrees */
156 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
|
157 xudong 1.27 //
|
158 mbobra 1.20 // Error analysis calculations are done in radians (since derivatives are only true in units of radians),
159 // and multiplied by (180./PI) at the end for consistency in units.
|
160 xudong 1.1
|
161 mbobra 1.9 int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
162 float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
|
163 xudong 1.1 {
|
164 mbobra 1.14 int nx = dims[0];
165 int ny = dims[1];
|
166 mbobra 1.15 int i = 0;
167 int j = 0;
168 int count_mask = 0;
169 double sum = 0.0;
170 double err = 0.0;
171 *mean_gamma_ptr = 0.0;
|
172 xudong 1.27
|
173 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
174 xudong 1.27
175 for (i = 0; i < nx; i++)
176 {
177 for (j = 0; j < ny; j++)
178 {
179 if (bh[j * nx + i] > 100)
180 {
181 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
182 if isnan(bz[j * nx + i]) continue;
183 if isnan(bz_err[j * nx + i]) continue;
184 if isnan(bh_err[j * nx + i]) continue;
185 if isnan(bh[j * nx + i]) continue;
186 if (bz[j * nx + i] == 0) continue;
187 sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
188 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])))) *
189 ( ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) +
190 ((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])) );
191 count_mask++;
192 }
193 }
194 }
195 xudong 1.27
196 *mean_gamma_ptr = sum/count_mask;
197 *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
198 //printf("MEANGAM=%f\n",*mean_gamma_ptr);
199 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
200 return 0;
|
201 xudong 1.1 }
202
203 /*===========================================*/
204 /* Example function 4: Calculate B_Total*/
205 // Native units of B_Total are in gauss
206
|
207 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)
|
208 xudong 1.1 {
|
209 xudong 1.27
|
210 mbobra 1.14 int nx = dims[0];
211 int ny = dims[1];
|
212 mbobra 1.15 int i = 0;
213 int j = 0;
214 int count_mask = 0;
|
215 xudong 1.1
216 if (nx <= 0 || ny <= 0) return 1;
|
217 xudong 1.27
218 for (i = 0; i < nx; i++)
219 {
220 for (j = 0; j < ny; j++)
221 {
222 if isnan(bx[j * nx + i])
223 {
224 bt[j * nx + i] = NAN;
225 bt_err[j * nx + i] = NAN;
226 continue;
227 }
228 if isnan(by[j * nx + i])
229 {
230 bt[j * nx + i] = NAN;
231 bt_err[j * nx + i] = NAN;
232 continue;
233 }
234 if isnan(bz[j * nx + i])
235 {
236 bt[j * nx + i] = NAN;
237 bt_err[j * nx + i] = NAN;
238 xudong 1.27 continue;
239 }
240 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]);
241 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];
242 }
243 }
244 return 0;
|
245 xudong 1.1 }
246
247 /*===========================================*/
248 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
249
|
250 mbobra 1.37 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_termAt, float *err_termBt)
|
251 xudong 1.1 {
|
252 xudong 1.27
|
253 mbobra 1.14 int nx = dims[0];
254 int ny = dims[1];
|
255 mbobra 1.15 int i = 0;
256 int j = 0;
257 int count_mask = 0;
|
258 xudong 1.27 double sum = 0.0;
|
259 mbobra 1.15 double err = 0.0;
|
260 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
|
261 xudong 1.27
|
262 xudong 1.1 if (nx <= 0 || ny <= 0) return 1;
|
263 xudong 1.27
264 /* brute force method of calculating the derivative (no consideration for edges) */
265 for (i = 1; i <= nx-2; i++)
266 {
|
267 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
268 xudong 1.27 {
|
269 mbobra 1.34 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
|
270 mbobra 1.37 err_termAt[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
|
271 xudong 1.27 }
272 }
273
274 /* brute force method of calculating the derivative (no consideration for edges) */
275 for (i = 0; i <= nx-1; i++)
276 {
|
277 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
278 xudong 1.27 {
|
279 mbobra 1.34 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
|
280 mbobra 1.37 err_termBt[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
|
281 xudong 1.27 }
282 }
283
|
284 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
285 ignore the edges for the error terms as those arrays have been initialized to zero.
286 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.*/
|
287 xudong 1.27 i=0;
288 for (j = 0; j <= ny-1; j++)
289 {
290 derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
291 }
292
293 i=nx-1;
294 for (j = 0; j <= ny-1; j++)
295 {
296 derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
297 }
298
299 j=0;
300 for (i = 0; i <= nx-1; i++)
301 {
302 dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
303 }
304
305 j=ny-1;
306 for (i = 0; i <= nx-1; i++)
307 {
308 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;
309 }
|
310 mbobra 1.34
311 // Calculate the sum only
|
312 xudong 1.27 for (i = 1; i <= nx-2; i++)
313 {
314 for (j = 1; j <= ny-2; j++)
315 {
316 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
317 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
318 if isnan(bt[j * nx + i]) continue;
319 if isnan(bt[(j+1) * nx + i]) continue;
320 if isnan(bt[(j-1) * nx + i]) continue;
321 if isnan(bt[j * nx + i-1]) continue;
322 if isnan(bt[j * nx + i+1]) continue;
323 if isnan(bt_err[j * nx + i]) continue;
324 if isnan(derx_bt[j * nx + i]) continue;
325 if isnan(dery_bt[j * nx + i]) continue;
326 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 */
|
327 mbobra 1.37 err += err_termBt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ))+
328 err_termAt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] )) ;
|
329 xudong 1.27 count_mask++;
330 }
331 }
332
333 *mean_derivative_btotal_ptr = (sum)/(count_mask);
334 *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask);
335 //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
336 //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
337
338 return 0;
|
339 xudong 1.1 }
340
341
342 /*===========================================*/
343 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
344
|
345 mbobra 1.37 int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_termAh, float *err_termBh)
|
346 xudong 1.1 {
|
347 xudong 1.27
348 int nx = dims[0];
349 int ny = dims[1];
350 int i = 0;
351 int j = 0;
352 int count_mask = 0;
353 double sum= 0.0;
354 double err =0.0;
355 *mean_derivative_bh_ptr = 0.0;
356
357 if (nx <= 0 || ny <= 0) return 1;
358
359 /* brute force method of calculating the derivative (no consideration for edges) */
360 for (i = 1; i <= nx-2; i++)
361 {
|
362 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
363 xudong 1.27 {
|
364 mbobra 1.34 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
|
365 mbobra 1.37 err_termAh[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
|
366 xudong 1.27 }
367 }
368
369 /* brute force method of calculating the derivative (no consideration for edges) */
370 for (i = 0; i <= nx-1; i++)
371 {
|
372 mbobra 1.34 for (j = 1; j <= ny-2; j++)
373 {
374 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
|
375 mbobra 1.37 err_termBh[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
|
376 mbobra 1.34 }
|
377 xudong 1.27 }
378
|
379 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
380 ignore the edges for the error terms as those arrays have been initialized to zero.
381 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.*/
|
382 xudong 1.27 i=0;
383 for (j = 0; j <= ny-1; j++)
384 {
385 derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
386 }
387
388 i=nx-1;
389 for (j = 0; j <= ny-1; j++)
390 {
391 derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
392 }
393
394 j=0;
395 for (i = 0; i <= nx-1; i++)
396 {
397 dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
398 }
399
400 j=ny-1;
401 for (i = 0; i <= nx-1; i++)
402 {
403 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;
404 }
405
406
407 for (i = 0; i <= nx-1; i++)
408 {
409 for (j = 0; j <= ny-1; j++)
410 {
411 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
412 if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
413 if isnan(bh[j * nx + i]) continue;
414 if isnan(bh[(j+1) * nx + i]) continue;
415 if isnan(bh[(j-1) * nx + i]) continue;
416 if isnan(bh[j * nx + i-1]) continue;
417 if isnan(bh[j * nx + i+1]) continue;
418 if isnan(bh_err[j * nx + i]) continue;
419 if isnan(derx_bh[j * nx + i]) continue;
420 if isnan(dery_bh[j * nx + i]) continue;
421 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 */
|
422 mbobra 1.37 err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ))+
423 err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] )) ;
|
424 xudong 1.27 count_mask++;
425 }
426 }
427
428 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
429 *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
430 //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
431 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
432
433 return 0;
|
434 xudong 1.1 }
435
436 /*===========================================*/
437 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
438
|
439 mbobra 1.37 int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_termA, float *err_termB)
|
440 xudong 1.1 {
|
441 xudong 1.27
442 int nx = dims[0];
443 int ny = dims[1];
444 int i = 0;
445 int j = 0;
446 int count_mask = 0;
|
447 mbobra 1.34 double sum = 0.0;
|
448 xudong 1.27 double err = 0.0;
|
449 mbobra 1.34 *mean_derivative_bz_ptr = 0.0;
|
450 xudong 1.27
|
451 mbobra 1.34 if (nx <= 0 || ny <= 0) return 1;
|
452 xudong 1.27
453 /* brute force method of calculating the derivative (no consideration for edges) */
454 for (i = 1; i <= nx-2; i++)
455 {
|
456 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
457 xudong 1.27 {
|
458 mbobra 1.34 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
|
459 mbobra 1.37 err_termA[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
|
460 xudong 1.27 }
461 }
462
463 /* brute force method of calculating the derivative (no consideration for edges) */
464 for (i = 0; i <= nx-1; i++)
465 {
|
466 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
467 xudong 1.27 {
|
468 mbobra 1.34 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
|
469 mbobra 1.37 err_termB[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
|
470 xudong 1.27 }
471 }
472
|
473 mbobra 1.34 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
|
474 mbobra 1.35 ignore the edges for the error terms as those arrays have been initialized to zero.
|
475 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.*/
|
476 xudong 1.27 i=0;
477 for (j = 0; j <= ny-1; j++)
478 {
479 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
480 }
481
482 i=nx-1;
483 for (j = 0; j <= ny-1; j++)
484 {
485 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
486 }
487
488 j=0;
489 for (i = 0; i <= nx-1; i++)
490 {
491 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
492 }
493
494 j=ny-1;
495 for (i = 0; i <= nx-1; i++)
496 {
497 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;
498 }
499
500
501 for (i = 0; i <= nx-1; i++)
502 {
503 for (j = 0; j <= ny-1; j++)
504 {
505 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
506 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
507 if isnan(bz[j * nx + i]) continue;
508 if isnan(bz[(j+1) * nx + i]) continue;
509 if isnan(bz[(j-1) * nx + i]) continue;
510 if isnan(bz[j * nx + i-1]) continue;
511 if isnan(bz[j * nx + i+1]) continue;
512 if isnan(bz_err[j * nx + i]) continue;
513 if isnan(derx_bz[j * nx + i]) continue;
514 if isnan(dery_bz[j * nx + i]) continue;
|
515 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 */
|
516 mbobra 1.37 err += err_termB[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) +
517 err_termA[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) ;
|
518 xudong 1.27 count_mask++;
519 }
520 }
521
|
522 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
|
523 xudong 1.27 *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
524 //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
525 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
526
|
527 xudong 1.1 return 0;
528 }
529
530 /*===========================================*/
531 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
532
533 // In discretized space like data pixels,
534 // the current (or curl of B) is calculated as
535 // the integration of the field Bx and By along
536 // the circumference of the data pixel divided by the area of the pixel.
537 // One form of differencing (a word for the differential operator
|
538 xudong 1.27 // in the discretized space) of the curl is expressed as
|
539 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
540 // +dy * (By(i+1,j)+By(i,j)) / 2
541 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
|
542 xudong 1.27 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
543 //
|
544 xudong 1.1 //
545 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
546 // one must perform the following unit conversions:
|
547 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
548 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
|
549 xudong 1.27 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
|
550 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
|
551 mbobra 1.4 //
|
552 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
553 // In that case, we would have the following:
554 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
555 // jz * (35.0)
556 //
557 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
|
558 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)
559 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
|
560 xudong 1.1
561
|
562 xudong 1.27 // Comment out random number generator, which can only run on solar3
|
563 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,
|
564 xudong 1.27 // int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
|
565 mbobra 1.9 // float *noiseby, float *noisebz)
566
567 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
|
568 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)
|
569 mbobra 1.9
|
570 xudong 1.1
|
571 xudong 1.27 {
572 int nx = dims[0];
573 int ny = dims[1];
574 int i = 0;
575 int j = 0;
576 int count_mask = 0;
577
578 if (nx <= 0 || ny <= 0) return 1;
579
580 /* Calculate the derivative*/
581 /* brute force method of calculating the derivative (no consideration for edges) */
582
583 for (i = 1; i <= nx-2; i++)
584 {
|
585 mbobra 1.34 for (j = 0; j <= ny-1; j++)
|
586 xudong 1.27 {
|
587 mbobra 1.34 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
588 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]);
|
589 xudong 1.27 }
590 }
591
592 for (i = 0; i <= nx-1; i++)
593 {
|
594 mbobra 1.34 for (j = 1; j <= ny-2; j++)
|
595 xudong 1.27 {
|
596 mbobra 1.34 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
597 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]);
|
598 xudong 1.27 }
599 }
|
600 mbobra 1.33
|
601 mbobra 1.35 /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
602 ignore the edges for the error terms as those arrays have been initialized to zero.
603 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.*/
604
|
605 xudong 1.27 i=0;
606 for (j = 0; j <= ny-1; j++)
607 {
|
608 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;
|
609 xudong 1.27 }
610
611 i=nx-1;
612 for (j = 0; j <= ny-1; j++)
613 {
|
614 mbobra 1.33 derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
|
615 xudong 1.27 }
616
617 j=0;
618 for (i = 0; i <= nx-1; i++)
619 {
|
620 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;
|
621 xudong 1.27 }
622
623 j=ny-1;
624 for (i = 0; i <= nx-1; i++)
625 {
|
626 mbobra 1.33 dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
|
627 xudong 1.27 }
|
628 mbobra 1.33
|
629 xudong 1.27
|
630 mbobra 1.33 for (i = 0; i <= nx-1; i++)
|
631 xudong 1.27 {
|
632 mbobra 1.33 for (j = 0; j <= ny-1; j++)
|
633 xudong 1.27 {
|
634 mbobra 1.33 // calculate jz at all points
|
635 xudong 1.27 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
|
636 mbobra 1.33 jz_err[j * nx + i] = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;
|
637 xudong 1.27 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
|
638 mbobra 1.33 count_mask++;
|
639 xudong 1.27 }
640 }
641 return 0;
642 }
|
643 mbobra 1.33
|
644 mbobra 1.5 /*===========================================*/
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 mbobra 1.36 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
|
835 xudong 1.1 // 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.32 double sum1=0.0;
|
851 xudong 1.27 double sum2=0.0;
852 double err=0.0;
|
853 mbobra 1.33 *totaljzptr=0.0;
|
854 xudong 1.27
|
855 mbobra 1.33 if (nx <= 0 || ny <= 0) return 1;
|
856 xudong 1.27
|
857 mbobra 1.33 for (i = 0; i < nx; i++)
|
858 xudong 1.27 {
|
859 mbobra 1.33 for (j = 0; j < ny; j++)
|
860 xudong 1.27 {
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.32 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are Amperes per arcsecond */
|
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 mbobra 1.33 return 0;
|
877 xudong 1.1 }
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.31 double sum = 0.0;
|
903 xudong 1.27 double sum1 = 0.0;
904 double err = 0.0;
905 *totpotptr = 0.0;
|
906 mbobra 1.31 *meanpotptr = 0.0;
|
907 xudong 1.27
|
908 mbobra 1.31 if (nx <= 0 || ny <= 0) return 1;
|
909 xudong 1.27
|
910 mbobra 1.31 for (i = 0; i < nx; i++)
|
911 xudong 1.27 {
|
912 mbobra 1.31 for (j = 0; j < ny; j++)
|
913 xudong 1.27 {
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 mbobra 1.31 return 0;
|
940 xudong 1.1 }
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.31 *meanshear_angleptr = 0.0;
|
971 xudong 1.1
|
972 mbobra 1.31 if (nx <= 0 || ny <= 0) return 1;
|
973 xudong 1.27
|
974 mbobra 1.31 for (i = 0; i < nx; i++)
|
975 xudong 1.27 {
|
976 mbobra 1.31 for (j = 0; j < ny; j++)
|
977 xudong 1.27 {
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 //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 xudong 1.27 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 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 xudong 1.27 (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 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 */
1039 //
1040 // Note that there is a restriction on the function fsample()
1041 // If the following occurs:
1042 // nx_out > floor((ny_in-1)/scale + 1)
1043 // ny_out > floor((ny_in-1)/scale + 1),
1044 // where n*_out are the dimensions of the output array and n*_in
1045 // are the dimensions of the input array, fsample() will usually result
1046 // in a segfault (though not always, depending on how the segfault was accessed.)
1047
|
1048 xudong 1.27 int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
1049 float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
|
1050 mbobra 1.30 float *pmap, int nx1, int ny1,
1051 int scale, float *p1pad, int nxp, int nyp, float *pmapn)
1052
|
1053 mbobra 1.29 {
|
1054 xudong 1.27 int nx = dims[0];
1055 int ny = dims[1];
1056 int i = 0;
1057 int j = 0;
|
1058 mbobra 1.30 int index, index1;
|
1059 xudong 1.27 double sum = 0.0;
1060 double err = 0.0;
1061 *Rparam = 0.0;
1062 struct fresize_struct fresboxcar, fresgauss;
|
1063 mbobra 1.30 struct fint_struct fints;
|
1064 xudong 1.27 float sigma = 10.0/2.3548;
1065
|
1066 mbobra 1.30 // set up convolution kernels
|
1067 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1);
1068 init_fresize_gaussian(&fresgauss,sigma,20,1);
|
1069 mbobra 1.29
|
1070 mbobra 1.30 // =============== [STEP 1] ===============
1071 // bin the line-of-sight magnetogram down by a factor of scale
|
1072 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
|
1073 mbobra 1.29
|
1074 mbobra 1.30 // =============== [STEP 2] ===============
1075 // identify positive and negative pixels greater than +/- 150 gauss
1076 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
|
1077 xudong 1.27 for (i = 0; i < nx1; i++)
1078 {
1079 for (j = 0; j < ny1; j++)
1080 {
1081 index = j * nx1 + i;
1082 if (rim[index] > 150)
1083 p1p0[index]=1.0;
1084 else
1085 p1p0[index]=0.0;
1086 if (rim[index] < -150)
1087 p1n0[index]=1.0;
1088 else
1089 p1n0[index]=0.0;
1090 }
1091 }
|
1092 mbobra 1.30
1093 // =============== [STEP 3] ===============
1094 // smooth each of the negative and positive pixel bitmaps
|
1095 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1096 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
|
1097 mbobra 1.30
1098 // =============== [STEP 4] ===============
1099 // find the pixels for which p1p and p1n are both equal to 1.
1100 // this defines the polarity inversion line
|
1101 xudong 1.27 for (i = 0; i < nx1; i++)
1102 {
1103 for (j = 0; j < ny1; j++)
1104 {
1105 index = j * nx1 + i;
|
1106 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
|
1107 xudong 1.27 p1[index]=1.0;
1108 else
1109 p1[index]=0.0;
1110 }
1111 }
|
1112 mbobra 1.30
1113 // pad p1 with zeroes so that the gaussian colvolution in step 5
1114 // does not cut off data within hwidth of the edge
1115
1116 // step i: zero p1pad
1117 for (i = 0; i < nxp; i++)
1118 {
1119 for (j = 0; j < nyp; j++)
1120 {
1121 index = j * nxp + i;
1122 p1pad[index]=0.0;
1123 }
1124 }
1125
1126 // step ii: place p1 at the center of p1pad
1127 for (i = 0; i < nx1; i++)
1128 {
1129 for (j = 0; j < ny1; j++)
1130 {
1131 index = j * nx1 + i;
1132 index1 = (j+20) * nxp + (i+20);
1133 mbobra 1.30 p1pad[index1]=p1[index];
1134 }
1135 }
1136
1137 // =============== [STEP 5] ===============
1138 // convolve the polarity inversion line map with a gaussian
1139 // to identify the region near the plarity inversion line
1140 // the resultant array is called pmap
1141 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);
1142
1143
1144 // select out the nx1 x ny1 non-padded array within pmap
1145 for (i = 0; i < nx1; i++)
1146 {
1147 for (j = 0; j < ny1; j++)
1148 {
1149 index = j * nx1 + i;
1150 index1 = (j+20) * nxp + (i+20);
1151 pmapn[index]=pmap[index1];
1152 }
1153 }
1154 mbobra 1.30
1155 // =============== [STEP 6] ===============
1156 // the R parameter is calculated
|
1157 xudong 1.27 for (i = 0; i < nx1; i++)
1158 {
1159 for (j = 0; j < ny1; j++)
1160 {
1161 index = j * nx1 + i;
|
1162 mbobra 1.32 if isnan(pmapn[index]) continue;
1163 if isnan(rim[index]) continue;
|
1164 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
|
1165 xudong 1.27 }
1166 }
1167
1168 if (sum < 1.0)
1169 *Rparam = 0.0;
1170 else
1171 *Rparam = log10(sum);
|
1172 mbobra 1.30
|
1173 mbobra 1.32 //printf("R_VALUE=%f\n",*Rparam);
1174
|
1175 xudong 1.27 free_fresize(&fresboxcar);
1176 free_fresize(&fresgauss);
|
1177 mbobra 1.30
|
1178 xudong 1.27 return 0;
|
1179 mbobra 1.30
|
1180 xudong 1.1 }
1181
|
1182 mbobra 1.31 /*===========================================*/
1183 /* Example function 16: Lorentz force as defined in Fisher, 2012 */
1184 //
1185 // This calculation is adapted from Xudong's code
1186 // at /proj/cgem/lorentz/apps/lorentz.c
1187
1188 int computeLorentz(float *bx, float *by, float *bz, float *fx, float *fy, float *fz, int *dims,
1189 float *totfx_ptr, float *totfy_ptr, float *totfz_ptr, float *totbsq_ptr,
1190 float *epsx_ptr, float *epsy_ptr, float *epsz_ptr, int *mask, int *bitmask,
1191 float cdelt1, double rsun_ref, double rsun_obs)
1192
1193 {
1194
1195 int nx = dims[0];
1196 int ny = dims[1];
1197 int nxny = nx*ny;
1198 int j = 0;
1199 int index;
1200 double totfx = 0, totfy = 0, totfz = 0;
1201 double bsq = 0, totbsq = 0;
1202 double epsx = 0, epsy = 0, epsz = 0;
1203 mbobra 1.31 double area = cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
1204 double k_h = -1.0 * area / (4. * PI) / 1.0e20;
1205 double k_z = area / (8. * PI) / 1.0e20;
1206
1207 if (nx <= 0 || ny <= 0) return 1;
1208
1209 for (int i = 0; i < nxny; i++)
1210 {
1211 if ( mask[i] < 70 || bitmask[i] < 30 ) continue;
1212 if isnan(bx[i]) continue;
1213 if isnan(by[i]) continue;
1214 if isnan(bz[i]) continue;
|
1215 mbobra 1.32 fx[i] = bx[i] * bz[i] * k_h;
1216 fy[i] = by[i] * bz[i] * k_h;
|
1217 mbobra 1.31 fz[i] = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
1218 bsq = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
1219 totfx += fx[i]; totfy += fy[i]; totfz += fz[i];
1220 totbsq += bsq;
1221 }
1222
1223 *totfx_ptr = totfx;
1224 *totfy_ptr = totfy;
1225 *totfz_ptr = totfz;
1226 *totbsq_ptr = totbsq;
1227 *epsx_ptr = (totfx / k_h) / totbsq;
1228 *epsy_ptr = (totfy / k_h) / totbsq;
1229 *epsz_ptr = (totfz / k_z) / totbsq;
1230
|
1231 mbobra 1.32 //printf("TOTBSQ=%f\n",*totbsq_ptr);
1232
|
1233 mbobra 1.31 return 0;
1234
1235 }
1236
|
1237 xudong 1.1 /*==================KEIJI'S CODE =========================*/
1238
1239 // #include <omp.h>
1240 #include <math.h>
1241
1242 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
1243 {
|
1244 xudong 1.27 /* local workings */
1245 int inx, iny, i, j, n;
1246 /* local array */
1247 float *pfpot, *rdist;
1248 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
1249 rdist=(float *)malloc(sizeof(float) *nnx*nny);
1250 float *bztmp;
1251 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
1252 /* make nan */
1253 // unsigned long long llnan = 0x7ff0000000000000;
1254 // float NAN = (float)(llnan);
1255
1256 // #pragma omp parallel for private (inx)
1257 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
1258 // #pragma omp parallel for private (inx)
1259 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
1260 // #pragma omp parallel for private (inx)
1261 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
1262 // #pragma omp parallel for private (inx)
1263 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
1264 // #pragma omp parallel for private (inx)
1265 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1266 {
1267 float val0 = bz[nnx*iny + inx];
1268 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
1269 }}
1270
1271 // dz is the monopole depth
1272 float dz = 0.001;
1273
1274 // #pragma omp parallel for private (inx)
1275 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
1276 {
1277 float rdd, rdd1, rdd2;
1278 float r;
1279 rdd1 = (float)(inx);
1280 rdd2 = (float)(iny);
1281 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
1282 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
1283 }}
1284
1285 int iwindow;
1286 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
1287 float rwindow;
1288 rwindow = (float)(iwindow);
1289 rwindow = rwindow * rwindow + 0.01; // must be of square
1290
1291 rwindow = 1.0e2; // limit the window size to be 10.
1292
1293 rwindow = sqrt(rwindow);
1294 iwindow = (int)(rwindow);
1295
1296 // #pragma omp parallel for private(inx)
1297 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
|
1298 xudong 1.1 {
|
1299 xudong 1.27 float val0 = bz[nnx*iny + inx];
1300 if (isnan(val0))
1301 {
1302 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
1303 }
1304 else
1305 {
1306 float sum;
1307 sum = 0.0;
1308 int j2, i2;
1309 int j2s, j2e, i2s, i2e;
1310 j2s = iny - iwindow;
1311 j2e = iny + iwindow;
1312 if (j2s < 0){j2s = 0;}
1313 if (j2e > nny){j2e = nny;}
1314 i2s = inx - iwindow;
1315 i2e = inx + iwindow;
1316 if (i2s < 0){i2s = 0;}
|
1317 mbobra 1.37 if (i2e > nnx){i2e = nnx;}
1318
1319 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
1320 {
|
1321 xudong 1.27 float val1 = bztmp[nnx*j2 + i2];
|
1322 mbobra 1.37 float rr, r1, r2;
1323 // r1 = (float)(i2 - inx);
1324 // r2 = (float)(j2 - iny);
1325 // rr = r1*r1 + r2*r2;
1326 // if (rr < rwindow)
1327 // {
|
1328 xudong 1.27 int di, dj;
1329 di = abs(i2 - inx);
1330 dj = abs(j2 - iny);
1331 sum = sum + val1 * rdist[nnx * dj + di] * dz;
|
1332 mbobra 1.37 // }
|
1333 xudong 1.27 } }
1334 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
1335 }
1336 } } // end of OpenMP parallelism
1337
1338 // #pragma omp parallel for private(inx)
1339 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
|
1340 xudong 1.1 {
|
1341 xudong 1.27 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
1342 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
1343 } } // end of OpenMP parallelism
|
1344 mbobra 1.37
|
1345 xudong 1.27 free(rdist);
1346 free(pfpot);
1347 free(bztmp);
|
1348 xudong 1.1 } // end of void func. greenpot
1349
1350
1351 /*===========END OF KEIJI'S CODE =========================*/
|
1352 mbobra 1.14
1353 char *sw_functions_version() // Returns CVS version of sw_functions.c
1354 {
|
1355 xudong 1.27 return strdup("$Id");
|
1356 mbobra 1.14 }
1357
|
1358 xudong 1.1 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps