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