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