Return to
sw_functions.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps
|
Return to
sw_functions.c
CVS log
|
Up to [Development] / JSOC / proj / sharp / apps
|
|
File: [Development] / JSOC / proj / sharp / apps / sw_functions.c
(download)
Revision: 1.10, Mon May 20 23:01:34 2013 UTC (10 years ago) by mbobra Branch: MAIN Changes since 1.9: +60 -295 lines Oops, last commit used a 6th-order finite difference method. I replace it with a second-order finite difference method. |
/*===========================================
The following 14 functions calculate the following spaceweather indices:
USFLUX Total unsigned flux in Maxwells
MEANGAM Mean inclination angle, gamma, in degrees
MEANGBT Mean value of the total field gradient, in Gauss/Mm
MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
MEANJZD Mean vertical current density, in mA/m2
TOTUSJZ Total unsigned vertical current, in Amperes
MEANALP Mean twist parameter, alpha, in 1/Mm
MEANJZH Mean current helicity in G2/m
TOTUSJH Total unsigned current helicity in G2/m
ABSNJZH Absolute value of the net current helicity in G2/m
SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
MEANSHR Mean shear angle (measured using Btotal) in degrees
The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
coordinate bitmaps are interpolated.
In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:
For conf_disambig:
50 : not all solutions agree (weak field method applied)
60 : not all solutions agree (weak field + annealed)
90 : all solutions agree (strong field + annealed)
0 : not disambiguated
For bitmap:
1 : weak field outside smooth bounding curve
2 : strong field outside smooth bounding curve
33 : weak field inside smooth bounding curve
34 : strong field inside smooth bounding curve
Written by Monica Bobra 15 August 2012
Potential Field code (appended) written by Keiji Hayashi
===========================================*/
#include <math.h>
#include <mkl.h>
#define PI (M_PI)
#define MUNAUGHT (0.0000012566370614) /* magnetic constant */
/*===========================================*/
/* Example function 1: Compute total unsigned flux in units of G/cm^2 */
// To compute the unsigned flux, we simply calculate
// flux = surface integral [(vector Bz) dot (normal vector)],
// = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
// However, since the field is radial, we will assume cos theta = 1.
// Therefore the pixels only need to be corrected for the projection.
// 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/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
// =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
// =(1.30501e15)Gauss*cm^2
// The disambig mask value selects only the pixels with values of 5 or 7 -- that is,
// 5: pixels for which the radial acute disambiguation solution was chosen
// 7: pixels for which the radial acute and NRWA disambiguation agree
int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
double sum,err=0.0;
if (nx <= 0 || ny <= 0) return 1;
*absFlux = 0.0;
*mean_vf_ptr =0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;
sum += (fabs(bz[j * nx + i]));
err += bz_err[j * nx + i]*bz_err[j * nx + i];
count_mask++;
}
}
*mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
*mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux
*count_mask_ptr = count_mask;
printf("CMASK=%g\n",*count_mask_ptr);
printf("USFLUX=%g\n",*mean_vf_ptr);
printf("sum=%f\n",sum);
printf("USFLUX_err=%g\n",*mean_vf_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
// Native units of Bh are Gauss
int computeBh(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
float *mean_hf_ptr, int *mask, int *bitmask)
{
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 (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if isnan(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) continue;
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];
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];
count_mask++;
}
}
*mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
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)
// Redo calculation in radians for error analysis (since derivatives are only true in units of radians).
int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
{
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,err,err_value=0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if (bh[j * nx + i] > 100)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;
if isnan(bz_err[j * nx + i]) continue;
if isnan(bh_err[j * nx + i]) continue;
if (bz[j * nx + i] == 0) continue;
sum += (atan(fabs(bz[j * nx + i]/bh[j * nx + i] )))*(180./PI);
err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i]))) * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);
count_mask++;
}
}
}
*mean_gamma_ptr = sum/count_mask;
*mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.); // error in the quantity (sum)/(count_mask)
printf("MEANGAM=%f\n",*mean_gamma_ptr);
printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 4: Calculate B_Total*/
// Native units of B_Total are in gauss
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)
{
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++)
{
if isnan(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;
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]);
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];
}
}
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 *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr)
{
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 sum, err = 0.0;
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)
{
for (j = 0; j <= ny-1; j++)
{
derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
}
}
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)
{
for (j = 1; j <= ny-2; j++)
{
dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
}
}
/* consider the edges */
i=0;
for (j = 0; j <= ny-1; j++)
{
derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
}
i=nx-1;
for (j = 0; j <= ny-1; j++)
{
derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
}
j=0;
for (i = 0; i <= nx-1; i++)
{
dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
}
j=ny-1;
for (i = 0; i <= nx-1; i++)
{
dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
}
for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(derx_bt[j * nx + i]) continue;
if isnan(dery_bt[j * nx + i]) continue;
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 */
err += (2.0)*bt_err[j * nx + i]*bt_err[j * nx + i];
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
*mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 6: Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
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)
{
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 sum,err = 0.0;
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)
{
for (j = 0; j <= ny-1; j++)
{
derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
}
}
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)
{
for (j = 1; j <= ny-2; j++)
{
dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
}
}
/* consider the edges */
i=0;
for (j = 0; j <= ny-1; j++)
{
derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
}
i=nx-1;
for (j = 0; j <= ny-1; j++)
{
derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
}
j=0;
for (i = 0; i <= nx-1; i++)
{
dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
}
j=ny-1;
for (i = 0; i <= nx-1; i++)
{
dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
}
for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(derx_bh[j * nx + i]) continue;
if isnan(dery_bh[j * nx + i]) continue;
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 */
err += (2.0)*bh_err[j * nx + i]*bh_err[j * nx + i];
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
*mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 7: Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
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)
{
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 sum,err = 0.0;
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)
{
for (j = 0; j <= ny-1; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
}
}
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)
{
for (j = 1; j <= ny-2; j++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
}
}
/* consider the edges */
i=0;
for (j = 0; j <= ny-1; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
}
i=nx-1;
for (j = 0; j <= ny-1; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
}
j=0;
for (i = 0; i <= nx-1; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
}
j=ny-1;
for (i = 0; i <= nx-1; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
}
for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
// if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;
//if isnan(bz_err[j * nx + i]) continue;
if isnan(derx_bz[j * nx + i]) continue;
if isnan(dery_bz[j * nx + i]) continue;
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 */
err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];
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
*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
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)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
// (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
// (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
// 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/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)
// = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
// Comment out random number generator, which can only run on solar3
//int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
// int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
// float *noiseby, float *noisebz)
int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
if (nx <= 0 || ny <= 0) return 1;
float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
/* Calculate the derivative*/
/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)
{
for (j = 0; j <= ny-1; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
}
}
for (i = 0; i <= nx-1; i++)
{
for (j = 1; j <= ny-2; j++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
}
}
// consider the edges
i=0;
for (j = 0; j <= ny-1; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
}
i=nx-1;
for (j = 0; j <= ny-1; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
}
j=0;
for (i = 0; i <= nx-1; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
}
j=ny-1;
for (i = 0; i <= nx-1; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
}
for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
// calculate jz at all points
jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); // jz is in units of Gauss/pix
jz_err[j * nx + i]=0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * 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)]) ) ;
jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);
count_mask++;
}
}
return 0;
}
/*===========================================*/
/* Example function 9: Compute quantities on smoothed Jz array */
// All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's
// fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels
// and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values
// of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course,
// give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels.
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,
float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
if (nx <= 0 || ny <= 0) return 1;
float curl,us_i,test_perimeter,mean_curl,err=0.0;
/* 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*/
for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
//printf("%f ",jz_smooth[j * nx + i]);
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(derx[j * nx + i]) continue;
if isnan(dery[j * nx + i]) continue;
//if isnan(jz_smooth[j * nx + i]) continue;
//curl += (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
//us_i += fabs(jz_smooth[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */
//jz_rms_err[j * nx + i] = sqrt(jz_err_squared_smooth[j * nx + i]);
//err += (jz_rms_err[j * nx + i]*jz_rms_err[j * nx + i]);
if isnan(jz[j * nx + i]) continue;
curl += (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A */
err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
count_mask++;
}
//printf("\n");
}
/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
*mean_jz_ptr = curl/(count_mask); /* mean_jz gets populated as MEANJZD */
*mean_jz_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // error in the quantity MEANJZD
*us_i_ptr = (us_i); /* us_i gets populated as TOTUSJZ */
*us_i_err_ptr = (sqrt(err))*fabs((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
printf("MEANJZD=%f\n",*mean_jz_ptr);
printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
printf("TOTUSJZ=%g\n",*us_i_ptr);
printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 10: Twist Parameter, alpha */
// The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
// for alpha is calculated in the following way (different from Leka and Barnes' approach):
// (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
// (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
// avg alpha = avg Jz / avg Bz
// The sign is assigned as follows:
// If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
// If this value is > 0, then alpha is > 0.
// If this value is < 0, then alpha is <0.
//
// If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
// If this value is > 0, then alpha is < 0.
// If this value is < 0, then alpha is > 0.
// The units of alpha are in 1/Mm
// The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
//
// Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
// = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
// = 1/Mm
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)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask, a,b,c,d=0;
if (nx <= 0 || ny <= 0) return 1;
float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6, sum_err=0.0;
for (i = 1; i < nx-1; i++)
{
for (j = 1; j < ny-1; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
//if isnan(jz_smooth[j * nx + i]) continue;
if isnan(jz[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;
//if (jz_smooth[j * nx + i] == 0) continue;
if (jz[j * nx + i] == 0.0) continue;
if (bz_err[j * nx + i] == 0.0) continue;
if (bz[j * nx + i] == 0.0) continue;
if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i] ); a++;
if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;
//if (bz[j * nx + i] > 0) sum3 += ( jz_smooth[j * nx + i]);
//if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
if (bz[j * nx + i] > 0) sum3 += ( jz[j * nx + i] ); c++;
if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;
sum5 += bz[j * nx + i];
/* sum_err is a fractional uncertainty */
sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs( ( (jz[j * nx + i]) / (bz[j * nx + i]) ) *(1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
count_mask++;
}
}
sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */
/* Determine the sign of alpha */
if ((sum5 > 0) && (sum3 > 0)) sum=sum;
if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
if ((sum5 < 0) && (sum4 > 0)) sum=-sum;
*mean_alpha_ptr = sum; /* Units are 1/Mm */
*mean_alpha_err_ptr = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent
printf("a=%d\n",a);
printf("b=%d\n",b);
printf("d=%d\n",d);
printf("c=%d\n",c);
printf("MEANALP=%f\n",*mean_alpha_ptr);
printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 11: Helicity (mean current helicty, total unsigned current helicity, absolute value of net 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)(pix/arcsec)(arcsec/meter)
// = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
// = G^2 / m.
int computeHelicity(float *jz_err, float *jz_rms_err, float *bz_err, float *bz, int *dims, float *jz, float *mean_ih_ptr,
float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
{
int nx = dims[0], ny = dims[1];
int i, j, count_mask=0;
if (nx <= 0 || ny <= 0) return 1;
float sum,sum2,sum_err=0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));
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 ; keyword is TOTUSJH */
*total_abs_ih_ptr = fabs(sum) ; /* Units are G^2 / m ; keyword is ABSNJZH */
*mean_ih_err_ptr = (sqrt(sum_err*sum_err)) / (count_mask*100.) ; // error in the quantity MEANJZH
*total_us_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.) ; // error in the quantity TOTUSJH
*total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.) ; // error in the quantity ABSNJZH
printf("MEANJZH=%f\n",*mean_ih_ptr);
printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
printf("TOTUSJH=%f\n",*total_us_ih_ptr);
printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 12: 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:
// Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
// = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
//
// The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
{
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,sum2,err=0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;
if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
count_mask++;
}
}
*totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */
*totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
printf("SAVNCPP=%g\n",*totaljzptr);
printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
return 0;
}
/*===========================================*/
/* Example function 13: 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^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
// = 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_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_ptr, int *mask, int *bitmask,
float cdelt1, double rsun_ref, double rsun_obs)
{
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,err=0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) continue;
sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) ) / 8.*PI;
err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
//err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];
count_mask++;
}
}
*meanpotptr = (sum) / (count_mask); /* Units are ergs per cubic centimeter */
*meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
//*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
/* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
*totpotptr = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/8.*PI)) ;
*totpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/8.*PI));
printf("MEANPOT=%g\n",*meanpotptr);
printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
printf("TOTPOT=%g\n",*totpotptr);
printf("TOTPOT_err=%g\n",*totpot_err_ptr);
return 0;
}
/*===========================================*/
/* 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 */
int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
{
int nx = dims[0], ny = dims[1];
int i, j;
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,err=0.0, sum = 0.0, count=0.0, count_mask=0.0;
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bpx[j * nx + i]) continue;
if isnan(bpy[j * nx + i]) continue;
if isnan(bpz[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;
if isnan(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) 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);
count ++;
sum += shear_angle ;
err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
if (shear_angle > 45) count_mask ++;
}
}
/* For mean 3D shear angle, area with shear greater than 45*/
*meanshear_angleptr = (sum)/(count); /* Units are degrees */
*meanshear_angle_err_ptr = (sqrt(err*err))/(count); // error in the quantity (sum)/(count_mask)
*area_w_shear_gt_45ptr = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */
printf("MEANSHR=%f\n",*meanshear_angleptr);
printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
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);
float *bztmp;
bztmp=(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;}}
// #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)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
}}
// dz is the monopole depth
float dz = 0.001;
// #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; // must be of square
rwindow = 1.0e2; // limit the window size to be 10.
rwindow = sqrt(rwindow);
iwindow = (int)(rwindow);
// #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] = 0.0; // hmmm.. NAN;
}
else
{
float sum;
sum = 0.0;
int j2, i2;
int j2s, j2e, i2s, i2e;
j2s = iny - iwindow;
j2e = iny + iwindow;
if (j2s < 0){j2s = 0;}
if (j2e > nny){j2e = nny;}
i2s = inx - iwindow;
i2e = inx + iwindow;
if (i2s < 0){i2s = 0;}
if (i2e > nnx){i2e = nnx;}
for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
{
float val1 = bztmp[nnx*j2 + i2];
float rr, r1, r2;
// r1 = (float)(i2 - inx);
// r2 = (float)(j2 - iny);
// rr = r1*r1 + r2*r2;
// if (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)]) * 0.5;
by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
} } // end of OpenMP parallelism
free(rdist);
free(pfpot);
free(bztmp);
} // end of void func. greenpot
/*===========END OF KEIJI'S CODE =========================*/
/* ---------------- end of this file ----------------*/
| Karen Tian |
Powered by ViewCVS 0.9.4 |