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