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