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

 File: [Development] / JSOC / proj / sharp / apps / sharp_functions.c (download) Revision: 1.1, Wed Apr 4 21:55:03 2012 UTC (10 years, 9 months ago) by mbobra Branch: MAIN CVS Tags: Ver_LATEST, Ver_9-5, Ver_9-41, Ver_9-4, Ver_9-3, Ver_9-2, Ver_9-1, Ver_9-0, Ver_8-8, Ver_8-7, Ver_8-6, Ver_8-5, Ver_8-4, Ver_8-3, Ver_8-2, Ver_8-12, Ver_8-11, Ver_8-10, Ver_8-1, Ver_8-0, Ver_7-1, Ver_7-0, Ver_6-4, Ver_6-3, Ver_6-2, HEAD ```functions for the space weather hmi active region patch (sharp) computations ```

```/* sharp_functions.c */
#define PI (3.141592653589793)

/*===========================================*/

/* Example function 1: Compute total unsigned flux in units of G/cm^2, Mean Vertical Field */

//  To convert G to G/cm^2, simply multiply by the number of square centimeters per pixel.
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  Gauss(CDELT1)(RSUN_REF/RSUN_OBS)(100.)
//  = Gauss(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
//  = Gauss(1.30501e15)

int computeAbsFlux(float *bz, int *dims, float *absFlux, float *mean_vf_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
double sum=0.0;

if (nx <= 0 || ny <= 0) return 1;

*absFlux = 0.0;
*mean_vf_ptr =0.0;

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)
{
if (mask[j * nx + i] <= 1) continue;
sum += (fabs(bz[j * nx + i]));
count_mask++;
}
}

printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum);
*mean_vf_ptr = sum*(1.30501e15);
return 0;
}

/*===========================================*/
/* Example function 2: Calculate Bh in units of Gauss */
// Native units of Bh are Gauss

int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_hf_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
float sum=0.0;
*mean_hf_ptr =0.0;

if (nx <= 0 || ny <= 0) return 1;

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)
{
bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
sum += bh[j * nx + i];
count_mask++;
}
}

*mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
printf("*mean_hf_ptr=%f\n",*mean_hf_ptr);
return 0;
}

/*===========================================*/

/* Example function 3: Calculate Gamma in units of degrees
// Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)

this is wrong, also the divide by nx*ny is wrong */

int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, float *mean_gamma_ptr, int *mask)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_gamma_ptr=0.0;
float sum=0.0;
int count=0;

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if (bh[j * nx + i] > 100)
{
if (mask[j * nx + i] <= 1) continue;
sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
count_mask++;
}
}
}

*mean_gamma_ptr = sum/count_mask;
printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);
return 0;
}

/*===========================================*/

/* Example function 4: Calculate B_Total*/
// Native units of B_Total are in gauss

int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
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]);
}
}
return 0;
}

/*===========================================*/

/* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */

int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_derivative_btotal_ptr = 0.0;
float derx, dery, sum = 0.0;

for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if (mask[j * nx + i] <= 1) continue;
/* Missing a (*0.5) */
derx = bt[j * nx + i+1] - bt[j * nx + i-1];
dery = bt[(j+1) * nx + i] - bt[(j-1) * nx + i];
sum += sqrt(derx*derx + dery*dery);
count_mask++;
}
}

*mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);
return 0;
}

/*===========================================*/
/* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */

int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_derivative_bh_ptr = 0.0;
float derx, dery, sum = 0.0;

for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if (mask[j * nx + i] <= 1) continue;
/* Missing a (*0.5) */
derx = bh[j * nx + i+1] - bh[j * nx + i-1];
//printf("derx=%f\n",derx);
dery = bh[(j+1) * nx + i] - bh[(j-1) * nx + i];
//                printf("dery=%f\n",dery);
sum += sqrt(derx*derx + dery*dery);
//printf("sum=%f\n",sum);
count_mask++;
//printf("count_mask=%d\n",count_mask);
}
}

*mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
printf("*mean_derivative_bh_ptr=%f,nx=%d,ny=%d,sum=%f\n",*mean_derivative_bh_ptr,nx,ny,sum);
return 0;
}

/*===========================================*/
/* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */

int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_derivative_bz_ptr = 0.0;
float derx, dery, sum = 0.0;

for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if (mask[j * nx + i] <= 1) continue;
/* Missing a (*0.5) */
derx = bz[j * nx + i+1] - bz[j * nx + i-1];
dery = bz[(j+1) * nx + i] - bz[(j-1) * nx + i];
sum += sqrt(derx*derx + dery*dery);
count_mask++;
}
}

*mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram

return 0;
}

/*===========================================*/
/* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */

//  In discretized space like data pixels,
//  the current (or curl of B) is calculated as
//  the integration of the field Bx and By along
//  the circumference of the data pixel divided by the area of the pixel.
//  One form of differencing (a word for the differential operator
//  in the discretized space) of the curl is expressed as
//  (dx * (Bx(i,j-1)+Bx(i,j)) / 2
//  +dy * (By(i+1,j)+By(i,j)) / 2
//  -dx * (Bx(i,j+1)+Bx(i,j)) / 2
//  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
//
//
//  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
//  one must perform the following unit conversions:
//  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
//  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
//  where a Tesla is represented as a Newton/Ampere*meter.
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  In that case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
//  jz * (35.0)
//
//  The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or
//  (Gauss/pix)(1.81*10^10)
//  (Gauss/pix)(18100000000.)

int computeJz(float *bx, float *by, int *dims, float *jz, float *mean_jz_ptr, float *us_i_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_jz_ptr = 0.0;
float derx, dery, curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;

for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if (mask[j * nx + i] <= 1) continue;
derx = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
dery = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
curl += (derx-dery)*(35.0);               /* curl is in units of mA/m^2 */
jz[j * nx + i] = (derx-dery);             /* jz is in units of Gauss per pixel */
us_i += fabs(derx-dery)*(18100000000.) ;  /* us_i is in units of A */
count_mask++;
}
}

mean_curl      = (curl/count_mask);
printf("mean_curl=%f\n",mean_curl);
*mean_jz_ptr     = curl/(count_mask);
printf("count_mask=%d\n",count_mask);
*us_i_ptr = (us_i);
return 0;
}

/*===========================================*/
/* Example function 9:  Twist Parameter, alpha */

// The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm
// The units of Jz are in G/pix; the units of Bz are in G.
// Therefore, the units of Jz/Bz = (pix/arcsec)(arcsec/meter)(meter/Mm), or
//                         Jz/Bz = (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6).
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  In that case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^6), or
//  (Gauss/pix)*2.7

int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_alpha_ptr = 0.0;
float aa, bb, cc, bznew, alpha2, sum=0.0;

for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if (mask[j * nx + i] <= 1) continue;
if isnan(jz[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;
if (bz[j * nx + i] == 0.0) continue;
sum += (jz[j * nx + i] / bz[j * nx + i])*2.7 ; /* the units for (jz) Gauss/pix; the units for bz are Gauss */
//printf("sum=%f\n",sum);
//printf("jz[j * nx + i]=%f\n",jz[j * nx + i]);
//printf("bz[j * nx + i]=%f\n",bz[j * nx + i]);
//printf("(jz[j * nx + i] / bz[j * nx + i])*2.7=%f\n",(jz[j * nx + i] / bz[j * nx + i])*2.7);
count_mask++;
//printf("count_mask=%d\n",count_mask);
}
}

printf("count_mask=%d\n",count_mask);
printf("sum=%f\n",sum);
*mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */
return 0;
}

/*===========================================*/

/* Example function 10:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */

//  The current helicity is defined as Bz*Jz and the units are G^2 / m
//  The units of Jz are in G/pix; the units of Bz are in G.
//  Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF) = G^2 / m.
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  In that case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500), or
//  (Gauss/pix)*0.000002768166

int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*mean_ih_ptr = 0.0;
float sum=0.0, sum2=0.0;

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)
{
if (mask[j * nx + i] <= 1) continue;
if isnan(jz[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;
if (bz[j * nx + i] == 0.0) continue;
if (jz[j * nx + i] == 0.0) continue;
sum  += (jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166);
sum2 += fabs(jz[j * nx + i]*(bz[j * nx + i]))*(0.000002768166);
count_mask++;
printf("(jz[j * nx + i])=%f\n",(jz[j * nx + i]));
printf("(bz[j * nx + i])=%f\n",(bz[j * nx + i]));
printf("(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166)=%f\n",(jz[j * nx + i])*(bz[j * nx + i])*(0.000002768166));
printf("sum=%f\n",sum);
printf("sum2=%f\n",sum2);
printf("count_mask=%d\n",count_mask);
}
}

printf("sum/count_mask=%f\n",sum/count_mask);
*mean_ih_ptr     = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
*total_us_ih_ptr = sum2;           /* Units are G^2 / m */
*total_abs_ih_ptr= fabs(sum);      /* Units are G^2 / m */

return 0;
}

/*===========================================*/

/* Example function 11:  Sum of Absolute Value per polarity */

//  The Sum of the Absolute Value per polarity is defined as the following:
//  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
//  The units of jz are in G/pix. In this case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(722500)(722500), or
//  (Gauss/pix)(1.81*10^10)
//  (Gauss/pix)(18100000000.)

int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int *mask)
{

int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*totaljzptr=0.0;
float sum1=0.0, sum2=0.0;

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if (mask[j * nx + i] <= 1) continue;
if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i]*18100000000.);
if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i]*18100000000.);
}
}

*totaljzptr = fabs(sum1)+fabs(sum2);  /* Units are A */
return 0;
}

/*===========================================*/

/* Example function 12:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
// The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
// automatically yields erg per cubic centimeter for an input B in Gauss.
//
// Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
// ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
// erg/cm^3(CDELT1)(RSUN_REF/RSUN_OBS)(100.)
// = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
// = erg/cm^3(1.30501e15)
// = erg/cm(1/pix^2)

int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, float *meanpotptr, float *totpotptr, int *mask)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*totpotptr=0.0;
*meanpotptr=0.0;
float sum=0.0;

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if (mask[j * nx + i] <= 1) continue;
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);
count_mask++;
}
}

*meanpotptr = (sum)/(count_mask);              /* Units are ergs per cubic centimeter */
*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 */
return 0;
}

/*===========================================*/
/* 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 */

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)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;

if (nx <= 0 || ny <= 0) return 1;

*area_w_shear_gt_45ptr=0.0;
*meanshear_angleptr=0.0;
float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0;
float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if (mask[j * nx + i] <= 1) continue;
/* For mean 3D shear angle, area with shear greater than 45*/
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]);
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]));
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]) );
shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
sum += shear_angle ;
if (shear_angle > 45) count++;

/* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */
dotproducth           = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]);
magnitude_potentialh  = sqrt((bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]));
magnitude_vectorh     = sqrt( (bx[j * nx + i]*bx[j * nx + i]) + (by[j * nx + i]*by[j * nx + i]) );
shear_angleh          = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
sum1 += shear_angleh ;
if (shear_angleh > 45) counth++;
count_mask++;
}
}

/* For mean 3D shear angle, area with shear greater than 45*/
*meanshear_angleptr = (sum)/(count_mask);              /* Units are degrees */
*area_w_shear_gt_45ptr = (count/(count_mask))*(100.);  /* The area here is a fractional area -- the % of the total area */

//        /* For mean horizontal shear angle, area with horizontal shear angle greater than 45 */
//        *meanshear_anglehptr = (sum1)/(count_mask);              /* Units are degrees */
//        *area_w_shear_gt_45hptr = (counth/(count_mask))*(100.);  /* The area here is a fractional area -- the % of the total area */
return 0;
}

/*==================KEIJI'S CODE =========================*/

// #include <omp.h>
#include <math.h>

void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
{
/* local workings */
int inx, iny, i, j, n;
/* local array */
float *pfpot, *rdist;
pfpot=(float *)malloc(sizeof(float) *nnx*nny);
rdist=(float *)malloc(sizeof(float) *nnx*nny);
/* make nan */
unsigned long long llnan = 0x7ff0000000000000;
//  float NAN = (float)(llnan);

// #pragma omp parallel for private (inx)
for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}

float dz;
dz = params_get_float(&cmdparams, "dzvalue");

// float dz = 0.0001;

// #pragma omp parallel for private (inx)
for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
{
float rdd, rdd1, rdd2;
float r;
rdd1 = (float)(inx);
rdd2 = (float)(iny);
rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
}}

int iwindow;
if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
float rwindow;
rwindow = (float)(iwindow);
rwindow = rwindow * rwindow + 0.01;
// #pragma omp parallel for private(inx)
for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
{
float val0 = bz[nnx*iny + inx];
if (isnan(val0))
{
pfpot[nnx*iny + inx] = NAN;
}
else
{
float sum;
sum = 0.0;
int j2, i2;
for (j2=0;j2<nny;j2++){for (i2=0;i2<nnx;i2++)
{
float val1 = bz[nnx*j2 + i2];
float rr, r1, r2;
r1 = (float)(i2 - inx);
r2 = (float)(j2 - iny);
rr = r1*r1 + r2*r2;
if ((!isnan(val1)) && (rr < rwindow))
{
int   di, dj;
di = abs(i2 - inx);
dj = abs(j2 - iny);
sum = sum + val1 * rdist[nnx * dj + di] * dz;
}
} }
pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
}
} } // end of OpenMP parallelism

// #pragma omp parallel for private(inx)
for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
{
bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) / 2.0;
by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) / 2.0;
} }
free(pfpot);
} // end of void func. greenpot

/*===========END OF KEIJI'S CODE =========================*/
/* ---------------- end of this file ----------------*/
```

 Karen Tian Powered byViewCVS 0.9.4