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

Diff for /JSOC/proj/sharp/apps/sw_functions.c between version 1.13 and 1.19

version 1.13, 2013/05/31 22:47:15 version 1.19, 2013/10/18 23:36:02
 Line 63
 Line 63
//  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
//  =(1.30501e15)Gauss*cm^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, int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                    int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{ {

int nx = dims[0], ny = dims[1];      int nx = dims[0];
int i, j, count_mask=0;      int ny = dims[1];
double sum,err=0.0;      int i = 0;
int j = 0;
if (nx <= 0 || ny <= 0) return 1;      int count_mask = 0;
double sum = 0.0;
double err = 0.0;
*absFlux = 0.0;     *absFlux = 0.0;
*mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;

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

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{         {
for (j = 0; j < ny; j++)                 for (j = 0; j < ny; j++)
 Line 89  int computeAbsFlux(float *bz_err, float
 Line 89  int computeAbsFlux(float *bz_err, float
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;                   if isnan(bz[j * nx + i]) continue;
sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
//printf("i,j,bz[j * nx + i]=%d,%d,%f\n",i,j,bz[j * nx + i]);
err += bz_err[j * nx + i]*bz_err[j * nx + i];                   err += bz_err[j * nx + i]*bz_err[j * nx + i];
}                 }
 Line 97  int computeAbsFlux(float *bz_err, float
 Line 98  int computeAbsFlux(float *bz_err, float
*mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;      *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      *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
printf("USFLUX=%g\n",*mean_vf_ptr);       //printf("rsun_obs=%f\n",rsun_obs);
printf("sum=%f\n",sum);       //printf("rsun_ref=%f\n",rsun_ref);
//printf("USFLUX=%g\n",*mean_vf_ptr);
//printf("sum=%f\n",sum);
//printf("USFLUX_err=%g\n",*mean_vf_err_ptr);
return 0;      return 0;
} }

 Line 113  int computeBh(float *bx_err, float *by_e
 Line 117  int computeBh(float *bx_err, float *by_e

{ {

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

if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
 Line 146  int computeBh(float *bx_err, float *by_e
 Line 153  int computeBh(float *bx_err, float *by_e
int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims, int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
{ {
int nx = dims[0], ny = dims[1];      int nx = dims[0];
int i, j, count_mask=0;      int ny = dims[1];
int i = 0;
if (nx <= 0 || ny <= 0) return 1;      int j = 0;
double sum = 0.0;
double err = 0.0;
*mean_gamma_ptr=0.0;     *mean_gamma_ptr=0.0;
float sum,err,err_value=0.0;

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

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{           {
 Line 166  int computeGamma(float *bz_err, float *b
 Line 175  int computeGamma(float *bz_err, float *b
if isnan(bz_err[j * nx + i]) continue;                     if isnan(bz_err[j * nx + i]) continue;
if isnan(bh_err[j * nx + i]) continue;                     if isnan(bh_err[j * nx + i]) continue;
if (bz[j * nx + i] == 0) continue;                     if (bz[j * nx + i] == 0) continue;
sum += (atan(fabs(bz[j * nx + i]/bh[j * nx + i] )))*(180./PI);                      sum += fabs(atan(bh[j * nx + i]/fabs(bz[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);                     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);
}                   }
 Line 174  int computeGamma(float *bz_err, float *b
 Line 183  int computeGamma(float *bz_err, float *b
}           }

printf("MEANGAM=%f\n",*mean_gamma_ptr);       //printf("MEANGAM=%f\n",*mean_gamma_ptr);
printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);       //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
return 0;      return 0;
} }

 Line 187  int computeGamma(float *bz_err, float *b
 Line 196  int computeGamma(float *bz_err, float *b
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 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 nx = dims[0];
int i, j, count_mask=0;      int ny = dims[1];
int i = 0;
int j = 0;

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

 Line 212  int computeB_total(float *bx_err, float
 Line 224  int computeB_total(float *bx_err, float
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 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 nx = dims[0];
int i, j, count_mask=0;      int ny = dims[1];
int i = 0;
if (nx <= 0 || ny <= 0) return 1;      int j = 0;
double sum = 0.0;
double err = 0.0;
*mean_derivative_btotal_ptr = 0.0;     *mean_derivative_btotal_ptr = 0.0;
float sum, err = 0.0;

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

/* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
 Line 281  int computeBtotalderivative(float *bt, i
 Line 295  int computeBtotalderivative(float *bt, i

*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_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);          //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);          //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
return 0;         return 0;
} }

 Line 293  int computeBtotalderivative(float *bt, i
 Line 307  int computeBtotalderivative(float *bt, i
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 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 nx = dims[0];
int i, j, count_mask=0;       int ny = dims[1];
int i = 0;
int j = 0;
double sum= 0.0;
double err =0.0;
*mean_derivative_bh_ptr = 0.0;

if (nx <= 0 || ny <= 0) return 1;         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) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
{           {
 Line 361  int computeBhderivative(float *bh, float
 Line 378  int computeBhderivative(float *bh, float

*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_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);          //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);          //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);

return 0;         return 0;
} }
 Line 373  int computeBhderivative(float *bh, float
 Line 390  int computeBhderivative(float *bh, float
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 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 nx = dims[0];
int i, j, count_mask=0;          int ny = dims[1];
int i = 0;
int j = 0;
double sum = 0.0;
double err = 0.0;
*mean_derivative_bz_ptr = 0.0;

if (nx <= 0 || ny <= 0) return 1;         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) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
{           {
 Line 450  int computeBzderivative(float *bz, float
 Line 470  int computeBzderivative(float *bz, float

*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_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);          //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);          //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);

return 0;         return 0;
} }
 Line 498  int computeJz(float *bx_err, float *by_e
 Line 518  int computeJz(float *bx_err, float *by_e

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

if (nx <= 0 || ny <= 0) return 1;         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*/         /* Calculate the derivative*/
/* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */

for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
 Line 556  int computeJz(float *bx_err, float *by_e
 Line 576  int computeJz(float *bx_err, float *by_e
dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
}           }

for (i = 1; i <= nx-2; i++)
for (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j <= ny-1; j++)              for (j = 1; j <= ny-2; j++)
{             {
// calculate jz at all points                // 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[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]) +                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)]) ) ;                                             (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]);                jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);

}             }
}           }

return 0;         return 0;
} }

 Line 585  int computeJzsmooth(float *bx, float *by
 Line 605  int computeJzsmooth(float *bx, float *by

{ {

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

if (nx <= 0 || ny <= 0) return 1;         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*/         /* 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 (i = 0; i <= nx-1; i++)
{           {
 Line 609  int computeJzsmooth(float *bx, float *by
 Line 632  int computeJzsmooth(float *bx, float *by
}             }
}           }

/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */          /* Calculate mean vertical current density (mean_jz) 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_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         *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_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         *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=%f\n",*mean_jz_ptr);
printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);          //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);

printf("TOTUSJZ=%g\n",*us_i_ptr);          //printf("TOTUSJZ=%g\n",*us_i_ptr);
printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);          //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);

return 0;         return 0;

 Line 631  int computeJzsmooth(float *bx, float *by
 Line 654  int computeJzsmooth(float *bx, float *by
/* Example function 10:  Twist Parameter, alpha */ /* Example function 10:  Twist Parameter, alpha */

// The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation // 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):  // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):

// (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz         // numerator   = sum of all Jz*Bz
// (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz         // denominator = sum of Bz*Bz
// avg alpha = avg Jz / avg Bz         // alpha       = numerator/denominator

// 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 alpha are in 1/Mm
// The units of Jz are in Gauss/pix; the units of Bz are in Gauss. // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
 Line 656  int computeJzsmooth(float *bx, float *by
 Line 670  int computeJzsmooth(float *bx, float *by
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 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 nx                     = dims[0];
int i, j, count_mask, a,b,c,d=0;          int ny                     = dims[1];
int i                      = 0;
int j                      = 0;
double alpha_total         = 0.0;
double C                   = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
double total               = 0.0;
double A                   = 0.0;
double B                   = 0.0;

if (nx <= 0 || ny <= 0) return 1;         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 (i = 1; i < nx-1; i++)
{           {
for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
{               {
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 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(jz[j * nx + i]) continue;
if isnan(bz[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 (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.0) continue;
if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i] ); a++;                  //if (jz_err[j * nx + i] > abs(jz[j * nx + i]) ) continue;
if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;                  //if (bz_err[j * nx + i] > abs(bz[j * nx + i]) ) continue;
//if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);                  A += jz[j*nx+i]*bz[j*nx+i];
//if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);                  B += bz[j*nx+i]*bz[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.));
}               }
}           }

sum     = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */          for (i = 1; i < nx-1; i++)
{
/* Determine the sign of alpha */              for (j = 1; j < ny-1; j++)
if ((sum5 > 0) && (sum3 >  0)) sum=sum;                {
if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if ((sum5 < 0) && (sum4 <= 0)) sum=sum;                  if isnan(jz[j * nx + i])   continue;
if ((sum5 < 0) && (sum4 >  0)) sum=-sum;                  if isnan(bz[j * nx + i])   continue;
if (jz[j * nx + i] == 0.0) continue;
*mean_alpha_ptr = sum; /* Units are 1/Mm */                  if (bz[j * nx + i] == 0.0) continue;
*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                  total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
}
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);          /* Determine the absolute value of alpha. The units for alpha are 1/Mm */
printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);          alpha_total              = ((A/B)*C);
*mean_alpha_ptr          = alpha_total;
*mean_alpha_err_ptr      = (C/B)*(sqrt(total));

return 0;         return 0;
} }
 Line 725  int computeHelicity(float *jz_err, float
 Line 733  int computeHelicity(float *jz_err, float

{ {

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

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

float sum,sum2,sum_err=0.0;

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{         {
for (j = 0; j < ny; j++)                 for (j = 0; j < ny; j++)
 Line 752  int computeHelicity(float *jz_err, float
 Line 764  int computeHelicity(float *jz_err, float
*total_us_ih_ptr      = sum2           ; /* Units are G^2 / m ; keyword is TOTUSJH */         *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 */         *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          *mean_ih_err_ptr      = (sqrt(sum_err*sum_err)) / (count_mask*100.0)    ;  // error in the quantity MEANJZH
*total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity TOTUSJH          *total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity TOTUSJH
*total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity ABSNJZH          *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity ABSNJZH

printf("MEANJZH=%f\n",*mean_ih_ptr);          //printf("MEANJZH=%f\n",*mean_ih_ptr);
printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);          //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);

printf("TOTUSJH=%f\n",*total_us_ih_ptr);          //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);          //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);

printf("ABSNJZH=%f\n",*total_abs_ih_ptr);          //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);          //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);

return 0;         return 0;
} }
 Line 783  int computeSumAbsPerPolarity(float *jz_e
 Line 795  int computeSumAbsPerPolarity(float *jz_e

{ {
int nx = dims[0], ny = dims[1];          int nx = dims[0];
int i, j, count_mask=0;          int ny = dims[1];
int i=0;
int j=0;
double sum1=0.0;
double sum2=0.0;
double err=0.0;
*totaljzptr=0.0;

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

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

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{           {
for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
{               {
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
if isnan(jz[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) 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);                 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]);                 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
 Line 806  int computeSumAbsPerPolarity(float *jz_e
 Line 823  int computeSumAbsPerPolarity(float *jz_e

*totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */
*totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));         *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
printf("SAVNCPP=%g\n",*totaljzptr);          //printf("SAVNCPP=%g\n",*totaljzptr);
printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);          //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);

return 0;         return 0;
} }
 Line 830  int computeFreeEnergy(float *bx_err, flo
 Line 847  int computeFreeEnergy(float *bx_err, flo
float cdelt1, double rsun_ref, double rsun_obs)                                           float cdelt1, double rsun_ref, double rsun_obs)

{ {
int nx = dims[0], ny = dims[1];          int nx = dims[0];
int i, j, count_mask=0;          int ny = dims[1];
int i = 0;
if (nx <= 0 || ny <= 0) return 1;          int j = 0;
double sum = 0.0;
double sum1 = 0.0;
double err = 0.0;
*totpotptr=0.0;         *totpotptr=0.0;
*meanpotptr=0.0;         *meanpotptr=0.0;
float sum,sum1,err=0.0;
if (nx <= 0 || ny <= 0) return 1;

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{           {
 Line 860  int computeFreeEnergy(float *bx_err, flo
 Line 881  int computeFreeEnergy(float *bx_err, flo
*totpotptr       = (sum)/(8.*PI);         *totpotptr       = (sum)/(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)));         *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=%g\n",*meanpotptr);
printf("MEANPOT_err=%g\n",*meanpot_err_ptr);          //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);

printf("TOTPOT=%g\n",*totpotptr);          //printf("TOTPOT=%g\n",*totpotptr);
printf("TOTPOT_err=%g\n",*totpot_err_ptr);          //printf("TOTPOT_err=%g\n",*totpot_err_ptr);

return 0;         return 0;
} }
 Line 875  int computeFreeEnergy(float *bx_err, flo
 Line 896  int computeFreeEnergy(float *bx_err, flo
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, 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,
{ {
int nx = dims[0], ny = dims[1];          int nx = dims[0];
int i, j;          int ny = dims[1];
int i = 0;
int j = 0;
double dotproduct = 0.0;
double magnitude_potential = 0.0;
double magnitude_vector = 0.0;
double shear_angle = 0.0;
double err = 0.0;
double sum = 0.0;
double count = 0.0;
*area_w_shear_gt_45ptr = 0.0;
*meanshear_angleptr = 0.0;

if (nx <= 0 || ny <= 0) return 1;         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 (i = 0; i < nx; i++)
{           {
for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
 Line 910  int computeShearAngle(float *bx_err, flo
 Line 939  int computeShearAngle(float *bx_err, flo
/* For mean 3D shear angle, area with shear greater than 45*/         /* For mean 3D shear angle, area with shear greater than 45*/
*meanshear_angleptr = (sum)/(count);                 /* Units are degrees */         *meanshear_angleptr = (sum)/(count);                 /* Units are degrees */
*meanshear_angle_err_ptr = (sqrt(err*err))/(count);  // error in the quantity (sum)/(count_mask)         *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 */          *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */

printf("MEANSHR=%f\n",*meanshear_angleptr);          //printf("MEANSHR=%f\n",*meanshear_angleptr);
printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);          //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);

return 0;         return 0;
} }
 Line 1034  void greenpot(float *bx, float *by, floa
 Line 1063  void greenpot(float *bx, float *by, floa

/*===========END OF KEIJI'S CODE =========================*/ /*===========END OF KEIJI'S CODE =========================*/

char *sw_functions_version() // Returns CVS version of sw_functions.c
{
return strdup("\$Id\$");
}

/* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

Legend:
 Removed from v.1.13 changed lines Added in v.1.19