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 int computeJzsmooth(float *bx, float *by, int *dims, float *jz_smooth,
547 float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
548 float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
549
550 {
551
552 int nx = dims[0], ny = dims[1];
553 int i, j, count_mask=0;
554 mbobra 1.5
555 if (nx <= 0 || ny <= 0) return 1;
556
557 *mean_jz_ptr = 0.0;
558 float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
559
560
561 /* 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*/
562 for (i = 0; i <= nx-1; i++)
563 {
564 for (j = 0; j <= ny-1; j++)
565 {
|
566 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
567 mbobra 1.4 if isnan(derx[j * nx + i]) continue;
568 if isnan(dery[j * nx + i]) continue;
|
569 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
570 //printf("%d,%d,%f\n",i,j,jz_smooth[j * nx + i]);
571 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 */
572 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 */
|
573 xudong 1.1 count_mask++;
574 }
575 }
576
|
577 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 */
|
578 xudong 1.1 mean_curl = (curl/count_mask);
579 printf("mean_curl=%f\n",mean_curl);
580 printf("cdelt1, what is it?=%f\n",cdelt1);
581 *mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
582 printf("count_mask=%d\n",count_mask);
|
583 mbobra 1.4 *us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
|
584 xudong 1.1 return 0;
585
586 }
587
|
588 mbobra 1.5 /*===========================================*/
589
590 /* Example function 10: Twist Parameter, alpha */
|
591 xudong 1.1
|
592 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
593 // for alpha is calculated in the following way (different from Leka and Barnes' approach):
594
595 // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
596 // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
597 // avg alpha = avg Jz / avg Bz
|
598 xudong 1.1
|
599 mbobra 1.5 // The units of alpha are in 1/Mm
|
600 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
601 //
602 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
603 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
604 // = 1/Mm
605
|
606 mbobra 1.5 int computeAlpha(float *bz, int *dims, float *jz_smooth, float *mean_alpha_ptr, int *mask, int *bitmask,
|
607 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
|
608 mbobra 1.5
|
609 xudong 1.1 {
610 int nx = dims[0], ny = dims[1];
|
611 mbobra 1.5 int i, j=0;
|
612 xudong 1.1
613 if (nx <= 0 || ny <= 0) return 1;
614
615 *mean_alpha_ptr = 0.0;
|
616 mbobra 1.5 float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=0.0;
|
617 xudong 1.1
618 for (i = 1; i < nx-1; i++)
619 {
620 for (j = 1; j < ny-1; j++)
621 {
|
622 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
623 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
|
624 xudong 1.1 if isnan(bz[j * nx + i]) continue;
625 if (bz[j * nx + i] == 0.0) continue;
|
626 mbobra 1.5 if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i]);
627 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i]);
628 if (bz[j * nx + i] > 0) sum3 += ( jz_smooth[j * nx + i]);
629 if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
630 sum5 += bz[j * nx + i];
|
631 xudong 1.1 }
632 }
|
633 mbobra 1.5
634 sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
|
635 xudong 1.1
|
636 mbobra 1.5 /* Determine the sign of alpha */
637 if ((sum5 > 0) && (sum3 > 0)) sum=sum;
638 if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
639 if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
640 if ((sum5 < 0) && (sum4 > 0)) sum=-sum;
641
642 *mean_alpha_ptr = sum; /* Units are 1/Mm */
|
643 xudong 1.1 return 0;
644 }
645
646 /*===========================================*/
|
647 mbobra 1.5 /* Example function 11: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
|
648 xudong 1.1
649 // The current helicity is defined as Bz*Jz and the units are G^2 / m
650 // The units of Jz are in G/pix; the units of Bz are in G.
651 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/m)
652 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
653 // = G^2 / m.
654
655
|
656 mbobra 1.5 int computeHelicity(float *bz, int *dims, float *jz_smooth, float *mean_ih_ptr, float *total_us_ih_ptr,
|
657 mbobra 1.3 float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
658 xudong 1.1
659 {
660
661 int nx = dims[0], ny = dims[1];
662 int i, j, count_mask=0;
663
664 if (nx <= 0 || ny <= 0) return 1;
665
666 *mean_ih_ptr = 0.0;
667 float sum=0.0, sum2=0.0;
668
|
669 mbobra 1.5 for (i = 0; i < nx; i++)
|
670 xudong 1.1 {
|
671 mbobra 1.5 for (j = 0; j < ny; j++)
|
672 xudong 1.1 {
|
673 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
674 mbobra 1.5 if isnan(jz_smooth[j * nx + i]) continue;
|
675 xudong 1.1 if isnan(bz[j * nx + i]) continue;
676 if (bz[j * nx + i] == 0.0) continue;
|
677 mbobra 1.5 if (jz_smooth[j * nx + i] == 0.0) continue;
678 sum += (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
679 sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
|
680 xudong 1.1 count_mask++;
681 }
682 }
683
684 printf("sum/count_mask=%f\n",sum/count_mask);
685 printf("(1/cdelt1)*(rsun_obs/rsun_ref)=%f\n",(1/cdelt1)*(rsun_obs/rsun_ref));
686 *mean_ih_ptr = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
687 *total_us_ih_ptr = sum2; /* Units are G^2 / m */
688 *total_abs_ih_ptr= fabs(sum); /* Units are G^2 / m */
689
690 return 0;
691 }
692
693 /*===========================================*/
|
694 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
|
695 xudong 1.1
696 // The Sum of the Absolute Value per polarity is defined as the following:
697 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
698 // The units of jz are in G/pix. In this case, we would have the following:
699 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
700 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
701
|
702 mbobra 1.5 int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr,
|
703 mbobra 1.3 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
|
704 xudong 1.1
705 {
706 int nx = dims[0], ny = dims[1];
707 int i, j, count_mask=0;
708
709 if (nx <= 0 || ny <= 0) return 1;
710
711 *totaljzptr=0.0;
712 float sum1=0.0, sum2=0.0;
713
714 for (i = 0; i < nx; i++)
715 {
716 for (j = 0; j < ny; j++)
717 {
|
718 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
719 mbobra 1.4 if isnan(bz[j * nx + i]) continue;
|
720 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);
721 if (bz[j * nx + i] <= 0) sum2 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
|
722 xudong 1.1 }
723 }
724
725 *totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
726 return 0;
727 }
728
729 /*===========================================*/
|
730 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
|
731 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
732 // automatically yields erg per cubic centimeter for an input B in Gauss.
733 //
734 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
735 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
736 // erg/cm^3(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.)^2
737 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
738 // = erg/cm^3(1.30501e15)
739 // = erg/cm(1/pix^2)
740
741 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims,
|
742 mbobra 1.3 float *meanpotptr, float *totpotptr, int *mask, int *bitmask,
|
743 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs)
744
745 {
746 int nx = dims[0], ny = dims[1];
747 int i, j, count_mask=0;
748
749 if (nx <= 0 || ny <= 0) return 1;
750
751 *totpotptr=0.0;
752 *meanpotptr=0.0;
753 float sum=0.0;
754
755 for (i = 0; i < nx; i++)
756 {
757 for (j = 0; j < ny; j++)
758 {
|
759 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
760 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
761 if isnan(by[j * nx + i]) continue;
|
762 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);
763 count_mask++;
764 }
765 }
766
767 *meanpotptr = (sum)/(count_mask); /* Units are ergs per cubic centimeter */
768 *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 */
769 return 0;
770 }
771
772 /*===========================================*/
|
773 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 */
|
774 xudong 1.1
775 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
776 float *meanshear_angleptr, float *area_w_shear_gt_45ptr,
777 float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,
|
778 mbobra 1.3 int *mask, int *bitmask)
|
779 xudong 1.1 {
780 int nx = dims[0], ny = dims[1];
781 int i, j;
782
783 if (nx <= 0 || ny <= 0) return 1;
784
785 *area_w_shear_gt_45ptr=0.0;
786 *meanshear_angleptr=0.0;
787 float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0, count_mask=0.0;
788 float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;
789
790 for (i = 0; i < nx; i++)
791 {
792 for (j = 0; j < ny; j++)
793 {
|
794 mbobra 1.3 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
|
795 xudong 1.1 if isnan(bpx[j * nx + i]) continue;
796 if isnan(bpy[j * nx + i]) continue;
797 if isnan(bpz[j * nx + i]) continue;
798 if isnan(bz[j * nx + i]) continue;
|
799 mbobra 1.4 if isnan(bx[j * nx + i]) continue;
800 if isnan(by[j * nx + i]) continue;
|
801 xudong 1.1 /* For mean 3D shear angle, area with shear greater than 45*/
802 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]);
803 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]));
804 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]) );
805 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
806 count ++;
807 sum += shear_angle ;
808 if (shear_angle > 45) count_mask ++;
809 }
810 }
811
812 /* For mean 3D shear angle, area with shear greater than 45*/
813 *meanshear_angleptr = (sum)/(count); /* Units are degrees */
814 printf("count_mask=%f\n",count_mask);
815 printf("count=%f\n",count);
816 *area_w_shear_gt_45ptr = (count_mask/(count))*(100.); /* The area here is a fractional area -- the % of the total area */
817
818 return 0;
819 }
820
821
822 xudong 1.1 /*==================KEIJI'S CODE =========================*/
823
824 // #include <omp.h>
825 #include <math.h>
826
827 void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
828 {
829 /* local workings */
830 int inx, iny, i, j, n;
831 /* local array */
832 float *pfpot, *rdist;
833 pfpot=(float *)malloc(sizeof(float) *nnx*nny);
834 rdist=(float *)malloc(sizeof(float) *nnx*nny);
835 float *bztmp;
836 bztmp=(float *)malloc(sizeof(float) *nnx*nny);
837 /* make nan */
838 // unsigned long long llnan = 0x7ff0000000000000;
839 // float NAN = (float)(llnan);
840
841 // #pragma omp parallel for private (inx)
842 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
843 xudong 1.1 // #pragma omp parallel for private (inx)
844 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
845 // #pragma omp parallel for private (inx)
846 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
847 // #pragma omp parallel for private (inx)
848 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
849 // #pragma omp parallel for private (inx)
850 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
851 {
852 float val0 = bz[nnx*iny + inx];
853 if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
854 }}
855
856 // dz is the monopole depth
857 float dz = 0.001;
858
859 // #pragma omp parallel for private (inx)
860 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
861 {
862 float rdd, rdd1, rdd2;
863 float r;
864 xudong 1.1 rdd1 = (float)(inx);
865 rdd2 = (float)(iny);
866 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
867 rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
868 }}
869
870 int iwindow;
871 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
872 float rwindow;
873 rwindow = (float)(iwindow);
874 rwindow = rwindow * rwindow + 0.01; // must be of square
875
876 rwindow = 1.0e2; // limit the window size to be 10.
877
878 rwindow = sqrt(rwindow);
879 iwindow = (int)(rwindow);
880
881 // #pragma omp parallel for private(inx)
882 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
883 {
884 float val0 = bz[nnx*iny + inx];
885 xudong 1.1 if (isnan(val0))
886 {
887 pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
888 }
889 else
890 {
891 float sum;
892 sum = 0.0;
893 int j2, i2;
894 int j2s, j2e, i2s, i2e;
895 j2s = iny - iwindow;
896 j2e = iny + iwindow;
897 if (j2s < 0){j2s = 0;}
898 if (j2e > nny){j2e = nny;}
899 i2s = inx - iwindow;
900 i2e = inx + iwindow;
901 if (i2s < 0){i2s = 0;}
902 if (i2e > nnx){i2e = nnx;}
903
904 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
905 {
906 xudong 1.1 float val1 = bztmp[nnx*j2 + i2];
907 float rr, r1, r2;
908 // r1 = (float)(i2 - inx);
909 // r2 = (float)(j2 - iny);
910 // rr = r1*r1 + r2*r2;
911 // if (rr < rwindow)
912 // {
913 int di, dj;
914 di = abs(i2 - inx);
915 dj = abs(j2 - iny);
916 sum = sum + val1 * rdist[nnx * dj + di] * dz;
917 // }
918 } }
919 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
920 }
921 } } // end of OpenMP parallelism
922
923 // #pragma omp parallel for private(inx)
924 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
925 {
926 bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
927 xudong 1.1 by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
928 } } // end of OpenMP parallelism
929
930 free(rdist);
931 free(pfpot);
932 free(bztmp);
933 } // end of void func. greenpot
934
935
936 /*===========END OF KEIJI'S CODE =========================*/
937 /* ---------------- end of this file ----------------*/
|
Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps