Return to sharp_functions.c CVS log Up to [Development] / JSOC / proj / sharp / apps

 ``` 1 mbobra 1.1 /* sharp_functions.c */ 2 #define PI (3.141592653589793) 3 4 /*===========================================*/ 5 6 /* Example function 1: Compute total unsigned flux in units of G/cm^2, Mean Vertical Field */ 7 8 // To convert G to G/cm^2, simply multiply by the number of square centimeters per pixel. 9 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 10 // Gauss(CDELT1)(RSUN_REF/RSUN_OBS)(100.) 11 // = Gauss(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 12 // = Gauss(1.30501e15) 13 14 int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask) 15 { 16 17 int nx = dims[0], ny = dims[1]; 18 int i, j, count_mask=0; 19 double sum=0.0; 20 21 if (nx <= 0 || ny <= 0) return 1; 22 mbobra 1.1 23 *absFlux = 0.0; 24 *mean_vf_ptr =0.0; 25 26 for (j = 0; j < ny; j++) 27 { 28 for (i = 0; i < nx; i++) 29 { 30 if (mask[j * nx + i] <= 1) continue; 31 sum += (fabs(bz[j * nx + i])); 32 count_mask++; 33 } 34 } 35 36 printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum); 37 *mean_vf_ptr = sum*(1.30501e15); 38 return 0; 39 } 40 41 /*===========================================*/ 42 /* Example function 2: Calculate Bh in units of Gauss */ 43 mbobra 1.1 // Native units of Bh are Gauss 44 45 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask) 46 { 47 48 int nx = dims[0], ny = dims[1]; 49 int i, j, count_mask=0; 50 float sum=0.0; 51 *mean_hf_ptr =0.0; 52 53 if (nx <= 0 || ny <= 0) return 1; 54 55 for (j = 0; j < ny; j++) 56 { 57 for (i = 0; i < nx; i++) 58 { 59 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] ); 60 sum += bh[j * nx + i]; 61 count_mask++; 62 } 63 } 64 mbobra 1.1 65 *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram 66 printf("*mean_hf_ptr=%f\n",*mean_hf_ptr); 67 return 0; 68 } 69 70 /*===========================================*/ 71 72 /* Example function 3: Calculate Gamma in units of degrees 73 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI) 74 75 this is wrong, also the divide by nx*ny is wrong */ 76 77 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_gamma_ptr, int *mask) 78 { 79 int nx = dims[0], ny = dims[1]; 80 int i, j, count_mask=0; 81 82 if (nx <= 0 || ny <= 0) return 1; 83 84 *mean_gamma_ptr=0.0; 85 mbobra 1.1 float sum=0.0; 86 int count=0; 87 88 for (i = 0; i < nx; i++) 89 { 90 for (j = 0; j < ny; j++) 91 { 92 if (bh[j * nx + i] > 100) 93 { 94 if (mask[j * nx + i] <= 1) continue; 95 sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI)); 96 count_mask++; 97 } 98 } 99 } 100 101 *mean_gamma_ptr = sum/count_mask; 102 printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr); 103 return 0; 104 } 105 106 mbobra 1.1 /*===========================================*/ 107 108 /* Example function 4: Calculate B_Total*/ 109 // Native units of B_Total are in gauss 110 111 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask) 112 { 113 114 int nx = dims[0], ny = dims[1]; 115 int i, j, count_mask=0; 116 117 if (nx <= 0 || ny <= 0) return 1; 118 119 for (i = 0; i < nx; i++) 120 { 121 for (j = 0; j < ny; j++) 122 { 123 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]); 124 } 125 } 126 return 0; 127 mbobra 1.1 } 128 129 /*===========================================*/ 130 131 /* Example function 5: Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ 132 133 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask) 134 { 135 136 int nx = dims[0], ny = dims[1]; 137 int i, j, count_mask=0; 138 139 if (nx <= 0 || ny <= 0) return 1; 140 141 *mean_derivative_btotal_ptr = 0.0; 142 float derx, dery, sum = 0.0; 143 144 for (i = 1; i < nx-1; i++) 145 { 146 for (j = 1; j < ny-1; j++) 147 { 148 mbobra 1.1 if (mask[j * nx + i] <= 1) continue; 149 /* Missing a (*0.5) */ 150 derx = bt[j * nx + i+1] - bt[j * nx + i-1]; 151 dery = bt[(j+1) * nx + i] - bt[(j-1) * nx + i]; 152 sum += sqrt(derx*derx + dery*dery); 153 count_mask++; 154 } 155 } 156 157 *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram 158 printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr); 159 return 0; 160 } 161 162 /*===========================================*/ 163 /* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ 164 165 int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask) 166 { 167 168 int nx = dims[0], ny = dims[1]; 169 mbobra 1.1 int i, j, count_mask=0; 170 171 if (nx <= 0 || ny <= 0) return 1; 172 173 *mean_derivative_bh_ptr = 0.0; 174 float derx, dery, sum = 0.0; 175 176 for (i = 1; i < nx-1; i++) 177 { 178 for (j = 1; j < ny-1; j++) 179 { 180 if (mask[j * nx + i] <= 1) continue; 181 /* Missing a (*0.5) */ 182 derx = bh[j * nx + i+1] - bh[j * nx + i-1]; 183 //printf("derx=%f\n",derx); 184 dery = bh[(j+1) * nx + i] - bh[(j-1) * nx + i]; 185 // printf("dery=%f\n",dery); 186 sum += sqrt(derx*derx + dery*dery); 187 //printf("sum=%f\n",sum); 188 count_mask++; 189 //printf("count_mask=%d\n",count_mask); 190 mbobra 1.1 } 191 } 192 193 *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram 194 printf("*mean_derivative_bh_ptr=%f,nx=%d,ny=%d,sum=%f\n",*mean_derivative_bh_ptr,nx,ny,sum); 195 return 0; 196 } 197 198 /*===========================================*/ 199 /* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ 200 201 int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask) 202 { 203 204 int nx = dims[0], ny = dims[1]; 205 int i, j, count_mask=0; 206 207 if (nx <= 0 || ny <= 0) return 1; 208 209 *mean_derivative_bz_ptr = 0.0; 210 float derx, dery, sum = 0.0; 211 mbobra 1.1 212 for (i = 1; i < nx-1; i++) 213 { 214 for (j = 1; j < ny-1; j++) 215 { 216 if (mask[j * nx + i] <= 1) continue; 217 /* Missing a (*0.5) */ 218 derx = bz[j * nx + i+1] - bz[j * nx + i-1]; 219 dery = bz[(j+1) * nx + i] - bz[(j-1) * nx + i]; 220 sum += sqrt(derx*derx + dery*dery); 221 count_mask++; 222 } 223 } 224 225 *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram 226 227 return 0; 228 } 229 230 /*===========================================*/ 231 /* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */ 232 mbobra 1.1 233 // In discretized space like data pixels, 234 // the current (or curl of B) is calculated as 235 // the integration of the field Bx and By along 236 // the circumference of the data pixel divided by the area of the pixel. 237 // One form of differencing (a word for the differential operator 238 // in the discretized space) of the curl is expressed as 239 // (dx * (Bx(i,j-1)+Bx(i,j)) / 2 240 // +dy * (By(i+1,j)+By(i,j)) / 2 241 // -dx * (Bx(i,j+1)+Bx(i,j)) / 2 242 // -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) 243 // 244 // 245 // To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), 246 // one must perform the following unit conversions: 247 // (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or 248 // (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), 249 // where a Tesla is represented as a Newton/Ampere*meter. 250 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 251 // In that case, we would have the following: 252 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or 253 mbobra 1.1 // jz * (35.0) 254 // 255 // The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following: 256 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or 257 // (Gauss/pix)(1.81*10^10) 258 // (Gauss/pix)(18100000000.) 259 260 int computeJz(float *bx, float *by, int *dims, float *jz, float *mean_jz_ptr, float *us_i_ptr, int *mask) 261 { 262 263 int nx = dims[0], ny = dims[1]; 264 int i, j, count_mask=0; 265 266 if (nx <= 0 || ny <= 0) return 1; 267 268 *mean_jz_ptr = 0.0; 269 float derx, dery, curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0; 270 271 for (i = 1; i < nx-1; i++) 272 { 273 for (j = 1; j < ny-1; j++) 274 mbobra 1.1 { 275 if (mask[j * nx + i] <= 1) continue; 276 derx = (by[j * nx + i+1] - by[j * nx + i-1])*0.5; 277 dery = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5; 278 curl += (derx-dery)*(35.0); /* curl is in units of mA/m^2 */ 279 jz[j * nx + i] = (derx-dery); /* jz is in units of Gauss per pixel */ 280 us_i += fabs(derx-dery)*(18100000000.) ; /* us_i is in units of A */ 281 count_mask++; 282 } 283 } 284 285 mean_curl = (curl/count_mask); 286 printf("mean_curl=%f\n",mean_curl); 287 *mean_jz_ptr = curl/(count_mask); 288 printf("count_mask=%d\n",count_mask); 289 *us_i_ptr = (us_i); 290 return 0; 291 } 292 293 294 295 mbobra 1.1 /*===========================================*/ 296 /* Example function 9: Twist Parameter, alpha */ 297 298 // The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm 299 // The units of Jz are in G/pix; the units of Bz are in G. 300 // Therefore, the units of Jz/Bz = (pix/arcsec)(arcsec/meter)(meter/Mm), or 301 // Jz/Bz = (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6). 302 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 303 // In that case, we would have the following: 304 // (Gauss/pix)(1/0.5)(1/722500)(10^6), or 305 // (Gauss/pix)*2.7 306 307 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask) 308 { 309 int nx = dims[0], ny = dims[1]; 310 int i, j, count_mask=0; 311 312 if (nx <= 0 || ny <= 0) return 1; 313 314 *mean_alpha_ptr = 0.0; 315 float aa, bb, cc, bznew, alpha2, sum=0.0; 316 mbobra 1.1 317 for (i = 1; i < nx-1; i++) 318 { 319 for (j = 1; j < ny-1; j++) 320 { 321 if (mask[j * nx + i] <= 1) continue; 322 if isnan(jz[j * nx + i]) continue; 323 if isnan(bz[j * nx + i]) continue; 324 if (bz[j * nx + i] == 0.0) continue; 325 sum += (jz[j * nx + i] / bz[j * nx + i])*2.7 ; /* the units for (jz) Gauss/pix; the units for bz are Gauss */ 326 //printf("sum=%f\n",sum); 327 //printf("jz[j * nx + i]=%f\n",jz[j * nx + i]); 328 //printf("bz[j * nx + i]=%f\n",bz[j * nx + i]); 329 //printf("(jz[j * nx + i] / bz[j * nx + i])*2.7=%f\n",(jz[j * nx + i] / bz[j * nx + i])*2.7); 330 count_mask++; 331 //printf("count_mask=%d\n",count_mask); 332 } 333 } 334 335 printf("count_mask=%d\n",count_mask); 336 printf("sum=%f\n",sum); 337 mbobra 1.1 *mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */ 338 return 0; 339 } 340 341 /*===========================================*/ 342 343 /* Example function 10: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */ 344 345 // The current helicity is defined as Bz*Jz and the units are G^2 / m 346 // The units of Jz are in G/pix; the units of Bz are in G. 347 // Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF) = G^2 / m. 348 // As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). 349 // In that case, we would have the following: 350 // (Gauss/pix)(1/0.5)(1/722500), or 351 // (Gauss/pix)*0.000002768166 352 353 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr, int *mask) 354 { 355 356 357 int nx = dims[0], ny = dims[1]; 358 mbobra 1.1 int i, j, count_mask=0; 359 360 if (nx <= 0 || ny <= 0) return 1; 361 362 *mean_ih_ptr = 0.0; 363 float sum=0.0, sum2=0.0; 364 365 for (j = 0; j < ny; j++) 366 { 367 for (i = 0; i < nx; i++) 368 { 369 if (mask[j * nx + i] <= 1) continue; 370 if isnan(jz[j * nx + i]) continue; 371 if isnan(bz[j * nx + i]) continue; 372 if (bz[j * nx + i] == 0.0) continue; 373 if (jz[j * nx + i] == 0.0) continue; 374 sum += (jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166); 375 sum2 += fabs(jz[j * nx + i]*(bz[j * nx + i]))*(0.000002768166); 376 count_mask++; 377 printf("(jz[j * nx + i])=%f\n",(jz[j * nx + i])); 378 printf("(bz[j * nx + i])=%f\n",(bz[j * nx + i])); 379 mbobra 1.1 printf("(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166)=%f\n",(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166)); 380 printf("sum=%f\n",sum); 381 printf("sum2=%f\n",sum2); 382 printf("count_mask=%d\n",count_mask); 383 } 384 } 385 386 printf("sum/count_mask=%f\n",sum/count_mask); 387 *mean_ih_ptr = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */ 388 *total_us_ih_ptr = sum2; /* Units are G^2 / m */ 389 *total_abs_ih_ptr= fabs(sum); /* Units are G^2 / m */ 390 391 return 0; 392 } 393 394 /*===========================================*/ 395 396 /* Example function 11: Sum of Absolute Value per polarity */ 397 398 // The Sum of the Absolute Value per polarity is defined as the following: 399 // fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes. 400 mbobra 1.1 // The units of jz are in G/pix. In this case, we would have the following: 401 // (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or 402 // (Gauss/pix)(1.81*10^10) 403 // (Gauss/pix)(18100000000.) 404 405 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int *mask) 406 { 407 408 409 int nx = dims[0], ny = dims[1]; 410 int i, j, count_mask=0; 411 412 if (nx <= 0 || ny <= 0) return 1; 413 414 *totaljzptr=0.0; 415 float sum1=0.0, sum2=0.0; 416 417 for (i = 0; i < nx; i++) 418 { 419 for (j = 0; j < ny; j++) 420 { 421 mbobra 1.1 if (mask[j * nx + i] <= 1) continue; 422 if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i]*18100000000.); 423 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i]*18100000000.); 424 } 425 } 426 427 *totaljzptr = fabs(sum1)+fabs(sum2); /* Units are A */ 428 return 0; 429 } 430 431 /*===========================================*/ 432 433 /* Example function 12: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */ 434 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV 435 // automatically yields erg per cubic centimeter for an input B in Gauss. 436 // 437 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert 438 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm: 439 // erg/cm^3(CDELT1)(RSUN_REF/RSUN_OBS)(100.) 440 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 441 // = erg/cm^3(1.30501e15) 442 mbobra 1.1 // = erg/cm(1/pix^2) 443 444 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, float *meanpotptr, float *totpotptr, int *mask) 445 { 446 int nx = dims[0], ny = dims[1]; 447 int i, j, count_mask=0; 448 449 if (nx <= 0 || ny <= 0) return 1; 450 451 *totpotptr=0.0; 452 *meanpotptr=0.0; 453 float sum=0.0; 454 455 for (i = 0; i < nx; i++) 456 { 457 for (j = 0; j < ny; j++) 458 { 459 if (mask[j * nx + i] <= 1) continue; 460 sum += (( ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) - ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i])) )/8.*PI); 461 count_mask++; 462 } 463 mbobra 1.1 } 464 465 *meanpotptr = (sum)/(count_mask); /* Units are ergs per cubic centimeter */ 466 *totpotptr = sum*(1.30501e15)*(count_mask); /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2, units of count_mask are pix^2; therefore, units of totpotptr are ergs per centimeter */ 467 return 0; 468 } 469 470 /*===========================================*/ 471 /* Example function 13: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */ 472 473 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, float *meanshear_angleptr, float *area_w_shear_gt_45ptr, float *meanshear_anglehptr, float *area_w_shear_gt_45hptr, int *mask) 474 { 475 int nx = dims[0], ny = dims[1]; 476 int i, j, count_mask=0; 477 478 if (nx <= 0 || ny <= 0) return 1; 479 480 *area_w_shear_gt_45ptr=0.0; 481 *meanshear_angleptr=0.0; 482 float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0; 483 float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0; 484 mbobra 1.1 485 for (i = 0; i < nx; i++) 486 { 487 for (j = 0; j < ny; j++) 488 { 489 if (mask[j * nx + i] <= 1) continue; 490 /* For mean 3D shear angle, area with shear greater than 45*/ 491 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]); 492 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])); 493 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]) ); 494 shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI); 495 sum += shear_angle ; 496 if (shear_angle > 45) count++; 497 498 /* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */ 499 dotproducth = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]); 500 magnitude_potentialh = sqrt((bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i])); 501 magnitude_vectorh = sqrt( (bx[j * nx + i]*bx[j * nx + i]) + (by[j * nx + i]*by[j * nx + i]) ); 502 shear_angleh = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI); 503 sum1 += shear_angleh ; 504 if (shear_angleh > 45) counth++; 505 mbobra 1.1 count_mask++; 506 } 507 } 508 509 /* For mean 3D shear angle, area with shear greater than 45*/ 510 *meanshear_angleptr = (sum)/(count_mask); /* Units are degrees */ 511 *area_w_shear_gt_45ptr = (count/(count_mask))*(100.); /* The area here is a fractional area -- the % of the total area */ 512 513 // /* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */ 514 // *meanshear_anglehptr = (sum1)/(count_mask); /* Units are degrees */ 515 // *area_w_shear_gt_45hptr = (counth/(count_mask))*(100.); /* The area here is a fractional area -- the % of the total area */ 516 return 0; 517 } 518 519 /*==================KEIJI'S CODE =========================*/ 520 521 // #include 522 #include 523 524 void greenpot(float *bx, float *by, float *bz, int nnx, int nny) 525 { 526 mbobra 1.1 /* local workings */ 527 int inx, iny, i, j, n; 528 /* local array */ 529 float *pfpot, *rdist; 530 pfpot=(float *)malloc(sizeof(float) *nnx*nny); 531 rdist=(float *)malloc(sizeof(float) *nnx*nny); 532 /* make nan */ 533 unsigned long long llnan = 0x7ff0000000000000; 534 // float NAN = (float)(llnan); 535 536 // #pragma omp parallel for private (inx) 537 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}} 538 // #pragma omp parallel for private (inx) 539 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}} 540 // #pragma omp parallel for private (inx) 541 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}} 542 // #pragma omp parallel for private (inx) 543 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}} 544 545 float dz; 546 dz = params_get_float(&cmdparams, "dzvalue"); 547 mbobra 1.1 548 // float dz = 0.0001; 549 550 // #pragma omp parallel for private (inx) 551 for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++) 552 { 553 float rdd, rdd1, rdd2; 554 float r; 555 rdd1 = (float)(inx); 556 rdd2 = (float)(iny); 557 rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz; 558 rdist[nnx*iny+inx] = 1.0/sqrt(rdd); 559 }} 560 561 int iwindow; 562 if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;} 563 float rwindow; 564 rwindow = (float)(iwindow); 565 rwindow = rwindow * rwindow + 0.01; 566 // #pragma omp parallel for private(inx) 567 for (iny=0;iny