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