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