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