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
46 #define PI (M_PI)
47 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
48
49 /*===========================================*/
50
51 /* Example function 1: Compute total unsigned flux in units of G/cm^2 */
52
53 // To compute the unsigned flux, we simply calculate
54 // flux = surface integral [(vector Bz) dot (normal vector)],
55 // = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
56 // However, since the field is radial, we will assume cos theta = 1.
57 // Therefore the pixels only need to be corrected for the projection.
58
59 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
60 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
61 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
62 // =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
63 // =(1.30501e15)Gauss*cm^2
64
65 // The disambig mask value selects only the pixels with values of 5 or 7 -- that is,
66 // 5: pixels for which the radial acute disambiguation solution was chosen
67 // 7: pixels for which the radial acute and NRWA disambiguation agree
68
69 int computeAbsFlux(float *bz, int *dims, float *absFlux,
|
70 mbobra 1.3 float *mean_vf_ptr, int *mask, int *bitmask,
|
71 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
72
73 {
74
75 int nx = dims[0], ny = dims[1];
76 int i, j, count_mask=0;
77 double sum=0.0;
78
79 if (nx <= 0 || ny <= 0) return 1;
80
81 *absFlux = 0.0;
82 *mean_vf_ptr =0.0;
83
|
84 mbobra 1.5 for (i = 0; i < nx; i++)
|
85 xudong 1.1 {
|
86 mbobra 1.5 for (j = 0; j < ny; j++)
|
87 xudong 1.1 {
|
88 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
89 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
90 xudong 1.1 sum += (fabs(bz[j * nx + i]));
91 count_mask++;
92 }
93 }
94
95 printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum);
96 printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);
97 *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
98 return 0;
99 }
100
101 /*===========================================*/
102 /* Example function 2: Calculate Bh in units of Gauss */
103 // Native units of Bh are Gauss
104
105 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims,
|
106 mbobra 1.3 float *mean_hf_ptr, int *mask, int *bitmask)
|
107 xudong 1.1
108 {
109
110 int nx = dims[0], ny = dims[1];
111 int i, j, count_mask=0;
112 float sum=0.0;
113 *mean_hf_ptr =0.0;
114
115 if (nx <= 0 || ny <= 0) return 1;
116
|
117 mbobra 1.5 for (i = 0; i < nx; i++)
|
118 xudong 1.1 {
|
119 mbobra 1.5 for (j = 0; j < ny; j++)
|
120 xudong 1.1 {
|
121 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
122 if isnan(by[j * nx + i]) continue;
|
123 xudong 1.1 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
124 sum += bh[j * nx + i];
125 count_mask++;
126 }
127 }
128
129 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
130 printf("*mean_hf_ptr=%f\n",*mean_hf_ptr);
131 return 0;
132 }
133
134 /*===========================================*/
135 /* Example function 3: Calculate Gamma in units of degrees */
136 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
137
138
139 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims,
|
140 mbobra 1.3 float *mean_gamma_ptr, int *mask, int *bitmask)
|
141 xudong 1.1 {
142 int nx = dims[0], ny = dims[1];
143 int i, j, count_mask=0;
144
145 if (nx <= 0 || ny <= 0) return 1;
146
147 *mean_gamma_ptr=0.0;
148 float sum=0.0;
149 int count=0;
150
151 for (i = 0; i < nx; i++)
152 {
153 for (j = 0; j < ny; j++)
154 {
155 if (bh[j * nx + i] > 100)
156 {
|
157 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
158 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
159 xudong 1.1 sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
160 count_mask++;
161 }
162 }
163 }
164
165 *mean_gamma_ptr = sum/count_mask;
166 printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);
167 return 0;
168 }
169
170 /*===========================================*/
171 /* Example function 4: Calculate B_Total*/
172 // Native units of B_Total are in gauss
173
|
174 mbobra 1.3 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask)
|
175 xudong 1.1 {
176
177 int nx = dims[0], ny = dims[1];
178 int i, j, count_mask=0;
179
180 if (nx <= 0 || ny <= 0) return 1;
181
182 for (i = 0; i < nx; i++)
183 {
184 for (j = 0; j < ny; j++)
185 {
|
186 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
187 if isnan(by[j * nx + i]) continue;
188 if isnan(bz[j * nx + i]) continue;
|
189 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]);
190 }
191 }
192 return 0;
193 }
194
195 /*===========================================*/
196 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
197
|
198 mbobra 1.3 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt)
|
199 xudong 1.1 {
200
201 int nx = dims[0], ny = dims[1];
202 int i, j, count_mask=0;
203
204 if (nx <= 0 || ny <= 0) return 1;
205
206 *mean_derivative_btotal_ptr = 0.0;
207 float sum = 0.0;
208
209
210 /* brute force method of calculating the derivative (no consideration for edges) */
211 for (i = 1; i <= nx-2; i++)
212 {
213 for (j = 0; j <= ny-1; j++)
214 {
215 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
216 }
217 }
218
219 /* brute force method of calculating the derivative (no consideration for edges) */
220 xudong 1.1 for (i = 0; i <= nx-1; i++)
221 {
222 for (j = 1; j <= ny-2; j++)
223 {
224 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
225 }
226 }
227
228
229 /* consider the edges */
230 i=0;
231 for (j = 0; j <= ny-1; j++)
232 {
233 derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
234 }
235
236 i=nx-1;
237 for (j = 0; j <= ny-1; j++)
238 {
239 derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
240 }
241 xudong 1.1
242 j=0;
243 for (i = 0; i <= nx-1; i++)
244 {
245 dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
246 }
247
248 j=ny-1;
249 for (i = 0; i <= nx-1; i++)
250 {
251 dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
252 }
253
254
255 for (i = 0; i <= nx-1; i++)
256 {
257 for (j = 0; j <= ny-1; j++)
258 {
|
259 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
260 mbobra 1.5 if isnan(derx_bt[j * nx + i]) continue;
261 if isnan(dery_bt[j * nx + i]) continue;
|
262 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 */
263 count_mask++;
264 }
265 }
266
267 *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
268 printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);
269 return 0;
270 }
271
272
273 /*===========================================*/
274 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
275
|
276 mbobra 1.3 int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)
|
277 xudong 1.1 {
278
279 int nx = dims[0], ny = dims[1];
280 int i, j, count_mask=0;
281
282 if (nx <= 0 || ny <= 0) return 1;
283
284 *mean_derivative_bh_ptr = 0.0;
285 float sum = 0.0;
286
287 /* brute force method of calculating the derivative (no consideration for edges) */
288 for (i = 1; i <= nx-2; i++)
289 {
290 for (j = 0; j <= ny-1; j++)
291 {
292 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
293 }
294 }
295
296 /* brute force method of calculating the derivative (no consideration for edges) */
297 for (i = 0; i <= nx-1; i++)
298 xudong 1.1 {
299 for (j = 1; j <= ny-2; j++)
300 {
301 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
302 }
303 }
304
305
306 /* consider the edges */
307 i=0;
308 for (j = 0; j <= ny-1; j++)
309 {
310 derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
311 }
312
313 i=nx-1;
314 for (j = 0; j <= ny-1; j++)
315 {
316 derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
317 }
318
319 xudong 1.1 j=0;
320 for (i = 0; i <= nx-1; i++)
321 {
322 dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
323 }
324
325 j=ny-1;
326 for (i = 0; i <= nx-1; i++)
327 {
328 dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
329 }
330
331
332 for (i = 0; i <= nx-1; i++)
333 {
334 for (j = 0; j <= ny-1; j++)
335 {
|
336 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
337 mbobra 1.5 if isnan(derx_bh[j * nx + i]) continue;
338 if isnan(dery_bh[j * nx + i]) continue;
|
339 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 */
340 count_mask++;
341 }
342 }
343
344 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
345 return 0;
346 }
347
348 /*===========================================*/
349 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
350
|
351 mbobra 1.3 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)
|
352 xudong 1.1 {
353
354 int nx = dims[0], ny = dims[1];
355 int i, j, count_mask=0;
356
357 if (nx <= 0 || ny <= 0) return 1;
358
359 *mean_derivative_bz_ptr = 0.0;
360 float sum = 0.0;
361
362 /* brute force method of calculating the derivative (no consideration for edges) */
363 for (i = 1; i <= nx-2; i++)
364 {
365 for (j = 0; j <= ny-1; j++)
366 {
|
367 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
368 xudong 1.1 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
369 }
370 }
371
372 /* brute force method of calculating the derivative (no consideration for edges) */
373 for (i = 0; i <= nx-1; i++)
374 {
375 for (j = 1; j <= ny-2; j++)
376 {
|
377 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
378 xudong 1.1 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
379 }
380 }
381
382
383 /* consider the edges */
384 i=0;
385 for (j = 0; j <= ny-1; j++)
386 {
|
387 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
388 xudong 1.1 derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
389 }
390
391 i=nx-1;
392 for (j = 0; j <= ny-1; j++)
393 {
|
394 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
395 xudong 1.1 derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
396 }
397
398 j=0;
399 for (i = 0; i <= nx-1; i++)
400 {
|
401 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
402 xudong 1.1 dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
403 }
404
405 j=ny-1;
406 for (i = 0; i <= nx-1; i++)
407 {
|
408 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
409 xudong 1.1 dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
410 }
411
412
413 for (i = 0; i <= nx-1; i++)
414 {
415 for (j = 0; j <= ny-1; j++)
416 {
417 // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
|
418 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
419 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
420 if isnan(derx_bz[j * nx + i]) continue;
421 if isnan(dery_bz[j * nx + i]) continue;
|
422 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 */
423 count_mask++;
424 }
425 }
426
427 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
428 return 0;
429 }
430
431 /*===========================================*/
432 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */
433
434 // In discretized space like data pixels,
435 // the current (or curl of B) is calculated as
436 // the integration of the field Bx and By along
437 // the circumference of the data pixel divided by the area of the pixel.
438 // One form of differencing (a word for the differential operator
439 // in the discretized space) of the curl is expressed as
440 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2
441 // +dy * (By(i+1,j)+By(i,j)) / 2
442 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
443 xudong 1.1 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
444 //
445 //
446 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
447 // one must perform the following unit conversions:
448 // (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
449 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
450 // where a Tesla is represented as a Newton/Ampere*meter.
|
451 mbobra 1.4 //
|
452 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
453 // In that case, we would have the following:
454 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
455 // jz * (35.0)
456 //
457 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
458 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(1000.)
459 // =(Gauss/pix)(1/CDELT1)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
460 // =(Gauss/pix)(1/0.5)(10^-4)(4*PI*10^7)(722500)(1000.)
461 // =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
462
|
463 mbobra 1.5 int computeJz(float *bx, float *by, int *dims, float *jz,
464 int *mask, int *bitmask,
465 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
|
466 xudong 1.1
467
468
469 {
470
471 int nx = dims[0], ny = dims[1];
472 int i, j, count_mask=0;
473
474 if (nx <= 0 || ny <= 0) return 1;
475 float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
476
477
478 /* brute force method of calculating the derivative (no consideration for edges) */
479 for (i = 1; i <= nx-2; i++)
480 {
481 for (j = 0; j <= ny-1; j++)
482 {
|
483 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
484 xudong 1.1 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
485 }
486 }
487
488 /* brute force method of calculating the derivative (no consideration for edges) */
489 for (i = 0; i <= nx-1; i++)
490 {
491 for (j = 1; j <= ny-2; j++)
492 {
|
493 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
494 xudong 1.1 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
495 }
496 }
497
498
499 /* consider the edges */
500 i=0;
501 for (j = 0; j <= ny-1; j++)
502 {
|
503 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
504 xudong 1.1 derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
505 }
506
507 i=nx-1;
508 for (j = 0; j <= ny-1; j++)
509 {
|
510 mbobra 1.4 if isnan(by[j * nx + i]) continue;
|
511 xudong 1.1 derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
512 }
513
514 j=0;
515 for (i = 0; i <= nx-1; i++)
516 {
|
517 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
518 xudong 1.1 dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
519 }
520
521 j=ny-1;
522 for (i = 0; i <= nx-1; i++)
523 {
|
524 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
|
525 xudong 1.1 dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
526 }
527
528
529 for (i = 0; i <= nx-1; i++)
530 {
531 for (j = 0; j <= ny-1; j++)
532 {
|
533 mbobra 1.5 /* calculate jz at all points */
534 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); /* jz is in units of Gauss/pix */
535 count_mask++;
536 }
537 }
538
539 return 0;
540 }
541
542 /*===========================================*/
543
544 /* Example function 9: Compute quantities on smoothed Jz array */
545
|
546 mbobra 1.6 // All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's
547 // fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels
548 // and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values
549 // of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course,
550 // give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels.
551
|
552 mbobra 1.5 int computeJzsmooth(float *bx, float *by, int *dims, float *jz_smooth,
553 float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
554 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
555
556 {
557
558 int nx = dims[0], ny = dims[1];
559 int i, j, count_mask=0;
560
561 if (nx <= 0 || ny <= 0) return 1;
562
563 *mean_jz_ptr = 0.0;
564 float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
565
566
567 /* 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*/
568 for (i = 0; i <= nx-1; i++)
569 {
570 for (j = 0; j <= ny-1; j++)
571 {
|
572 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
573 mbobra 1.4 if isnan(derx[j * nx + i]) continue;
574 if isnan(dery[j * nx + i]) continue;
|
575 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
576 //printf("%d,%d,%f\n",i,j,jz_smooth[j * nx + i]);
577 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 */
578 us_i += fabs(jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A / m^2 */
|
579 xudong 1.1 count_mask++;
580 }
581 }
582
|
583 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 */
|
584 xudong 1.1 mean_curl = (curl/count_mask);
585 printf("mean_curl=%f\n",mean_curl);
586 printf("cdelt1, what is it?=%f\n",cdelt1);
587 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
588 printf("count_mask=%d\n",count_mask);
|
589 mbobra 1.4 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
|
590 xudong 1.1 return 0;
591
592 }
593
|
594 mbobra 1.5 /*===========================================*/
595
596 /* Example function 10: Twist Parameter, alpha */
|
597 xudong 1.1
|
598 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
599 // for alpha is calculated in the following way (different from Leka and Barnes' approach):
600
601 // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
602 // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
603 // avg alpha = avg Jz / avg Bz
|
604 xudong 1.1
|
605 mbobra 1.6 // The sign is assigned as follows:
606 // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
607 // If this value is > 0, then alpha is > 0.
608 // If this value is < 0, then alpha is <0.
609 //
610 // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
611 // If this value is > 0, then alpha is < 0.
612 // If this value is < 0, then alpha is > 0.
613
|
614 mbobra 1.5 // The units of alpha are in 1/Mm
|
615 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
616 //
617 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
618 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
619 // = 1/Mm
620
|
621 mbobra 1.5 int computeAlpha(float *bz, int *dims, float *jz_smooth, float *mean_alpha_ptr, int *mask, int *bitmask,
|
622 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
|
623 mbobra 1.5
|
624 xudong 1.1 {
625 int nx = dims[0], ny = dims[1];
|
626 mbobra 1.5 int i, j=0;
|
627 xudong 1.1
628 if (nx <= 0 || ny <= 0) return 1;
629
630 *mean_alpha_ptr = 0.0;
|
631 mbobra 1.5 float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=0.0;
|
632 xudong 1.1
633 for (i = 1; i < nx-1; i++)
634 {
635 for (j = 1; j < ny-1; j++)
636 {
|
637 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
638 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
|
639 xudong 1.1 if isnan(bz[j * nx + i]) continue;
640 if (bz[j * nx + i] == 0.0) continue;
|
641 mbobra 1.5 if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i]);
642 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i]);
643 if (bz[j * nx + i] > 0) sum3 += ( jz_smooth[j * nx + i]);
644 if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
645 sum5 += bz[j * nx + i];
|
646 xudong 1.1 }
647 }
|
648 mbobra 1.5
649 sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
|
650 xudong 1.1
|
651 mbobra 1.5 /* Determine the sign of alpha */
652 if ((sum5 > 0) && (sum3 > 0)) sum=sum;
653 if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
654 if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
655 if ((sum5 < 0) && (sum4 > 0)) sum=-sum;
656
657 *mean_alpha_ptr = sum; /* Units are 1/Mm */
|
658 xudong 1.1 return 0;
659 }
660
661 /*===========================================*/
|
662 mbobra 1.5 /* Example function 11: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
|
663 xudong 1.1
664 // The current helicity is defined as Bz*Jz and the units are G^2 / m
665 // The units of Jz are in G/pix; the units of Bz are in G.
666 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/m)
667 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
668 // = G^2 / m.
669
670
|
671 mbobra 1.5 int computeHelicity(float *bz, int *dims, float *jz_smooth, float *mean_ih_ptr, float *total_us_ih_ptr,
|
672 mbobra 1.3 float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
673 xudong 1.1
674 {
675
676 int nx = dims[0], ny = dims[1];
677 int i, j, count_mask=0;
678
679 if (nx <= 0 || ny <= 0) return 1;
680
681 *mean_ih_ptr = 0.0;
682 float sum=0.0, sum2=0.0;
683
|
684 mbobra 1.5 for (i = 0; i < nx; i++)
|
685 xudong 1.1 {
|
686 mbobra 1.5 for (j = 0; j < ny; j++)
|
687 xudong 1.1 {
|
688 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
689 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
|
690 xudong 1.1 if isnan(bz[j * nx + i]) continue;
691 if (bz[j * nx + i] == 0.0) continue;
|
692 mbobra 1.5 if (jz_smooth[j * nx + i] == 0.0) continue;
693 sum += (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
694 sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
|
695 xudong 1.1 count_mask++;
696 }
697 }
698
699 printf("sum/count_mask=%f\n",sum/count_mask);
700 printf("(1/cdelt1)*(rsun_obs/rsun_ref)=%f\n",(1/cdelt1)*(rsun_obs/rsun_ref));
701 *mean_ih_ptr = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
702 *total_us_ih_ptr = sum2; /* Units are G^2 / m */
703 *total_abs_ih_ptr= fabs(sum); /* Units are G^2 / m */
704
705 return 0;
706 }
707
708 /*===========================================*/
|
709 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
|
710 xudong 1.1
711 // The Sum of the Absolute Value per polarity is defined as the following:
712 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
713 // The units of jz are in G/pix. In this case, we would have the following:
714 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
715 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
716
|
717 mbobra 1.5 int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr,
|
718 mbobra 1.3 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
719 xudong 1.1
720 {
721 int nx = dims[0], ny = dims[1];
722 int i, j, count_mask=0;
723
724 if (nx <= 0 || ny <= 0) return 1;
725
726 *totaljzptr=0.0;
727 float sum1=0.0, sum2=0.0;
728
729 for (i = 0; i < nx; i++)
730 {
731 for (j = 0; j < ny; j++)
732 {
|
733 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
734 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
735 mbobra 1.5 if (bz[j * nx + i] > 0) sum1 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
736 if (bz[j * nx + i] <= 0) sum2 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
|
737 xudong 1.1 }
738 }
739
740 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
741 return 0;
742 }
743
744 /*===========================================*/
|
745 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
746 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
747 // automatically yields erg per cubic centimeter for an input B in Gauss.
748 //
749 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
750 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
751 // erg/cm^3(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.)^2
752 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
753 // = erg/cm^3(1.30501e15)
754 // = erg/cm(1/pix^2)
755
756 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims,
|
757 mbobra 1.3 float *meanpotptr, float *totpotptr, int *mask, int *bitmask,
|
758 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
759
760 {
761 int nx = dims[0], ny = dims[1];
762 int i, j, count_mask=0;
763
764 if (nx <= 0 || ny <= 0) return 1;
765
766 *totpotptr=0.0;
767 *meanpotptr=0.0;
768 float sum=0.0;
769
770 for (i = 0; i < nx; i++)
771 {
772 for (j = 0; j < ny; j++)
773 {
|
774 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
775 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
776 if isnan(by[j * nx + i]) continue;
|
777 xudong 1.1 sum += (( ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) - ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i])) )/8.*PI);
778 count_mask++;
779 }
780 }
781
782 *meanpotptr = (sum)/(count_mask); /* Units are ergs per cubic centimeter */
783 *totpotptr = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0)*(count_mask); /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2, units of count_mask are pix^2; therefore, units of totpotptr are ergs per centimeter */
784 return 0;
785 }
786
787 /*===========================================*/
|
788 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 */
|
789 xudong 1.1
790 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
791 float *meanshear_angleptr, float *area_w_shear_gt_45ptr,
792 float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,
|
793 mbobra 1.3 int *mask, int *bitmask)
|
794 xudong 1.1 {
795 int nx = dims[0], ny = dims[1];
796 int i, j;
797
798 if (nx <= 0 || ny <= 0) return 1;
799
800 *area_w_shear_gt_45ptr=0.0;
801 *meanshear_angleptr=0.0;
802 float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0, count_mask=0.0;
803 float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;
804
805 for (i = 0; i < nx; i++)
806 {
807 for (j = 0; j < ny; j++)
808 {
|
809 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
810 xudong 1.1 if isnan(bpx[j * nx + i]) continue;
811 if isnan(bpy[j * nx + i]) continue;
812 if isnan(bpz[j * nx + i]) continue;
813 if isnan(bz[j * nx + i]) continue;
|
814 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
815 if isnan(by[j * nx + i]) continue;
|
816 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
817 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]);
818 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]));
819 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]) );
820 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
821 count ++;
822 sum += shear_angle ;
823 if (shear_angle > 45) count_mask ++;
824 }
825 }
826
827 /* For mean 3D shear angle, area with shear greater than 45*/
828 *meanshear_angleptr = (sum)/(count); /* Units are degrees */
829 printf("count_mask=%f\n",count_mask);
830 printf("count=%f\n",count);
831 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.); /* The area here is a fractional area -- the % of the total area */
832
833 return 0;
834 }
835
836
837 xudong 1.1 /*==================KEIJI'S CODE =========================*/
838
839 // #include <omp.h>
840 #include <math.h>
841
842 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
843 {
844 /* local workings */
845 int inx, iny, i, j, n;
846 /* local array */
847 float *pfpot, *rdist;
848 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
849 rdist=(float *)malloc(sizeof(float) *nnx*nny);
850 float *bztmp;
851 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
852 /* make nan */
853 // unsigned long long llnan = 0x7ff0000000000000;
854 // float NAN = (float)(llnan);
855
856 // #pragma omp parallel for private (inx)
857 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
858 xudong 1.1 // #pragma omp parallel for private (inx)
859 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
860 // #pragma omp parallel for private (inx)
861 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
862 // #pragma omp parallel for private (inx)
863 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
864 // #pragma omp parallel for private (inx)
865 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
866 {
867 float val0 = bz[nnx*iny + inx];
868 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
869 }}
870
871 // dz is the monopole depth
872 float dz = 0.001;
873
874 // #pragma omp parallel for private (inx)
875 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
876 {
877 float rdd, rdd1, rdd2;
878 float r;
879 xudong 1.1 rdd1 = (float)(inx);
880 rdd2 = (float)(iny);
881 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
882 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
883 }}
884
885 int iwindow;
886 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
887 float rwindow;
888 rwindow = (float)(iwindow);
889 rwindow = rwindow * rwindow + 0.01; // must be of square
890
891 rwindow = 1.0e2; // limit the window size to be 10.
892
893 rwindow = sqrt(rwindow);
894 iwindow = (int)(rwindow);
895
896 // #pragma omp parallel for private(inx)
897 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
898 {
899 float val0 = bz[nnx*iny + inx];
900 xudong 1.1 if (isnan(val0))
901 {
902 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
903 }
904 else
905 {
906 float sum;
907 sum = 0.0;
908 int j2, i2;
909 int j2s, j2e, i2s, i2e;
910 j2s = iny - iwindow;
911 j2e = iny + iwindow;
912 if (j2s < 0){j2s = 0;}
913 if (j2e > nny){j2e = nny;}
914 i2s = inx - iwindow;
915 i2e = inx + iwindow;
916 if (i2s < 0){i2s = 0;}
917 if (i2e > nnx){i2e = nnx;}
918
919 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
920 {
921 xudong 1.1 float val1 = bztmp[nnx*j2 + i2];
922 float rr, r1, r2;
923 // r1 = (float)(i2 - inx);
924 // r2 = (float)(j2 - iny);
925 // rr = r1*r1 + r2*r2;
926 // if (rr < rwindow)
927 // {
928 int di, dj;
929 di = abs(i2 - inx);
930 dj = abs(j2 - iny);
931 sum = sum + val1 * rdist[nnx * dj + di] * dz;
932 // }
933 } }
934 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
935 }
936 } } // end of OpenMP parallelism
937
938 // #pragma omp parallel for private(inx)
939 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
940 {
941 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
942 xudong 1.1 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
943 } } // end of OpenMP parallelism
944
945 free(rdist);
946 free(pfpot);
947 free(bztmp);
948 } // end of void func. greenpot
949
950
951 /*===========END OF KEIJI'S CODE =========================*/
952 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps