JSOC/proj/sharp/apps/sw_functions.c - annotate

Return to sw_functions.c CVS log

Up to [Development] / JSOC / proj / sharp / apps

1 xudong 1.1 /*===========================================
2 xudong 1.27 3 The following 14 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 mbobra 1.29 R_VALUE Karel Schrijver's R parameter
21 xudong 1.27 22 The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and 23 pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD 24 coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data 25 prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI 26 contain nearest-neighbor bitmaps). 27 28 In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig 29 and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values: 30 31 For conf_disambig: 32 50 : not all solutions agree (weak field method applied) 33 60 : not all solutions agree (weak field + annealed) 34 90 : all solutions agree (strong field + annealed) 35 0 : not disambiguated 36 37 For bitmap: 38 1 : weak field outside smooth bounding curve 39 2 : strong field outside smooth bounding curve 40 33 : weak field inside smooth bounding curve 41 34 : strong field inside smooth bounding curve 42 xudong 1.27 43 Written by Monica Bobra 15 August 2012 44 Potential Field code (appended) written by Keiji Hayashi 45 Error analysis modification 21 October 2013 46 47 ===========================================*/
48 xudong 1.1 #include <math.h>
49 mbobra 1.9 #include <mkl.h>
50 xudong 1.1
51 mbobra 1.9 #define PI (M_PI)
52 xudong 1.1 #define MUNAUGHT (0.0000012566370614) /* magnetic constant / 53 54 /===========================================/ 55 56 / Example function 1: Compute total unsigned flux in units of G/cm^2 */ 57
58 xudong 1.27 // To compute the unsigned flux, we simply calculate
59 xudong 1.1 // flux = surface integral [(vector Bz) dot (normal vector)], 60 // = surface integral [(magnitude Bz)(magnitude normal)(cos theta)]. 61 // However, since the field is radial, we will assume cos theta = 1. 62 // Therefore the pixels only need to be corrected for the projection. 63 64 // To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel. 65 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 66 // (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
67 mbobra 1.20 // =Gauss*cm^2
68 xudong 1.1
69 xudong 1.27 int computeAbsFlux(float bz_err, float bz, int dims, float absFlux, 70 float mean_vf_ptr, float mean_vf_err_ptr, float count_mask_ptr, int mask,
71 mbobra 1.9 int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
72 xudong 1.1 73 {
74 xudong 1.27
75 mbobra 1.14 int nx = dims[0]; 76 int ny = dims[1];
77 mbobra 1.15 int i = 0; 78 int j = 0; 79 int count_mask = 0; 80 double sum = 0.0; 81 double err = 0.0;
82 mbobra 1.14 absFlux = 0.0; 83 mean_vf_ptr = 0.0;
84 xudong 1.27 85
86 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
87 xudong 1.27 88 for (i = 0; i < nx; i++)
89 xudong 1.1 {
90 xudong 1.27 for (j = 0; j < ny; j++)
91 xudong 1.1 {
92 xudong 1.27 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 93 if isnan(bz[j * nx + i]) continue; 94 sum += (fabs(bz[j * nx + i])); 95 err += bz_err[j * nx + i]bz_err[j nx + i]; 96 count_mask++; 97 }
98 xudong 1.1 }
99 xudong 1.27 100 mean_vf_ptr = sumcdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 101 mean_vf_err_ptr = (sqrt(err))fabs(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0); // error in the unsigned flux 102 *count_mask_ptr = count_mask; 103 return 0;
104 xudong 1.1 } 105 106 /===========================================/
107 mbobra 1.9 /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
108 xudong 1.1 // Native units of Bh are Gauss 109
110 xudong 1.27 int computeBh(float bx_err, float by_err, float bh_err, float bx, float by, float bz, float bh, int dims,
111 mbobra 1.3 float mean_hf_ptr, int mask, int *bitmask)
112 xudong 1.1 113 {
114 xudong 1.27
115 mbobra 1.14 int nx = dims[0]; 116 int ny = dims[1];
117 mbobra 1.15 int i = 0;
118 xudong 1.27 int j = 0;
119 mbobra 1.15 int count_mask = 0;
120 xudong 1.27 double sum = 0.0;
121 mbobra 1.9 *mean_hf_ptr = 0.0;
122 xudong 1.27
123 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
124 xudong 1.27 125 for (i = 0; i < nx; i++) 126 {
127 mbobra 1.5 for (j = 0; j < ny; j++)
128 xudong 1.27 { 129 if isnan(bx[j * nx + i]) 130 { 131 bh[j * nx + i] = NAN; 132 bh_err[j * nx + i] = NAN; 133 continue; 134 } 135 if isnan(by[j * nx + i]) 136 { 137 bh[j * nx + i] = NAN; 138 bh_err[j * nx + i] = NAN; 139 continue; 140 } 141 bh[j * nx + i] = sqrt( bx[j * nx + i]bx[j nx + i] + by[j * nx + i]by[j nx + i] ); 142 sum += bh[j * nx + i]; 143 bh_err[j * nx + i]=sqrt( bx[j * nx + i]bx[j nx + i]bx_err[j nx + i]bx_err[j nx + i] + by[j * nx + i]by[j nx + i]by_err[j nx + i]by_err[j nx + i])/ bh[j * nx + i]; 144 count_mask++; 145 } 146 } 147
148 xudong 1.1 mean_hf_ptr = sum/(count_mask); // would be divided by nxny if shape of count_mask = shape of magnetogram
149 xudong 1.27
150 xudong 1.1 return 0; 151 } 152 153 /===========================================/ 154 /* Example function 3: Calculate Gamma in units of degrees */ 155 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
156 xudong 1.27 //
157 mbobra 1.20 // Error analysis calculations are done in radians (since derivatives are only true in units of radians), 158 // and multiplied by (180./PI) at the end for consistency in units.
159 xudong 1.1
160 mbobra 1.9 int computeGamma(float bz_err, float bh_err, float bx, float by, float bz, float bh, int dims, 161 float mean_gamma_ptr, float mean_gamma_err_ptr, int mask, int *bitmask)
162 xudong 1.1 {
163 mbobra 1.14 int nx = dims[0]; 164 int ny = dims[1];
165 mbobra 1.15 int i = 0; 166 int j = 0; 167 int count_mask = 0; 168 double sum = 0.0; 169 double err = 0.0; 170 *mean_gamma_ptr = 0.0;
171 xudong 1.27
172 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
173 xudong 1.27 174 for (i = 0; i < nx; i++) 175 { 176 for (j = 0; j < ny; j++) 177 { 178 if (bh[j * nx + i] > 100) 179 { 180 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 181 if isnan(bz[j * nx + i]) continue; 182 if isnan(bz_err[j * nx + i]) continue; 183 if isnan(bh_err[j * nx + i]) continue; 184 if isnan(bh[j * nx + i]) continue; 185 if (bz[j * nx + i] == 0) continue; 186 sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))(180./PI); 187 err += (1/(1+((bh[j nx + i]bh[j nx + i])/(bz[j * nx + i]bz[j nx + i]))))(1/(1+((bh[j nx + i]bh[j nx + i])/(bz[j * nx + i]bz[j nx + i])))) * 188 ( ((bh_err[j * nx + i]bh_err[j nx + i])/(bz[j * nx + i]bz[j nx + i])) + 189 ((bh[j * nx + i]bh[j nx + i]bz_err[j nx + i]bz_err[j nx + i])/(bz[j * nx + i]bz[j nx + i]bz[j nx + i]bz[j nx + i])) ); 190 count_mask++; 191 } 192 } 193 } 194 xudong 1.27 195 mean_gamma_ptr = sum/count_mask; 196 mean_gamma_err_ptr = (sqrt(err)/(count_mask))(180./PI); 197 //printf("MEANGAM=%f\n",mean_gamma_ptr); 198 //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr); 199 return 0;
200 xudong 1.1 } 201 202 /===========================================/ 203 /* Example function 4: Calculate B_Total*/ 204 // Native units of B_Total are in gauss 205
206 mbobra 1.9 int computeB_total(float bx_err, float by_err, float bz_err, float bt_err, float bx, float by, float bz, float bt, int dims, int mask, int *bitmask)
207 xudong 1.1 {
208 xudong 1.27
209 mbobra 1.14 int nx = dims[0]; 210 int ny = dims[1];
211 mbobra 1.15 int i = 0; 212 int j = 0; 213 int count_mask = 0;
214 xudong 1.1 215 if (nx <= 0 \|\| ny <= 0) return 1;
216 xudong 1.27 217 for (i = 0; i < nx; i++) 218 { 219 for (j = 0; j < ny; j++) 220 { 221 if isnan(bx[j * nx + i]) 222 { 223 bt[j * nx + i] = NAN; 224 bt_err[j * nx + i] = NAN; 225 continue; 226 } 227 if isnan(by[j * nx + i]) 228 { 229 bt[j * nx + i] = NAN; 230 bt_err[j * nx + i] = NAN; 231 continue; 232 } 233 if isnan(bz[j * nx + i]) 234 { 235 bt[j * nx + i] = NAN; 236 bt_err[j * nx + i] = NAN; 237 xudong 1.27 continue; 238 } 239 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]); 240 bt_err[j * nx + i]=sqrt(bx[j * nx + i]bx[j nx + i]bx_err[j nx + i]bx_err[j nx + i] + by[j * nx + i]by[j nx + i]by_err[j nx + i]by_err[j nx + i] + bz[j * nx + i]bz[j nx + i]bz_err[j nx + i]bz_err[j nx + i] ) / bt[j * nx + i]; 241 } 242 } 243 return 0;
244 xudong 1.1 } 245 246 /===========================================/ 247 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ 248
249 mbobra 1.9 int computeBtotalderivative(float bt, int dims, float mean_derivative_btotal_ptr, int mask, int bitmask, float derx_bt, float dery_bt, float bt_err, float *mean_derivative_btotal_err_ptr)
250 xudong 1.1 {
251 xudong 1.27
252 mbobra 1.14 int nx = dims[0]; 253 int ny = dims[1];
254 mbobra 1.15 int i = 0; 255 int j = 0; 256 int count_mask = 0;
257 xudong 1.27 double sum = 0.0;
258 mbobra 1.15 double err = 0.0;
259 mbobra 1.14 *mean_derivative_btotal_ptr = 0.0;
260 xudong 1.27
261 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
262 xudong 1.27 263 /* brute force method of calculating the derivative (no consideration for edges) / 264 for (i = 1; i <= nx-2; i++) 265 { 266 for (j = 0; j <= ny-1; j++) 267 { 268 derx_bt[j nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])0.5; 269 } 270 } 271 272 / brute force method of calculating the derivative (no consideration for edges) / 273 for (i = 0; i <= nx-1; i++) 274 { 275 for (j = 1; j <= ny-2; j++) 276 { 277 dery_bt[j nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])0.5; 278 } 279 } 280 281 282 / consider the edges / 283 xudong 1.27 i=0; 284 for (j = 0; j <= ny-1; j++) 285 { 286 derx_bt[j nx + i] = ( (-3bt[j nx + i]) + (4bt[j nx + (i+1)]) - (bt[j * nx + (i+2)]) )0.5; 287 } 288 289 i=nx-1; 290 for (j = 0; j <= ny-1; j++) 291 { 292 derx_bt[j nx + i] = ( (3bt[j nx + i]) + (-4bt[j nx + (i-1)]) - (-bt[j * nx + (i-2)]) )0.5; 293 } 294 295 j=0; 296 for (i = 0; i <= nx-1; i++) 297 { 298 dery_bt[j nx + i] = ( (-3bt[jnx + i]) + (4bt[(j+1) nx + i]) - (bt[(j+2) * nx + i]) )0.5; 299 } 300 301 j=ny-1; 302 for (i = 0; i <= nx-1; i++) 303 { 304 xudong 1.27 dery_bt[j nx + i] = ( (3bt[j nx + i]) + (-4bt[(j-1) nx + i]) - (-bt[(j-2) * nx + i]) )0.5; 305 } 306 307 308 for (i = 1; i <= nx-2; i++) 309 { 310 for (j = 1; j <= ny-2; j++) 311 { 312 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 313 if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue; 314 if isnan(bt[j * nx + i]) continue; 315 if isnan(bt[(j+1) * nx + i]) continue; 316 if isnan(bt[(j-1) * nx + i]) continue; 317 if isnan(bt[j * nx + i-1]) continue; 318 if isnan(bt[j * nx + i+1]) continue; 319 if isnan(bt_err[j * nx + i]) continue; 320 if isnan(derx_bt[j * nx + i]) continue; 321 if isnan(dery_bt[j * nx + i]) continue; 322 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 / 323 err += (((bt[(j+1) nx + i]-bt[(j-1) * nx + i])(bt[(j+1) nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]bt_err[(j+1) nx + i] + bt_err[(j-1) * nx + i]bt_err[(j-1) nx + i])) / (16.0( derx_bt[j nx + i]derx_bt[j nx + i] + dery_bt[j * nx + i]dery_bt[j nx + i] ))+ 324 (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])(bt[j nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]bt_err[j nx + (i+1)] + bt_err[j * nx + (i-1)]bt_err[j nx + (i-1)])) / (16.0( derx_bt[j nx + i]derx_bt[j nx + i] + dery_bt[j * nx + i]dery_bt[j nx + i] )) ; 325 xudong 1.27 count_mask++; 326 } 327 } 328 329 mean_derivative_btotal_ptr = (sum)/(count_mask); 330 mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); 331 //printf("MEANGBT=%f\n",mean_derivative_btotal_ptr); 332 //printf("MEANGBT_err=%f\n",mean_derivative_btotal_err_ptr); 333 334 return 0;
335 xudong 1.1 } 336 337 338 /===========================================/ 339 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ 340
341 mbobra 1.9 int computeBhderivative(float bh, float bh_err, int dims, float mean_derivative_bh_ptr, float mean_derivative_bh_err_ptr, int mask, int bitmask, float derx_bh, float *dery_bh)
342 xudong 1.1 {
343 xudong 1.27 344 int nx = dims[0]; 345 int ny = dims[1]; 346 int i = 0; 347 int j = 0; 348 int count_mask = 0; 349 double sum= 0.0; 350 double err =0.0; 351 mean_derivative_bh_ptr = 0.0; 352 353 if (nx <= 0 \|\| ny <= 0) return 1; 354 355 / brute force method of calculating the derivative (no consideration for edges) / 356 for (i = 1; i <= nx-2; i++) 357 { 358 for (j = 0; j <= ny-1; j++) 359 { 360 derx_bh[j nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])0.5; 361 } 362 } 363 364 xudong 1.27 / brute force method of calculating the derivative (no consideration for edges) / 365 for (i = 0; i <= nx-1; i++) 366 { 367 for (j = 1; j <= ny-2; j++) 368 { 369 dery_bh[j nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])0.5; 370 } 371 } 372 373 374 / consider the edges / 375 i=0; 376 for (j = 0; j <= ny-1; j++) 377 { 378 derx_bh[j nx + i] = ( (-3bh[j nx + i]) + (4bh[j nx + (i+1)]) - (bh[j * nx + (i+2)]) )0.5; 379 } 380 381 i=nx-1; 382 for (j = 0; j <= ny-1; j++) 383 { 384 derx_bh[j nx + i] = ( (3bh[j nx + i]) + (-4bh[j nx + (i-1)]) - (-bh[j * nx + (i-2)]) )0.5; 385 xudong 1.27 } 386 387 j=0; 388 for (i = 0; i <= nx-1; i++) 389 { 390 dery_bh[j nx + i] = ( (-3bh[jnx + i]) + (4bh[(j+1) nx + i]) - (bh[(j+2) * nx + i]) )0.5; 391 } 392 393 j=ny-1; 394 for (i = 0; i <= nx-1; i++) 395 { 396 dery_bh[j nx + i] = ( (3bh[j nx + i]) + (-4bh[(j-1) nx + i]) - (-bh[(j-2) * nx + i]) )0.5; 397 } 398 399 400 for (i = 0; i <= nx-1; i++) 401 { 402 for (j = 0; j <= ny-1; j++) 403 { 404 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 405 if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue; 406 xudong 1.27 if isnan(bh[j * nx + i]) continue; 407 if isnan(bh[(j+1) * nx + i]) continue; 408 if isnan(bh[(j-1) * nx + i]) continue; 409 if isnan(bh[j * nx + i-1]) continue; 410 if isnan(bh[j * nx + i+1]) continue; 411 if isnan(bh_err[j * nx + i]) continue; 412 if isnan(derx_bh[j * nx + i]) continue; 413 if isnan(dery_bh[j * nx + i]) continue; 414 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 / 415 err += (((bh[(j+1) nx + i]-bh[(j-1) * nx + i])(bh[(j+1) nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]bh_err[(j+1) nx + i] + bh_err[(j-1) * nx + i]bh_err[(j-1) nx + i])) / (16.0( derx_bh[j nx + i]derx_bh[j nx + i] + dery_bh[j * nx + i]dery_bh[j nx + i] ))+ 416 (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])(bh[j nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]bh_err[j nx + (i+1)] + bh_err[j * nx + (i-1)]bh_err[j nx + (i-1)])) / (16.0( derx_bh[j nx + i]derx_bh[j nx + i] + dery_bh[j * nx + i]dery_bh[j nx + i] )) ; 417 count_mask++; 418 } 419 } 420 421 mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram 422 mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask) 423 //printf("MEANGBH=%f\n",mean_derivative_bh_ptr); 424 //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr); 425 426 return 0;
427 xudong 1.1 } 428 429 /===========================================/ 430 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ 431
432 mbobra 1.9 int computeBzderivative(float bz, float bz_err, int dims, float mean_derivative_bz_ptr, float mean_derivative_bz_err_ptr, int mask, int bitmask, float derx_bz, float *dery_bz)
433 xudong 1.1 {
434 xudong 1.27 435 int nx = dims[0]; 436 int ny = dims[1]; 437 int i = 0; 438 int j = 0; 439 int count_mask = 0;
440 mbobra 1.15 double sum = 0.0;
441 xudong 1.27 double err = 0.0;
442 mbobra 1.14 *mean_derivative_bz_ptr = 0.0;
443 xudong 1.27
444 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
445 xudong 1.27 446 /* brute force method of calculating the derivative (no consideration for edges) / 447 for (i = 1; i <= nx-2; i++) 448 { 449 for (j = 0; j <= ny-1; j++) 450 { 451 if isnan(bz[j nx + i]) continue; 452 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])0.5; 453 } 454 } 455 456 / brute force method of calculating the derivative (no consideration for edges) / 457 for (i = 0; i <= nx-1; i++) 458 { 459 for (j = 1; j <= ny-2; j++) 460 { 461 if isnan(bz[j nx + i]) continue; 462 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])0.5; 463 } 464 } 465 466 xudong 1.27 467 / consider the edges / 468 i=0; 469 for (j = 0; j <= ny-1; j++) 470 { 471 if isnan(bz[j nx + i]) continue; 472 derx_bz[j * nx + i] = ( (-3bz[j nx + i]) + (4bz[j nx + (i+1)]) - (bz[j * nx + (i+2)]) )0.5; 473 } 474 475 i=nx-1; 476 for (j = 0; j <= ny-1; j++) 477 { 478 if isnan(bz[j nx + i]) continue; 479 derx_bz[j * nx + i] = ( (3bz[j nx + i]) + (-4bz[j nx + (i-1)]) - (-bz[j * nx + (i-2)]) )0.5; 480 } 481 482 j=0; 483 for (i = 0; i <= nx-1; i++) 484 { 485 if isnan(bz[j nx + i]) continue; 486 dery_bz[j * nx + i] = ( (-3bz[jnx + i]) + (4bz[(j+1) nx + i]) - (bz[(j+2) * nx + i]) )0.5; 487 xudong 1.27 } 488 489 j=ny-1; 490 for (i = 0; i <= nx-1; i++) 491 { 492 if isnan(bz[j nx + i]) continue; 493 dery_bz[j * nx + i] = ( (3bz[j nx + i]) + (-4bz[(j-1) nx + i]) - (-bz[(j-2) * nx + i]) )0.5; 494 } 495 496 497 for (i = 0; i <= nx-1; i++) 498 { 499 for (j = 0; j <= ny-1; j++) 500 { 501 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 502 if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue; 503 if isnan(bz[j * nx + i]) continue; 504 if isnan(bz[(j+1) * nx + i]) continue; 505 if isnan(bz[(j-1) * nx + i]) continue; 506 if isnan(bz[j * nx + i-1]) continue; 507 if isnan(bz[j * nx + i+1]) continue; 508 xudong 1.27 if isnan(bz_err[j * nx + i]) continue; 509 if isnan(derx_bz[j * nx + i]) continue; 510 if isnan(dery_bz[j * nx + i]) continue; 511 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 / 512 err += (((bz[(j+1) nx + i]-bz[(j-1) * nx + i])(bz[(j+1) nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]bz_err[(j+1) nx + i] + bz_err[(j-1) * nx + i]bz_err[(j-1) nx + i])) / 513 (16.0( derx_bz[j nx + i]derx_bz[j nx + i] + dery_bz[j * nx + i]dery_bz[j nx + i] )) + 514 (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])(bz[j nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]bz_err[j nx + (i+1)] + bz_err[j * nx + (i-1)]bz_err[j nx + (i-1)])) / 515 (16.0( derx_bz[j nx + i]derx_bz[j nx + i] + dery_bz[j * nx + i]dery_bz[j nx + i] )) ; 516 count_mask++; 517 } 518 } 519
520 xudong 1.1 mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)(ny-2)) if shape of count_mask = shape of magnetogram
521 xudong 1.27 mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask) 522 //printf("MEANGBZ=%f\n",mean_derivative_bz_ptr); 523 //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr); 524
525 xudong 1.1 return 0; 526 } 527 528 /===========================================/ 529 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */ 530 531 // In discretized space like data pixels, 532 // the current (or curl of B) is calculated as 533 // the integration of the field Bx and By along 534 // the circumference of the data pixel divided by the area of the pixel. 535 // One form of differencing (a word for the differential operator
536 xudong 1.27 // in the discretized space) of the curl is expressed as
537 xudong 1.1 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2 538 // +dy * (By(i+1,j)+By(i,j)) / 2 539 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2
540 xudong 1.27 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) 541 //
542 xudong 1.1 // 543 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), 544 // one must perform the following unit conversions:
545 mbobra 1.8 // (Gauss)(1/arcsec)(arcsec/meter)(Newton/GaussAmperemeter)(Ampere^2/Newton)(milliAmpere/Ampere), or 546 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4PI10^-7)( 10^3 milliAmpere/Ampere), or
547 xudong 1.27 // (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
548 xudong 1.1 // where a Tesla is represented as a Newton/Ampere*meter.
549 mbobra 1.4 //
550 xudong 1.1 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 551 // In that case, we would have the following: 552 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4PI10^7)(10^3), or 553 // jz * (35.0) 554 // 555 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
556 mbobra 1.8 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS) 557 // = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
558 xudong 1.1 559
560 xudong 1.27 // Comment out random number generator, which can only run on solar3
561 mbobra 1.9 //int computeJz(float bx_err, float by_err, float bx, float by, int dims, float jz, float jz_err, float jz_err_squared,
562 xudong 1.27 // int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery, float *noisebx,
563 mbobra 1.9 // float noiseby, float noisebz) 564 565 int computeJz(float bx_err, float by_err, float bx, float by, int dims, float jz, float jz_err, float jz_err_squared,
566 xudong 1.27 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery)
567 mbobra 1.9
568 xudong 1.1
569 xudong 1.27 { 570 int nx = dims[0]; 571 int ny = dims[1]; 572 int i = 0; 573 int j = 0; 574 int count_mask = 0; 575 576 if (nx <= 0 \|\| ny <= 0) return 1; 577 578 /* Calculate the derivative/ 579 / brute force method of calculating the derivative (no consideration for edges) / 580 581 for (i = 1; i <= nx-2; i++) 582 { 583 for (j = 0; j <= ny-1; j++) 584 { 585 if isnan(by[j nx + i]) continue; 586 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])0.5; 587 } 588 } 589 590 xudong 1.27 for (i = 0; i <= nx-1; i++) 591 { 592 for (j = 1; j <= ny-2; j++) 593 { 594 if isnan(bx[j nx + i]) continue; 595 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])0.5; 596 } 597 } 598 599 // consider the edges 600 i=0; 601 for (j = 0; j <= ny-1; j++) 602 { 603 if isnan(by[j nx + i]) continue; 604 derx[j * nx + i] = ( (-3by[j nx + i]) + (4by[j nx + (i+1)]) - (by[j * nx + (i+2)]) )0.5; 605 } 606 607 i=nx-1; 608 for (j = 0; j <= ny-1; j++) 609 { 610 if isnan(by[j nx + i]) continue; 611 xudong 1.27 derx[j * nx + i] = ( (3by[j nx + i]) + (-4by[j nx + (i-1)]) - (-by[j * nx + (i-2)]) )0.5; 612 } 613 614 j=0; 615 for (i = 0; i <= nx-1; i++) 616 { 617 if isnan(bx[j nx + i]) continue; 618 dery[j * nx + i] = ( (-3bx[jnx + i]) + (4bx[(j+1) nx + i]) - (bx[(j+2) * nx + i]) )0.5; 619 } 620 621 j=ny-1; 622 for (i = 0; i <= nx-1; i++) 623 { 624 if isnan(bx[j nx + i]) continue; 625 dery[j * nx + i] = ( (3bx[j nx + i]) + (-4bx[(j-1) nx + i]) - (-bx[(j-2) * nx + i]) )0.5; 626 } 627 628 for (i = 1; i <= nx-2; i++) 629 { 630 for (j = 1; j <= ny-2; j++) 631 { 632 xudong 1.27 // calculate jz at all points 633 634 jz[j nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix 635 jz_err[j * nx + i] = 0.5sqrt( (bx_err[(j+1) nx + i]bx_err[(j+1) nx + i]) + (bx_err[(j-1) * nx + i]bx_err[(j-1) nx + i]) + 636 (by_err[j * nx + (i+1)]by_err[j nx + (i+1)]) + (by_err[j * nx + (i-1)]by_err[j nx + (i-1)]) ) ; 637 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]jz_err[j nx + i]); 638 count_mask++; 639 640 } 641 } 642 return 0; 643 }
644 mbobra 1.5 645 /===========================================/ 646
647 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
648 xudong 1.27 // Compute mean and total current on Jz array.
649 mbobra 1.6
650 mbobra 1.9 int computeJzsmooth(float bx, float by, int dims, float jz, float jz_smooth, float jz_err, float jz_rms_err, float jz_err_squared_smooth,
651 xudong 1.27 float mean_jz_ptr, float mean_jz_err_ptr, float us_i_ptr, float us_i_err_ptr, int mask, int bitmask,
652 mbobra 1.9 float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery)
653 mbobra 1.5 654 {
655 xudong 1.27 656 int nx = dims[0]; 657 int ny = dims[1]; 658 int i = 0; 659 int j = 0; 660 int count_mask = 0;
661 mbobra 1.15 double curl = 0.0;
662 xudong 1.27 double us_i = 0.0; 663 double err = 0.0; 664
665 mbobra 1.5 if (nx <= 0 \|\| ny <= 0) return 1;
666 xudong 1.27 667 /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions/ 668 for (i = 0; i <= nx-1; i++) 669 { 670 for (j = 0; j <= ny-1; j++) 671 { 672 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 673 if isnan(derx[j * nx + i]) continue; 674 if isnan(dery[j * nx + i]) continue; 675 if isnan(jz[j * nx + i]) continue; 676 curl += (jz[j * nx + i])(1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.); / curl is in units of mA / m^2 / 677 us_i += fabs(jz[j nx + i])(cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT); /* us_i is in units of A / 678 err += (jz_err[j nx + i]jz_err[j nx + i]); 679 count_mask++; 680 } 681 } 682 683 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above / 684 mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD / 685 mean_jz_err_ptr = (sqrt(err)/count_mask)((1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.)); // error in the quantity MEANJZD 686 687 xudong 1.27 us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ / 688 us_i_err_ptr = (sqrt(err))((cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT)); // error in the quantity TOTUSJZ 689 690 //printf("MEANJZD=%f\n",mean_jz_ptr); 691 //printf("MEANJZD_err=%f\n",mean_jz_err_ptr); 692 693 //printf("TOTUSJZ=%g\n",us_i_ptr); 694 //printf("TOTUSJZ_err=%g\n",us_i_err_ptr); 695
696 xudong 1.1 return 0;
697 xudong 1.27
698 xudong 1.1 } 699
700 mbobra 1.5 /===========================================/ 701 702 /* Example function 10: Twist Parameter, alpha */
703 xudong 1.1
704 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
705 mbobra 1.19 // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
706 xudong 1.27 707 // numerator = sum of all JzBz 708 // denominator = sum of BzBz 709 // alpha = numerator/denominator
710 mbobra 1.6
711 mbobra 1.5 // The units of alpha are in 1/Mm
712 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. 713 //
714 xudong 1.27 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
715 xudong 1.1 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) 716 // = 1/Mm 717
718 mbobra 1.9 int computeAlpha(float jz_err, float bz_err, float bz, int dims, float jz, float jz_smooth, float mean_alpha_ptr, float mean_alpha_err_ptr, int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
719 mbobra 1.5
720 xudong 1.1 {
721 xudong 1.27 int nx = dims[0]; 722 int ny = dims[1]; 723 int i = 0; 724 int j = 0;
725 mbobra 1.19 double alpha_total = 0.0;
726 xudong 1.27 double C = ((1/cdelt1)(rsun_obs/rsun_ref)(1000000.)); 727 double total = 0.0; 728 double A = 0.0; 729 double B = 0.0; 730
731 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
732 xudong 1.27 for (i = 1; i < nx-1; i++) 733 { 734 for (j = 1; j < ny-1; j++) 735 { 736 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 737 if isnan(jz[j * nx + i]) continue; 738 if isnan(bz[j * nx + i]) continue; 739 if (jz[j * nx + i] == 0.0) continue; 740 if (bz[j * nx + i] == 0.0) continue; 741 A += jz[jnx+i]bz[jnx+i]; 742 B += bz[jnx+i]bz[jnx+i]; 743 } 744 } 745 746 for (i = 1; i < nx-1; i++) 747 { 748 for (j = 1; j < ny-1; j++) 749 { 750 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 751 if isnan(jz[j * nx + i]) continue; 752 if isnan(bz[j * nx + i]) continue; 753 xudong 1.27 if (jz[j * nx + i] == 0.0) continue; 754 if (bz[j * nx + i] == 0.0) continue; 755 total += bz[jnx+i]bz[jnx+i]jz_err[jnx+i]jz_err[jnx+i] + (jz[jnx+i]-2bz[jnx+i]A/B)(jz[jnx+i]-2bz[jnx+i]A/B)bz_err[jnx+i]bz_err[jnx+i]; 756 } 757 } 758 759 /* Determine the absolute value of alpha. The units for alpha are 1/Mm / 760 alpha_total = ((A/B)C); 761 mean_alpha_ptr = alpha_total; 762 mean_alpha_err_ptr = (C/B)*(sqrt(total)); 763
764 xudong 1.1 return 0; 765 } 766 767 /===========================================/
768 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
769 xudong 1.1 770 // The current helicity is defined as Bz*Jz and the units are G^2 / m 771 // The units of Jz are in G/pix; the units of Bz are in G.
772 xudong 1.27 // Therefore, the units of BzJz = (Gauss)(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter) 773 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
774 mbobra 1.9 // = G^2 / m.
775 xudong 1.1
776 xudong 1.27 int computeHelicity(float jz_err, float jz_rms_err, float bz_err, float bz, int dims, float jz, float mean_ih_ptr, 777 float mean_ih_err_ptr, float total_us_ih_ptr, float total_abs_ih_ptr,
778 mbobra 1.9 float total_us_ih_err_ptr, float total_abs_ih_err_ptr, int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
779 xudong 1.1 780 {
781 xudong 1.27 782 int nx = dims[0]; 783 int ny = dims[1]; 784 int i = 0; 785 int j = 0; 786 int count_mask = 0;
787 mbobra 1.20 double sum = 0.0; 788 double sum2 = 0.0; 789 double err = 0.0;
790 xudong 1.1 791 if (nx <= 0 \|\| ny <= 0) return 1;
792 xudong 1.27 793 for (i = 0; i < nx; i++)
794 xudong 1.1 {
795 xudong 1.27 for (j = 0; j < ny; j++)
796 xudong 1.1 {
797 xudong 1.27 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 798 if isnan(jz[j * nx + i]) continue; 799 if isnan(bz[j * nx + i]) continue; 800 if isnan(jz_err[j * nx + i]) continue; 801 if isnan(bz_err[j * nx + i]) continue; 802 if (bz[j * nx + i] == 0.0) continue; 803 if (jz[j * nx + i] == 0.0) continue; 804 sum += (jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH 805 sum2 += fabs(jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to TOTUSJH 806 err += (jz_err[j * nx + i]jz_err[j nx + i]bz[j nx + i]bz[j nx + i]) + (bz_err[j * nx + i]bz_err[j nx + i]jz[j nx + i]jz[j nx + i]); 807 count_mask++; 808 } 809 } 810 811 mean_ih_ptr = sum/count_mask ; / Units are G^2 / m ; keyword is MEANJZH */
812 mbobra 1.9 total_us_ih_ptr = sum2 ; / Units are G^2 / m ; keyword is TOTUSJH / 813 total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
814 xudong 1.27 815 mean_ih_err_ptr = (sqrt(err)/count_mask)(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity MEANJZH 816 total_us_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity TOTUSJH 817 total_abs_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity ABSNJZH 818 819 //printf("MEANJZH=%f\n",mean_ih_ptr); 820 //printf("MEANJZH_err=%f\n",mean_ih_err_ptr); 821 822 //printf("TOTUSJH=%f\n",total_us_ih_ptr); 823 //printf("TOTUSJH_err=%f\n",total_us_ih_err_ptr); 824 825 //printf("ABSNJZH=%f\n",total_abs_ih_ptr); 826 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr); 827
828 xudong 1.1 return 0; 829 } 830 831 /===========================================/
832 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
833 xudong 1.1 834 // The Sum of the Absolute Value per polarity is defined as the following: 835 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes. 836 // The units of jz are in G/pix. In this case, we would have the following: 837 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), 838 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
839 mbobra 1.9 // 840 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
841 xudong 1.1
842 xudong 1.27 int computeSumAbsPerPolarity(float jz_err, float bz_err, float bz, float jz, int dims, float totaljzptr, float *totaljz_err_ptr,
843 mbobra 1.3 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
844 xudong 1.1
845 xudong 1.27 { 846 int nx = dims[0]; 847 int ny = dims[1]; 848 int i=0; 849 int j=0; 850 int count_mask=0;
851 mbobra 1.15 double sum1=0.0;
852 xudong 1.27 double sum2=0.0; 853 double err=0.0;
854 mbobra 1.14 *totaljzptr=0.0;
855 xudong 1.27
856 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
857 xudong 1.27 858 for (i = 0; i < nx; i++) 859 { 860 for (j = 0; j < ny; j++) 861 { 862 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 863 if isnan(bz[j * nx + i]) continue; 864 if isnan(jz[j * nx + i]) continue; 865 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 866 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 867 err += (jz_err[j * nx + i]jz_err[j nx + i]); 868 count_mask++; 869 } 870 }
871 xudong 1.1
872 mbobra 1.9 totaljzptr = fabs(sum1) + fabs(sum2); / Units are A */
873 xudong 1.27 totaljz_err_ptr = sqrt(err)(1/cdelt1)fabs((0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs)); 874 //printf("SAVNCPP=%g\n",totaljzptr); 875 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr); 876
877 xudong 1.1 return 0; 878 } 879 880 /===========================================/
881 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
882 xudong 1.1 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
883 xudong 1.27 // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8PI can come out of the integral; thus, 884 // the integral is over B^2 dV and the 8PI is divided at the end.
885 xudong 1.1 // 886 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert 887 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
888 mbobra 1.9 // erg/cm^3(CDELT1^2)(RSUN_REF/RSUN_OBS ^2)(100.^2) 889 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 890 // = erg/cm^3*(1.30501e15)
891 xudong 1.1 // = erg/cm(1/pix^2) 892
893 xudong 1.27 int computeFreeEnergy(float bx_err, float by_err, float bx, float by, float bpx, float bpy, int dims, 894 float meanpotptr, float meanpot_err_ptr, float totpotptr, float totpot_err_ptr, int mask, int *bitmask,
895 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs) 896 897 {
898 xudong 1.27 int nx = dims[0]; 899 int ny = dims[1]; 900 int i = 0; 901 int j = 0; 902 int count_mask = 0;
903 mbobra 1.15 double sum = 0.0;
904 xudong 1.27 double sum1 = 0.0; 905 double err = 0.0; 906 *totpotptr = 0.0;
907 mbobra 1.15 *meanpotptr = 0.0;
908 xudong 1.27
909 mbobra 1.14 if (nx <= 0 \|\| ny <= 0) return 1;
910 xudong 1.27 911 for (i = 0; i < nx; i++) 912 { 913 for (j = 0; j < ny; j++) 914 { 915 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 916 if isnan(bx[j * nx + i]) continue; 917 if isnan(by[j * nx + i]) continue; 918 sum += ( ((bx[j * nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])) )(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0); 919 sum1 += ( ((bx[j * nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])) ); 920 err += 4.0(bx[j nx + i] - bpx[j * nx + i])(bx[j nx + i] - bpx[j * nx + i])(bx_err[j nx + i]bx_err[j nx + i]) + 921 4.0(by[j nx + i] - bpy[j * nx + i])(by[j nx + i] - bpy[j * nx + i])(by_err[j nx + i]by_err[j nx + i]); 922 count_mask++; 923 } 924 } 925 926 /* Units of meanpotptr are ergs per centimeter */
927 mbobra 1.20 meanpotptr = (sum1) / (count_mask8.PI) ; / Units are ergs per cubic centimeter */
928 xudong 1.27 meanpot_err_ptr = (sqrt(err)) / (count_mask8.PI); // error in the quantity (sum)/(count_mask) 929 930 / Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter / 931 totpotptr = (sum)/(8.PI); 932 totpot_err_ptr = (sqrt(err))fabs(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0(1/(8.PI))); 933 934 //printf("MEANPOT=%g\n",meanpotptr); 935 //printf("MEANPOT_err=%g\n",meanpot_err_ptr); 936 937 //printf("TOTPOT=%g\n",totpotptr); 938 //printf("TOTPOT_err=%g\n",totpot_err_ptr); 939
940 xudong 1.1 return 0; 941 } 942 943 /===========================================/
944 mbobra 1.5 /* Example function 14: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
945 xudong 1.1
946 mbobra 1.20 int computeShearAngle(float bx_err, float by_err, float bz_err, float bx, float by, float bz, float bpx, float bpy, float bpz, int dims,
947 mbobra 1.9 float meanshear_angleptr, float meanshear_angle_err_ptr, float area_w_shear_gt_45ptr, int mask, int *bitmask)
948 mbobra 1.21 949
950 xudong 1.27 { 951 int nx = dims[0]; 952 int ny = dims[1]; 953 int i = 0; 954 int j = 0; 955 float count_mask = 0; 956 float count = 0; 957 double dotproduct = 0.0; 958 double magnitude_potential = 0.0; 959 double magnitude_vector = 0.0; 960 double shear_angle = 0.0; 961 double denominator = 0.0; 962 double term1 = 0.0; 963 double term2 = 0.0; 964 double term3 = 0.0; 965 double sumsum = 0.0; 966 double err = 0.0; 967 double part1 = 0.0; 968 double part2 = 0.0; 969 double part3 = 0.0; 970 *area_w_shear_gt_45ptr = 0.0;
971 mbobra 1.15 *meanshear_angleptr = 0.0;
972 xudong 1.1 973 if (nx <= 0 \|\| ny <= 0) return 1;
974 xudong 1.27 975 for (i = 0; i < nx; i++) 976 { 977 for (j = 0; j < ny; j++) 978 { 979 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 980 if isnan(bpx[j * nx + i]) continue; 981 if isnan(bpy[j * nx + i]) continue; 982 if isnan(bpz[j * nx + i]) continue; 983 if isnan(bz[j * nx + i]) continue; 984 if isnan(bx[j * nx + i]) continue; 985 if isnan(by[j * nx + i]) continue; 986 if isnan(bx_err[j * nx + i]) continue; 987 if isnan(by_err[j * nx + i]) continue; 988 if isnan(bz_err[j * nx + i]) continue; 989 990 /* For mean 3D shear angle, percentage with shear greater than 45/ 991 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]); 992 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])); 993 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]) ); 994 //printf("dotproduct=%f\n",dotproduct); 995 xudong 1.27 //printf("magnitude_potential=%f\n",magnitude_potential); 996 //printf("magnitude_vector=%f\n",magnitude_vector); 997 998 shear_angle = acos(dotproduct/(magnitude_potentialmagnitude_vector))(180./PI); 999 sumsum += shear_angle; 1000 //printf("shear_angle=%f\n",shear_angle); 1001 count ++; 1002 1003 if (shear_angle > 45) count_mask ++; 1004 1005 // For the error analysis 1006 1007 term1 = bx[j * nx + i]by[j nx + i]bpy[j nx + i] - by[j * nx + i]by[j nx + i]bpx[j nx + i] + bz[j * nx + i]bx[j nx + i]bpz[j nx + i] - bz[j * nx + i]bz[j nx + i]bpx[j nx + i]; 1008 term2 = bx[j * nx + i]bx[j nx + i]bpy[j nx + i] - bx[j * nx + i]by[j nx + i]bpx[j nx + i] + bx[j * nx + i]bz[j nx + i]bpy[j nx + i] - bz[j * nx + i]by[j nx + i]bpz[j nx + i]; 1009 term3 = bx[j * nx + i]bx[j nx + i]bpz[j nx + i] - bx[j * nx + i]bz[j nx + i]bpx[j nx + i] + by[j * nx + i]by[j nx + i]bpz[j nx + i] - by[j * nx + i]bz[j nx + i]bpy[j nx + i]; 1010 1011 part1 = bx[j * nx + i]bx[j nx + i] + by[j * nx + i]by[j nx + i] + bz[j * nx + i]bz[j nx + i]; 1012 part2 = bpx[j * nx + i]bpx[j nx + i] + bpy[j * nx + i]bpy[j nx + i] + bpz[j * nx + i]bpz[j nx + i]; 1013 part3 = bx[j * nx + i]bpx[j nx + i] + by[j * nx + i]bpy[j nx + i] + bz[j * nx + i]bpz[j nx + i]; 1014 1015 // denominator is squared 1016 xudong 1.27 denominator = part1part1part1part2(1.0-((part3part3)/(part1part2))); 1017 1018 err = (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1019 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1020 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) ; 1021 1022 } 1023 } 1024 /* For mean 3D shear angle, area with shear greater than 45/ 1025 meanshear_angleptr = (sumsum)/(count); /* Units are degrees / 1026 meanshear_angle_err_ptr = (sqrt(err)/count_mask)(180./PI); 1027 1028 / The area here is a fractional area -- the % of the total area. This has no error associated with it. / 1029 area_w_shear_gt_45ptr = (count_mask/(count))(100.0); 1030 1031 //printf("MEANSHR=%f\n",meanshear_angleptr); 1032 //printf("MEANSHR_err=%f\n",meanshear_angle_err_ptr); 1033 //printf("SHRGT45=%f\n",area_w_shear_gt_45ptr); 1034 1035 return 0; 1036 }
1037 xudong 1.1
1038 xudong 1.27 /===========================================/
1039 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 / 1040 // 1041 // Note that there is a restriction on the function fsample() 1042 // If the following occurs: 1043 // nx_out > floor((ny_in-1)/scale + 1) 1044 // ny_out > floor((ny_in-1)/scale + 1), 1045 // where n_out are the dimensions of the output array and n*_in 1046 // are the dimensions of the input array, fsample() will usually result 1047 // in a segfault (though not always, depending on how the segfault was accessed.) 1048
1049 xudong 1.27 int computeR(float bz_err, float los, int dims, float Rparam, float cdelt1, 1050 float rim, float p1p0, float p1n0, float p1p, float p1n, float p1,
1051 mbobra 1.30 float pmap, int nx1, int ny1, 1052 int scale, float p1pad, int nxp, int nyp, float *pmapn) 1053
1054 mbobra 1.29 {
1055 xudong 1.27 int nx = dims[0]; 1056 int ny = dims[1]; 1057 int i = 0; 1058 int j = 0;
1059 mbobra 1.30 int index, index1;
1060 xudong 1.27 double sum = 0.0; 1061 double err = 0.0; 1062 *Rparam = 0.0; 1063 struct fresize_struct fresboxcar, fresgauss;
1064 mbobra 1.30 struct fint_struct fints;
1065 xudong 1.27 float sigma = 10.0/2.3548; 1066
1067 mbobra 1.30 // set up convolution kernels
1068 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1); 1069 init_fresize_gaussian(&fresgauss,sigma,20,1);
1070 mbobra 1.29
1071 mbobra 1.30 // =============== [STEP 1] =============== 1072 // bin the line-of-sight magnetogram down by a factor of scale
1073 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
1074 mbobra 1.29
1075 mbobra 1.30 // =============== [STEP 2] =============== 1076 // identify positive and negative pixels greater than +/- 150 gauss 1077 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
1078 xudong 1.27 for (i = 0; i < nx1; i++) 1079 { 1080 for (j = 0; j < ny1; j++) 1081 { 1082 index = j * nx1 + i; 1083 if (rim[index] > 150) 1084 p1p0[index]=1.0; 1085 else 1086 p1p0[index]=0.0; 1087 if (rim[index] < -150) 1088 p1n0[index]=1.0; 1089 else 1090 p1n0[index]=0.0; 1091 } 1092 }
1093 mbobra 1.30 1094 // =============== [STEP 3] =============== 1095 // smooth each of the negative and positive pixel bitmaps
1096 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 1097 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1098 mbobra 1.30 1099 // =============== [STEP 4] =============== 1100 // find the pixels for which p1p and p1n are both equal to 1. 1101 // this defines the polarity inversion line
1102 xudong 1.27 for (i = 0; i < nx1; i++) 1103 { 1104 for (j = 0; j < ny1; j++) 1105 { 1106 index = j * nx1 + i;
1107 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
1108 xudong 1.27 p1[index]=1.0; 1109 else 1110 p1[index]=0.0; 1111 } 1112 }
1113 mbobra 1.30 1114 // pad p1 with zeroes so that the gaussian colvolution in step 5 1115 // does not cut off data within hwidth of the edge 1116 1117 // step i: zero p1pad 1118 for (i = 0; i < nxp; i++) 1119 { 1120 for (j = 0; j < nyp; j++) 1121 { 1122 index = j * nxp + i; 1123 p1pad[index]=0.0; 1124 } 1125 } 1126 1127 // step ii: place p1 at the center of p1pad 1128 for (i = 0; i < nx1; i++) 1129 { 1130 for (j = 0; j < ny1; j++) 1131 { 1132 index = j * nx1 + i; 1133 index1 = (j+20) * nxp + (i+20); 1134 mbobra 1.30 p1pad[index1]=p1[index]; 1135 } 1136 } 1137 1138 // =============== [STEP 5] =============== 1139 // convolve the polarity inversion line map with a gaussian 1140 // to identify the region near the plarity inversion line 1141 // the resultant array is called pmap 1142 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0); 1143 1144 1145 // select out the nx1 x ny1 non-padded array within pmap 1146 for (i = 0; i < nx1; i++) 1147 { 1148 for (j = 0; j < ny1; j++) 1149 { 1150 index = j * nx1 + i; 1151 index1 = (j+20) * nxp + (i+20); 1152 pmapn[index]=pmap[index1]; 1153 } 1154 } 1155 mbobra 1.30 1156 // =============== [STEP 6] =============== 1157 // the R parameter is calculated
1158 xudong 1.27 for (i = 0; i < nx1; i++) 1159 { 1160 for (j = 0; j < ny1; j++) 1161 { 1162 index = j * nx1 + i;
1163 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
1164 xudong 1.27 } 1165 } 1166 1167 if (sum < 1.0) 1168 Rparam = 0.0; 1169 else 1170 Rparam = log10(sum);
1171 mbobra 1.30
1172 xudong 1.27 free_fresize(&fresboxcar); 1173 free_fresize(&fresgauss);
1174 mbobra 1.30
1175 xudong 1.27 return 0;
1176 mbobra 1.30
1177 xudong 1.1 } 1178 1179 /==================KEIJI'S CODE =========================/ 1180 1181 // #include <omp.h> 1182 #include <math.h> 1183 1184 void greenpot(float bx, float by, float *bz, int nnx, int nny) 1185 {
1186 xudong 1.27 /* local workings / 1187 int inx, iny, i, j, n; 1188 / local array / 1189 float pfpot, rdist; 1190 pfpot=(float )malloc(sizeof(float) nnxnny); 1191 rdist=(float )malloc(sizeof(float) nnxnny); 1192 float bztmp; 1193 bztmp=(float )malloc(sizeof(float) nnxnny); 1194 / make nan / 1195 // unsigned long long llnan = 0x7ff0000000000000; 1196 // float NAN = (float)(llnan); 1197 1198 // #pragma omp parallel for private (inx) 1199 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnxiny+inx] = 0.0;}} 1200 // #pragma omp parallel for private (inx) 1201 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnxiny+inx] = 0.0;}} 1202 // #pragma omp parallel for private (inx) 1203 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnxiny+inx] = 0.0;}} 1204 // #pragma omp parallel for private (inx) 1205 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnxiny+inx] = 0.0;}} 1206 // #pragma omp parallel for private (inx) 1207 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1208 { 1209 float val0 = bz[nnxiny + inx]; 1210 if (isnan(val0)){bztmp[nnxiny + inx] = 0.0;}else{bztmp[nnxiny + inx] = val0;} 1211 }} 1212 1213 // dz is the monopole depth 1214 float dz = 0.001; 1215 1216 // #pragma omp parallel for private (inx) 1217 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1218 { 1219 float rdd, rdd1, rdd2; 1220 float r; 1221 rdd1 = (float)(inx); 1222 rdd2 = (float)(iny); 1223 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz; 1224 rdist[nnxiny+inx] = 1.0/sqrt(rdd); 1225 }} 1226 1227 int iwindow; 1228 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;} 1229 float rwindow; 1230 rwindow = (float)(iwindow); 1231 rwindow = rwindow rwindow + 0.01; // must be of square 1232 1233 rwindow = 1.0e2; // limit the window size to be 10. 1234 1235 rwindow = sqrt(rwindow); 1236 iwindow = (int)(rwindow); 1237 1238 // #pragma omp parallel for private(inx) 1239 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1240 xudong 1.1 {
1241 xudong 1.27 float val0 = bz[nnxiny + inx]; 1242 if (isnan(val0)) 1243 { 1244 pfpot[nnxiny + inx] = 0.0; // hmmm.. NAN; 1245 } 1246 else 1247 { 1248 float sum; 1249 sum = 0.0; 1250 int j2, i2; 1251 int j2s, j2e, i2s, i2e; 1252 j2s = iny - iwindow; 1253 j2e = iny + iwindow; 1254 if (j2s < 0){j2s = 0;} 1255 if (j2e > nny){j2e = nny;} 1256 i2s = inx - iwindow; 1257 i2e = inx + iwindow; 1258 if (i2s < 0){i2s = 0;} 1259 if (i2e > nnx){i2e = nnx;} 1260 1261 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++) 1262 xudong 1.27 { 1263 float val1 = bztmp[nnxj2 + i2]; 1264 float rr, r1, r2; 1265 // r1 = (float)(i2 - inx); 1266 // r2 = (float)(j2 - iny); 1267 // rr = r1r1 + r2r2; 1268 // if (rr < rwindow) 1269 // { 1270 int di, dj; 1271 di = abs(i2 - inx); 1272 dj = abs(j2 - iny); 1273 sum = sum + val1 rdist[nnx * dj + di] * dz; 1274 // } 1275 } } 1276 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition. 1277 } 1278 } } // end of OpenMP parallelism 1279 1280 // #pragma omp parallel for private(inx) 1281 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1282 xudong 1.1 {
1283 xudong 1.27 bx[nnxiny + inx] = -(pfpot[nnxiny + (inx+1)]-pfpot[nnxiny + (inx-1)]) 0.5; 1284 by[nnxiny + inx] = -(pfpot[nnx(iny+1) + inx]-pfpot[nnx(iny-1) + inx]) 0.5; 1285 } } // end of OpenMP parallelism 1286 1287 free(rdist); 1288 free(pfpot); 1289 free(bztmp);
1290 xudong 1.1 } // end of void func. greenpot 1291 1292 1293 /===========END OF KEIJI'S CODE =========================/
1294 mbobra 1.14 1295 char *sw_functions_version() // Returns CVS version of sw_functions.c 1296 {
1297 xudong 1.27 return strdup("$Id");
1298 mbobra 1.14 } 1299
1300 xudong 1.1 /* ---------------- end of this file ----------------*/

Karen Tian