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.33 // 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 mbobra 1.33 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery, float err_term1, float err_term2)
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 mbobra 1.33 derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])0.5; 587 err_term1[j nx + i] = (by_err[j * nx + i+1])(by_err[j nx + i+1]) + (by_err[j * nx + i-1])(by_err[j nx + i-1]);
588 xudong 1.27 } 589 } 590 591 for (i = 0; i <= nx-1; i++) 592 { 593 for (j = 1; j <= ny-2; j++) 594 { 595 if isnan(bx[j * nx + i]) continue;
596 mbobra 1.33 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])0.5; 597 err_term2[j nx + i] = (bx_err[(j+1) * nx + i])(bx_err[(j+1) nx + i]) + (bx_err[(j-1) * nx + i])(bx_err[(j-1) nx + i]);
598 xudong 1.27 } 599 }
600 mbobra 1.33
601 xudong 1.27 // consider the edges 602 i=0; 603 for (j = 0; j <= ny-1; j++) 604 { 605 if isnan(by[j * nx + i]) continue;
606 mbobra 1.33 derx[j * nx + i] = ( (-3by[j nx + i]) + (4by[j nx + (i+1)]) - (by[j * nx + (i+2)]) )0.5; 607 err_term1[j nx + i] = ( (3by_err[j nx + i])(3by_err[j * nx + i]) + (4by_err[j nx + (i+1)])(4by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])(by_err[j nx + (i+2)]) );
608 xudong 1.27 } 609 610 i=nx-1; 611 for (j = 0; j <= ny-1; j++) 612 { 613 if isnan(by[j * nx + i]) continue;
614 mbobra 1.33 derx[j * nx + i] = ( (3by[j nx + i]) + (-4by[j nx + (i-1)]) - (-by[j * nx + (i-2)]) )0.5; 615 err_term1[j nx + i] = ( (3by_err[j nx + i])(3by_err[j * nx + i]) + (4by_err[j nx + (i+1)])(4by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])(by_err[j nx + (i+2)]) );
616 xudong 1.27 } 617 618 j=0; 619 for (i = 0; i <= nx-1; i++) 620 { 621 if isnan(bx[j * nx + i]) continue;
622 mbobra 1.33 dery[j * nx + i] = ( (-3bx[jnx + i]) + (4bx[(j+1) nx + i]) - (bx[(j+2) * nx + i]) )0.5; 623 err_term2[j nx + i] = ( (3bx_err[jnx + i])(3bx_err[jnx + i]) + (4bx_err[(j+1) * nx + i])(4bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])(bx_err[(j+2) nx + i]) );
624 xudong 1.27 } 625 626 j=ny-1; 627 for (i = 0; i <= nx-1; i++) 628 { 629 if isnan(bx[j * nx + i]) continue;
630 mbobra 1.33 dery[j * nx + i] = ( (3bx[j nx + i]) + (-4bx[(j-1) nx + i]) - (-bx[(j-2) * nx + i]) )0.5; 631 err_term2[j nx + i] = ( (3bx_err[jnx + i])(3bx_err[jnx + i]) + (4bx_err[(j+1) * nx + i])(4bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])(bx_err[(j+2) nx + i]) ); 632
633 xudong 1.27 }
634 mbobra 1.33
635 xudong 1.27
636 mbobra 1.33 for (i = 0; i <= nx-1; i++)
637 xudong 1.27 {
638 mbobra 1.33 for (j = 0; j <= ny-1; j++)
639 xudong 1.27 {
640 mbobra 1.33 // calculate jz at all points
641 xudong 1.27 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
642 mbobra 1.33 jz_err[j * nx + i] = 0.5sqrt( err_term1[j nx + i] + err_term2[j * nx + i] ) ;
643 xudong 1.27 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]jz_err[j nx + i]);
644 mbobra 1.33 count_mask++;
645 xudong 1.27 } 646 } 647 return 0; 648 }
649 mbobra 1.33
650 mbobra 1.5 /===========================================/ 651
652 mbobra 1.11 /* Example function 9: Compute quantities on Jz array */
653 xudong 1.27 // Compute mean and total current on Jz array.
654 mbobra 1.6
655 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,
656 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,
657 mbobra 1.9 float cdelt1, double rsun_ref, double rsun_obs,float derx, float dery)
658 mbobra 1.5 659 {
660 xudong 1.27 661 int nx = dims[0]; 662 int ny = dims[1]; 663 int i = 0; 664 int j = 0; 665 int count_mask = 0;
666 mbobra 1.15 double curl = 0.0;
667 xudong 1.27 double us_i = 0.0; 668 double err = 0.0; 669
670 mbobra 1.5 if (nx <= 0 \|\| ny <= 0) return 1;
671 xudong 1.27 672 /* 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/ 673 for (i = 0; i <= nx-1; i++) 674 { 675 for (j = 0; j <= ny-1; j++) 676 { 677 if ( mask[j nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 678 if isnan(derx[j * nx + i]) continue; 679 if isnan(dery[j * nx + i]) continue; 680 if isnan(jz[j * nx + i]) continue; 681 curl += (jz[j * nx + i])(1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.); / curl is in units of mA / m^2 / 682 us_i += fabs(jz[j nx + i])(cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT); /* us_i is in units of A / 683 err += (jz_err[j nx + i]jz_err[j nx + i]); 684 count_mask++; 685 } 686 } 687 688 /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above / 689 mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD / 690 mean_jz_err_ptr = (sqrt(err)/count_mask)((1/cdelt1)(rsun_obs/rsun_ref)(0.00010)(1/MUNAUGHT)(1000.)); // error in the quantity MEANJZD 691 692 xudong 1.27 us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ / 693 us_i_err_ptr = (sqrt(err))((cdelt1/1)(rsun_ref/rsun_obs)(0.00010)(1/MUNAUGHT)); // error in the quantity TOTUSJZ 694 695 //printf("MEANJZD=%f\n",mean_jz_ptr); 696 //printf("MEANJZD_err=%f\n",mean_jz_err_ptr); 697 698 //printf("TOTUSJZ=%g\n",us_i_ptr); 699 //printf("TOTUSJZ_err=%g\n",us_i_err_ptr); 700
701 xudong 1.1 return 0;
702 xudong 1.27
703 xudong 1.1 } 704
705 mbobra 1.5 /===========================================/ 706 707 /* Example function 10: Twist Parameter, alpha */
708 xudong 1.1
709 mbobra 1.5 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
710 mbobra 1.19 // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
711 xudong 1.27 712 // numerator = sum of all JzBz 713 // denominator = sum of BzBz 714 // alpha = numerator/denominator
715 mbobra 1.6
716 mbobra 1.5 // The units of alpha are in 1/Mm
717 xudong 1.1 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. 718 //
719 xudong 1.27 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
720 xudong 1.1 // = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) 721 // = 1/Mm 722
723 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)
724 mbobra 1.5
725 xudong 1.1 {
726 xudong 1.27 int nx = dims[0]; 727 int ny = dims[1]; 728 int i = 0; 729 int j = 0;
730 mbobra 1.19 double alpha_total = 0.0;
731 xudong 1.27 double C = ((1/cdelt1)(rsun_obs/rsun_ref)(1000000.)); 732 double total = 0.0; 733 double A = 0.0; 734 double B = 0.0; 735
736 xudong 1.1 if (nx <= 0 \|\| ny <= 0) return 1;
737 xudong 1.27 for (i = 1; i < nx-1; i++) 738 { 739 for (j = 1; j < ny-1; j++) 740 { 741 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 742 if isnan(jz[j * nx + i]) continue; 743 if isnan(bz[j * nx + i]) continue; 744 if (jz[j * nx + i] == 0.0) continue; 745 if (bz[j * nx + i] == 0.0) continue; 746 A += jz[jnx+i]bz[jnx+i]; 747 B += bz[jnx+i]bz[jnx+i]; 748 } 749 } 750 751 for (i = 1; i < nx-1; i++) 752 { 753 for (j = 1; j < ny-1; j++) 754 { 755 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 756 if isnan(jz[j * nx + i]) continue; 757 if isnan(bz[j * nx + i]) continue; 758 xudong 1.27 if (jz[j * nx + i] == 0.0) continue; 759 if (bz[j * nx + i] == 0.0) continue; 760 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]; 761 } 762 } 763 764 /* Determine the absolute value of alpha. The units for alpha are 1/Mm / 765 alpha_total = ((A/B)C); 766 mean_alpha_ptr = alpha_total; 767 mean_alpha_err_ptr = (C/B)*(sqrt(total)); 768
769 xudong 1.1 return 0; 770 } 771 772 /===========================================/
773 mbobra 1.9 /* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net current helicity) */
774 xudong 1.1 775 // The current helicity is defined as Bz*Jz and the units are G^2 / m 776 // The units of Jz are in G/pix; the units of Bz are in G.
777 xudong 1.27 // Therefore, the units of BzJz = (Gauss)(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter) 778 // = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
779 mbobra 1.9 // = G^2 / m.
780 xudong 1.1
781 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, 782 float mean_ih_err_ptr, float total_us_ih_ptr, float total_abs_ih_ptr,
783 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)
784 xudong 1.1 785 {
786 xudong 1.27 787 int nx = dims[0]; 788 int ny = dims[1]; 789 int i = 0; 790 int j = 0; 791 int count_mask = 0;
792 mbobra 1.20 double sum = 0.0; 793 double sum2 = 0.0; 794 double err = 0.0;
795 xudong 1.1 796 if (nx <= 0 \|\| ny <= 0) return 1;
797 xudong 1.27 798 for (i = 0; i < nx; i++)
799 xudong 1.1 {
800 xudong 1.27 for (j = 0; j < ny; j++)
801 xudong 1.1 {
802 xudong 1.27 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 803 if isnan(jz[j * nx + i]) continue; 804 if isnan(bz[j * nx + i]) continue; 805 if isnan(jz_err[j * nx + i]) continue; 806 if isnan(bz_err[j * nx + i]) continue; 807 if (bz[j * nx + i] == 0.0) continue; 808 if (jz[j * nx + i] == 0.0) continue; 809 sum += (jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH 810 sum2 += fabs(jz[j * nx + i]bz[j nx + i])(1/cdelt1)(rsun_obs/rsun_ref); // contributes to TOTUSJH 811 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]); 812 count_mask++; 813 } 814 } 815 816 mean_ih_ptr = sum/count_mask ; / Units are G^2 / m ; keyword is MEANJZH */
817 mbobra 1.9 total_us_ih_ptr = sum2 ; / Units are G^2 / m ; keyword is TOTUSJH / 818 total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
819 xudong 1.27 820 mean_ih_err_ptr = (sqrt(err)/count_mask)(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity MEANJZH 821 total_us_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity TOTUSJH 822 total_abs_ih_err_ptr = (sqrt(err))(1/cdelt1)(rsun_obs/rsun_ref) ; // error in the quantity ABSNJZH 823 824 //printf("MEANJZH=%f\n",mean_ih_ptr); 825 //printf("MEANJZH_err=%f\n",mean_ih_err_ptr); 826 827 //printf("TOTUSJH=%f\n",total_us_ih_ptr); 828 //printf("TOTUSJH_err=%f\n",total_us_ih_err_ptr); 829 830 //printf("ABSNJZH=%f\n",total_abs_ih_ptr); 831 //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr); 832
833 xudong 1.1 return 0; 834 } 835 836 /===========================================/
837 mbobra 1.5 /* Example function 12: Sum of Absolute Value per polarity */
838 xudong 1.1 839 // The Sum of the Absolute Value per polarity is defined as the following:
840 mbobra 1.32 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per arcsecond.
841 xudong 1.1 // The units of jz are in G/pix. In this case, we would have the following: 842 // Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), 843 // = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
844 mbobra 1.9 // 845 // The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
846 xudong 1.1
847 xudong 1.27 int computeSumAbsPerPolarity(float jz_err, float bz_err, float bz, float jz, int dims, float totaljzptr, float *totaljz_err_ptr,
848 mbobra 1.3 int mask, int bitmask, float cdelt1, double rsun_ref, double rsun_obs)
849 xudong 1.1
850 xudong 1.27 { 851 int nx = dims[0]; 852 int ny = dims[1]; 853 int i=0; 854 int j=0; 855 int count_mask=0;
856 mbobra 1.32 double sum1=0.0;
857 xudong 1.27 double sum2=0.0; 858 double err=0.0;
859 mbobra 1.33 *totaljzptr=0.0;
860 xudong 1.27
861 mbobra 1.33 if (nx <= 0 \|\| ny <= 0) return 1;
862 xudong 1.27
863 mbobra 1.33 for (i = 0; i < nx; i++)
864 xudong 1.27 {
865 mbobra 1.33 for (j = 0; j < ny; j++)
866 xudong 1.27 { 867 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 868 if isnan(bz[j * nx + i]) continue; 869 if isnan(jz[j * nx + i]) continue; 870 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 871 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])(1/cdelt1)(0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs); 872 err += (jz_err[j * nx + i]jz_err[j nx + i]); 873 count_mask++; 874 } 875 }
876 xudong 1.1
877 mbobra 1.32 totaljzptr = fabs(sum1) + fabs(sum2); / Units are Amperes per arcsecond */
878 xudong 1.27 totaljz_err_ptr = sqrt(err)(1/cdelt1)fabs((0.00010)(1/MUNAUGHT)(rsun_ref/rsun_obs)); 879 //printf("SAVNCPP=%g\n",totaljzptr); 880 //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr); 881
882 mbobra 1.33 return 0;
883 xudong 1.1 } 884 885 /===========================================/
886 mbobra 1.5 /* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
887 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
888 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, 889 // the integral is over B^2 dV and the 8PI is divided at the end.
890 xudong 1.1 // 891 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert 892 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
893 mbobra 1.9 // erg/cm^3(CDELT1^2)(RSUN_REF/RSUN_OBS ^2)(100.^2) 894 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 895 // = erg/cm^3*(1.30501e15)
896 xudong 1.1 // = erg/cm(1/pix^2) 897
898 xudong 1.27 int computeFreeEnergy(float bx_err, float by_err, float bx, float by, float bpx, float bpy, int dims, 899 float meanpotptr, float meanpot_err_ptr, float totpotptr, float totpot_err_ptr, int mask, int *bitmask,
900 xudong 1.1 float cdelt1, double rsun_ref, double rsun_obs) 901 902 {
903 xudong 1.27 int nx = dims[0]; 904 int ny = dims[1]; 905 int i = 0; 906 int j = 0; 907 int count_mask = 0;
908 mbobra 1.31 double sum = 0.0;
909 xudong 1.27 double sum1 = 0.0; 910 double err = 0.0; 911 *totpotptr = 0.0;
912 mbobra 1.31 *meanpotptr = 0.0;
913 xudong 1.27
914 mbobra 1.31 if (nx <= 0 \|\| ny <= 0) return 1;
915 xudong 1.27
916 mbobra 1.31 for (i = 0; i < nx; i++)
917 xudong 1.27 {
918 mbobra 1.31 for (j = 0; j < ny; j++)
919 xudong 1.27 { 920 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 921 if isnan(bx[j * nx + i]) continue; 922 if isnan(by[j * nx + i]) continue; 923 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); 924 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])) ); 925 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]) + 926 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]); 927 count_mask++; 928 } 929 } 930 931 /* Units of meanpotptr are ergs per centimeter */
932 mbobra 1.20 meanpotptr = (sum1) / (count_mask8.PI) ; / Units are ergs per cubic centimeter */
933 xudong 1.27 meanpot_err_ptr = (sqrt(err)) / (count_mask8.PI); // error in the quantity (sum)/(count_mask) 934 935 / Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter / 936 totpotptr = (sum)/(8.PI); 937 totpot_err_ptr = (sqrt(err))fabs(cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0(1/(8.PI))); 938 939 //printf("MEANPOT=%g\n",meanpotptr); 940 //printf("MEANPOT_err=%g\n",meanpot_err_ptr); 941 942 //printf("TOTPOT=%g\n",totpotptr); 943 //printf("TOTPOT_err=%g\n",totpot_err_ptr); 944
945 mbobra 1.31 return 0;
946 xudong 1.1 } 947 948 /===========================================/
949 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 */
950 xudong 1.1
951 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,
952 mbobra 1.9 float meanshear_angleptr, float meanshear_angle_err_ptr, float area_w_shear_gt_45ptr, int mask, int *bitmask)
953 mbobra 1.21 954
955 xudong 1.27 { 956 int nx = dims[0]; 957 int ny = dims[1]; 958 int i = 0; 959 int j = 0; 960 float count_mask = 0; 961 float count = 0; 962 double dotproduct = 0.0; 963 double magnitude_potential = 0.0; 964 double magnitude_vector = 0.0; 965 double shear_angle = 0.0; 966 double denominator = 0.0; 967 double term1 = 0.0; 968 double term2 = 0.0; 969 double term3 = 0.0; 970 double sumsum = 0.0; 971 double err = 0.0; 972 double part1 = 0.0; 973 double part2 = 0.0; 974 double part3 = 0.0; 975 *area_w_shear_gt_45ptr = 0.0;
976 mbobra 1.31 *meanshear_angleptr = 0.0;
977 xudong 1.1
978 mbobra 1.31 if (nx <= 0 \|\| ny <= 0) return 1;
979 xudong 1.27
980 mbobra 1.31 for (i = 0; i < nx; i++)
981 xudong 1.27 {
982 mbobra 1.31 for (j = 0; j < ny; j++)
983 xudong 1.27 { 984 if ( mask[j * nx + i] < 70 \|\| bitmask[j * nx + i] < 30 ) continue; 985 if isnan(bpx[j * nx + i]) continue; 986 if isnan(bpy[j * nx + i]) continue; 987 if isnan(bpz[j * nx + i]) continue; 988 if isnan(bz[j * nx + i]) continue; 989 if isnan(bx[j * nx + i]) continue; 990 if isnan(by[j * nx + i]) continue; 991 if isnan(bx_err[j * nx + i]) continue; 992 if isnan(by_err[j * nx + i]) continue; 993 if isnan(bz_err[j * nx + i]) continue; 994 995 /* For mean 3D shear angle, percentage with shear greater than 45/ 996 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]); 997 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])); 998 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]) ); 999 //printf("dotproduct=%f\n",dotproduct); 1000 //printf("magnitude_potential=%f\n",magnitude_potential); 1001 //printf("magnitude_vector=%f\n",magnitude_vector); 1002 1003 shear_angle = acos(dotproduct/(magnitude_potentialmagnitude_vector))(180./PI); 1004 xudong 1.27 sumsum += shear_angle; 1005 //printf("shear_angle=%f\n",shear_angle); 1006 count ++; 1007 1008 if (shear_angle > 45) count_mask ++; 1009 1010 // For the error analysis 1011 1012 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]; 1013 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]; 1014 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]; 1015 1016 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]; 1017 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]; 1018 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]; 1019 1020 // denominator is squared 1021 denominator = part1part1part1part2(1.0-((part3part3)/(part1part2))); 1022 1023 err = (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1024 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) + 1025 xudong 1.27 (term1term1bx_err[j * nx + i]bx_err[j nx + i])/(denominator) ; 1026 1027 } 1028 } 1029 /* For mean 3D shear angle, area with shear greater than 45/ 1030 meanshear_angleptr = (sumsum)/(count); /* Units are degrees / 1031 meanshear_angle_err_ptr = (sqrt(err)/count_mask)(180./PI); 1032 1033 / The area here is a fractional area -- the % of the total area. This has no error associated with it. / 1034 area_w_shear_gt_45ptr = (count_mask/(count))(100.0); 1035 1036 //printf("MEANSHR=%f\n",meanshear_angleptr); 1037 //printf("MEANSHR_err=%f\n",meanshear_angle_err_ptr); 1038 //printf("SHRGT45=%f\n",area_w_shear_gt_45ptr); 1039 1040 return 0; 1041 }
1042 xudong 1.1
1043 xudong 1.27 /===========================================/
1044 mbobra 1.29 /* Example function 15: R parameter as defined in Schrijver, 2007 / 1045 // 1046 // Note that there is a restriction on the function fsample() 1047 // If the following occurs: 1048 // nx_out > floor((ny_in-1)/scale + 1) 1049 // ny_out > floor((ny_in-1)/scale + 1), 1050 // where n_out are the dimensions of the output array and n*_in 1051 // are the dimensions of the input array, fsample() will usually result 1052 // in a segfault (though not always, depending on how the segfault was accessed.) 1053
1054 xudong 1.27 int computeR(float bz_err, float los, int dims, float Rparam, float cdelt1, 1055 float rim, float p1p0, float p1n0, float p1p, float p1n, float p1,
1056 mbobra 1.30 float pmap, int nx1, int ny1, 1057 int scale, float p1pad, int nxp, int nyp, float *pmapn) 1058
1059 mbobra 1.29 {
1060 xudong 1.27 int nx = dims[0]; 1061 int ny = dims[1]; 1062 int i = 0; 1063 int j = 0;
1064 mbobra 1.30 int index, index1;
1065 xudong 1.27 double sum = 0.0; 1066 double err = 0.0; 1067 *Rparam = 0.0; 1068 struct fresize_struct fresboxcar, fresgauss;
1069 mbobra 1.30 struct fint_struct fints;
1070 xudong 1.27 float sigma = 10.0/2.3548; 1071
1072 mbobra 1.30 // set up convolution kernels
1073 xudong 1.27 init_fresize_boxcar(&fresboxcar,1,1); 1074 init_fresize_gaussian(&fresgauss,sigma,20,1);
1075 mbobra 1.29
1076 mbobra 1.30 // =============== [STEP 1] =============== 1077 // bin the line-of-sight magnetogram down by a factor of scale
1078 xudong 1.27 fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
1079 mbobra 1.29
1080 mbobra 1.30 // =============== [STEP 2] =============== 1081 // identify positive and negative pixels greater than +/- 150 gauss 1082 // and label those pixels with a 1.0 in arrays p1p0 and p1n0
1083 xudong 1.27 for (i = 0; i < nx1; i++) 1084 { 1085 for (j = 0; j < ny1; j++) 1086 { 1087 index = j * nx1 + i; 1088 if (rim[index] > 150) 1089 p1p0[index]=1.0; 1090 else 1091 p1p0[index]=0.0; 1092 if (rim[index] < -150) 1093 p1n0[index]=1.0; 1094 else 1095 p1n0[index]=0.0; 1096 } 1097 }
1098 mbobra 1.30 1099 // =============== [STEP 3] =============== 1100 // smooth each of the negative and positive pixel bitmaps
1101 xudong 1.27 fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0); 1102 fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
1103 mbobra 1.30 1104 // =============== [STEP 4] =============== 1105 // find the pixels for which p1p and p1n are both equal to 1. 1106 // this defines the polarity inversion line
1107 xudong 1.27 for (i = 0; i < nx1; i++) 1108 { 1109 for (j = 0; j < ny1; j++) 1110 { 1111 index = j * nx1 + i;
1112 mbobra 1.30 if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
1113 xudong 1.27 p1[index]=1.0; 1114 else 1115 p1[index]=0.0; 1116 } 1117 }
1118 mbobra 1.30 1119 // pad p1 with zeroes so that the gaussian colvolution in step 5 1120 // does not cut off data within hwidth of the edge 1121 1122 // step i: zero p1pad 1123 for (i = 0; i < nxp; i++) 1124 { 1125 for (j = 0; j < nyp; j++) 1126 { 1127 index = j * nxp + i; 1128 p1pad[index]=0.0; 1129 } 1130 } 1131 1132 // step ii: place p1 at the center of p1pad 1133 for (i = 0; i < nx1; i++) 1134 { 1135 for (j = 0; j < ny1; j++) 1136 { 1137 index = j * nx1 + i; 1138 index1 = (j+20) * nxp + (i+20); 1139 mbobra 1.30 p1pad[index1]=p1[index]; 1140 } 1141 } 1142 1143 // =============== [STEP 5] =============== 1144 // convolve the polarity inversion line map with a gaussian 1145 // to identify the region near the plarity inversion line 1146 // the resultant array is called pmap 1147 fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0); 1148 1149 1150 // select out the nx1 x ny1 non-padded array within pmap 1151 for (i = 0; i < nx1; i++) 1152 { 1153 for (j = 0; j < ny1; j++) 1154 { 1155 index = j * nx1 + i; 1156 index1 = (j+20) * nxp + (i+20); 1157 pmapn[index]=pmap[index1]; 1158 } 1159 } 1160 mbobra 1.30 1161 // =============== [STEP 6] =============== 1162 // the R parameter is calculated
1163 xudong 1.27 for (i = 0; i < nx1; i++) 1164 { 1165 for (j = 0; j < ny1; j++) 1166 { 1167 index = j * nx1 + i;
1168 mbobra 1.32 if isnan(pmapn[index]) continue; 1169 if isnan(rim[index]) continue;
1170 mbobra 1.30 sum += pmapn[index]*abs(rim[index]);
1171 xudong 1.27 } 1172 } 1173 1174 if (sum < 1.0) 1175 Rparam = 0.0; 1176 else 1177 Rparam = log10(sum);
1178 mbobra 1.30
1179 mbobra 1.32 //printf("R_VALUE=%f\n",*Rparam); 1180
1181 xudong 1.27 free_fresize(&fresboxcar); 1182 free_fresize(&fresgauss);
1183 mbobra 1.30
1184 xudong 1.27 return 0;
1185 mbobra 1.30
1186 xudong 1.1 } 1187
1188 mbobra 1.31 /===========================================/ 1189 /* Example function 16: Lorentz force as defined in Fisher, 2012 / 1190 // 1191 // This calculation is adapted from Xudong's code 1192 // at /proj/cgem/lorentz/apps/lorentz.c 1193 1194 int computeLorentz(float bx, float by, float bz, float fx, float fy, float fz, int dims, 1195 float totfx_ptr, float totfy_ptr, float totfz_ptr, float totbsq_ptr, 1196 float epsx_ptr, float epsy_ptr, float epsz_ptr, int mask, int bitmask, 1197 float cdelt1, double rsun_ref, double rsun_obs) 1198 1199 { 1200 1201 int nx = dims[0]; 1202 int ny = dims[1]; 1203 int nxny = nxny; 1204 int j = 0; 1205 int index; 1206 double totfx = 0, totfy = 0, totfz = 0; 1207 double bsq = 0, totbsq = 0; 1208 double epsx = 0, epsy = 0, epsz = 0; 1209 mbobra 1.31 double area = cdelt1cdelt1(rsun_ref/rsun_obs)(rsun_ref/rsun_obs)100.0100.0; 1210 double k_h = -1.0 area / (4. * PI) / 1.0e20; 1211 double k_z = area / (8. * PI) / 1.0e20; 1212 1213 if (nx <= 0 \|\| ny <= 0) return 1; 1214 1215 for (int i = 0; i < nxny; i++) 1216 { 1217 if ( mask[i] < 70 \|\| bitmask[i] < 30 ) continue; 1218 if isnan(bx[i]) continue; 1219 if isnan(by[i]) continue; 1220 if isnan(bz[i]) continue;
1221 mbobra 1.32 fx[i] = bx[i] * bz[i] * k_h; 1222 fy[i] = by[i] * bz[i] * k_h;
1223 mbobra 1.31 fz[i] = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z; 1224 bsq = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i]; 1225 totfx += fx[i]; totfy += fy[i]; totfz += fz[i]; 1226 totbsq += bsq; 1227 } 1228 1229 totfx_ptr = totfx; 1230 totfy_ptr = totfy; 1231 totfz_ptr = totfz; 1232 totbsq_ptr = totbsq; 1233 epsx_ptr = (totfx / k_h) / totbsq; 1234 epsy_ptr = (totfy / k_h) / totbsq; 1235 *epsz_ptr = (totfz / k_z) / totbsq; 1236
1237 mbobra 1.32 //printf("TOTBSQ=%f\n",*totbsq_ptr); 1238
1239 mbobra 1.31 return 0; 1240 1241 } 1242
1243 xudong 1.1 /==================KEIJI'S CODE =========================/ 1244 1245 // #include <omp.h> 1246 #include <math.h> 1247 1248 void greenpot(float bx, float by, float *bz, int nnx, int nny) 1249 {
1250 xudong 1.27 /* local workings / 1251 int inx, iny, i, j, n; 1252 / local array / 1253 float pfpot, rdist; 1254 pfpot=(float )malloc(sizeof(float) nnxnny); 1255 rdist=(float )malloc(sizeof(float) nnxnny); 1256 float bztmp; 1257 bztmp=(float )malloc(sizeof(float) nnxnny); 1258 / make nan / 1259 // unsigned long long llnan = 0x7ff0000000000000; 1260 // float NAN = (float)(llnan); 1261 1262 // #pragma omp parallel for private (inx) 1263 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnxiny+inx] = 0.0;}} 1264 // #pragma omp parallel for private (inx) 1265 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnxiny+inx] = 0.0;}} 1266 // #pragma omp parallel for private (inx) 1267 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnxiny+inx] = 0.0;}} 1268 // #pragma omp parallel for private (inx) 1269 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnxiny+inx] = 0.0;}} 1270 // #pragma omp parallel for private (inx) 1271 xudong 1.27 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1272 { 1273 float val0 = bz[nnxiny + inx]; 1274 if (isnan(val0)){bztmp[nnxiny + inx] = 0.0;}else{bztmp[nnxiny + inx] = val0;} 1275 }} 1276 1277 // dz is the monopole depth 1278 float dz = 0.001; 1279 1280 // #pragma omp parallel for private (inx) 1281 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 1282 { 1283 float rdd, rdd1, rdd2; 1284 float r; 1285 rdd1 = (float)(inx); 1286 rdd2 = (float)(iny); 1287 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz; 1288 rdist[nnxiny+inx] = 1.0/sqrt(rdd); 1289 }} 1290 1291 int iwindow; 1292 xudong 1.27 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;} 1293 float rwindow; 1294 rwindow = (float)(iwindow); 1295 rwindow = rwindow rwindow + 0.01; // must be of square 1296 1297 rwindow = 1.0e2; // limit the window size to be 10. 1298 1299 rwindow = sqrt(rwindow); 1300 iwindow = (int)(rwindow); 1301 1302 // #pragma omp parallel for private(inx) 1303 for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
1304 xudong 1.1 {
1305 xudong 1.27 float val0 = bz[nnxiny + inx]; 1306 if (isnan(val0)) 1307 { 1308 pfpot[nnxiny + inx] = 0.0; // hmmm.. NAN; 1309 } 1310 else 1311 { 1312 float sum; 1313 sum = 0.0; 1314 int j2, i2; 1315 int j2s, j2e, i2s, i2e; 1316 j2s = iny - iwindow; 1317 j2e = iny + iwindow; 1318 if (j2s < 0){j2s = 0;} 1319 if (j2e > nny){j2e = nny;} 1320 i2s = inx - iwindow; 1321 i2e = inx + iwindow; 1322 if (i2s < 0){i2s = 0;} 1323 if (i2e > nnx){i2e = nnx;} 1324 1325 for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++) 1326 xudong 1.27 { 1327 float val1 = bztmp[nnxj2 + i2]; 1328 float rr, r1, r2; 1329 // r1 = (float)(i2 - inx); 1330 // r2 = (float)(j2 - iny); 1331 // rr = r1r1 + r2r2; 1332 // if (rr < rwindow) 1333 // { 1334 int di, dj; 1335 di = abs(i2 - inx); 1336 dj = abs(j2 - iny); 1337 sum = sum + val1 rdist[nnx * dj + di] * dz; 1338 // } 1339 } } 1340 pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition. 1341 } 1342 } } // end of OpenMP parallelism 1343 1344 // #pragma omp parallel for private(inx) 1345 for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
1346 xudong 1.1 {
1347 xudong 1.27 bx[nnxiny + inx] = -(pfpot[nnxiny + (inx+1)]-pfpot[nnxiny + (inx-1)]) 0.5; 1348 by[nnxiny + inx] = -(pfpot[nnx(iny+1) + inx]-pfpot[nnx(iny-1) + inx]) 0.5; 1349 } } // end of OpenMP parallelism 1350 1351 free(rdist); 1352 free(pfpot); 1353 free(bztmp);
1354 xudong 1.1 } // end of void func. greenpot 1355 1356 1357 /===========END OF KEIJI'S CODE =========================/
1358 mbobra 1.14 1359 char *sw_functions_version() // Returns CVS version of sw_functions.c 1360 {
1361 xudong 1.27 return strdup("$Id");
1362 mbobra 1.14 } 1363
1364 xudong 1.1 /* ---------------- end of this file ----------------*/

Karen Tian