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.33 and 1.36

version 1.33, 2015/01/23 01:18:13 version 1.36, 2015/10/30 18:53:02
 Line 196  int computeGamma(float *bz_err, float *b
 Line 196  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);
printf("bh[200,200],bh_err[200,200]=%f,%f\n",bh[200,200],bh_err[200,200]);

return 0;     return 0;
} }

 Line 246  int computeB_total(float *bx_err, float
 Line 249  int computeB_total(float *bx_err, float
/*===========================================*/ /*===========================================*/
/* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ /* 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 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, float *err_term1, float *err_term2)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 266  int computeBtotalderivative(float *bt, i
 Line 269  int computeBtotalderivative(float *bt, i
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{         {
derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;             derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
}         }
}     }

 Line 275  int computeBtotalderivative(float *bt, i
 Line 279  int computeBtotalderivative(float *bt, i
for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
{         {
dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;             dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
}         }
}     }

/* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
/* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
i=0;     i=0;
for (j = 0; j <= ny-1; j++)     for (j = 0; j <= ny-1; j++)
{     {
 Line 304  int computeBtotalderivative(float *bt, i
 Line 310  int computeBtotalderivative(float *bt, i
dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;         dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
}     }

// Calculate the sum only
for (i = 1; i <= nx-2; i++)     for (i = 1; i <= nx-2; i++)
{     {
for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
 Line 320  int computeBtotalderivative(float *bt, i
 Line 326  int computeBtotalderivative(float *bt, i
if isnan(derx_bt[j * nx + i]) continue;             if isnan(derx_bt[j * nx + i]) continue;
if isnan(dery_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 */             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 += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+              err += err_term2[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+
(((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;                     err_term1[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
}         }
}     }
 Line 338  int computeBtotalderivative(float *bt, i
 Line 344  int computeBtotalderivative(float *bt, i
/*===========================================*/ /*===========================================*/
/* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ /* 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 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, float *err_term1, float *err_term2)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 358  int computeBhderivative(float *bh, float
 Line 364  int computeBhderivative(float *bh, float
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{         {
derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;             derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
}         }
}     }

 Line 367  int computeBhderivative(float *bh, float
 Line 374  int computeBhderivative(float *bh, float
for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
{         {
dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;             dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
}         }
}     }

/* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
/* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
i=0;     i=0;
for (j = 0; j <= ny-1; j++)     for (j = 0; j <= ny-1; j++)
{     {
 Line 412  int computeBhderivative(float *bh, float
 Line 421  int computeBhderivative(float *bh, float
if isnan(derx_bh[j * nx + i]) continue;             if isnan(derx_bh[j * nx + i]) continue;
if isnan(dery_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 */             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 += (((bh[(j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+              err += err_term2[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+
(((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;                     err_term1[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;
}         }
}     }
 Line 429  int computeBhderivative(float *bh, float
 Line 438  int computeBhderivative(float *bh, float
/*===========================================*/ /*===========================================*/
/* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ /* 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 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, float *err_term1, float *err_term2)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 448  int computeBzderivative(float *bz, float
 Line 457  int computeBzderivative(float *bz, float
{     {
for (j = 0; j <= ny-1; j++)             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;             derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
}         }
}     }

 Line 458  int computeBzderivative(float *bz, float
 Line 467  int computeBzderivative(float *bz, float
{     {
for (j = 1; j <= ny-2; j++)             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;             dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
}         }
}     }

/* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
/* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
i=0;     i=0;
for (j = 0; j <= ny-1; j++)     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;         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;     i=nx-1;
for (j = 0; j <= ny-1; j++)     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;         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;     j=0;
for (i = 0; i <= nx-1; i++)     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;         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;     j=ny-1;
for (i = 0; i <= nx-1; i++)     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;         dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
}     }

 Line 509  int computeBzderivative(float *bz, float
 Line 515  int computeBzderivative(float *bz, float
if isnan(derx_bz[j * nx + i]) continue;             if isnan(derx_bz[j * nx + i]) continue;
if isnan(dery_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 */             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 += (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i])) /              err += err_term2[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +
(16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +                     err_term1[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
(((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)])) /
(16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
}         }
}     }
 Line 582  int computeJz(float *bx_err, float *by_e
 Line 586  int computeJz(float *bx_err, float *by_e
{     {
for (j = 0; j <= ny-1; j++)             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;             derx[j * nx + i]      = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);             err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);
}         }
 Line 592  int computeJz(float *bx_err, float *by_e
 Line 595  int computeJz(float *bx_err, float *by_e
{     {
for (j = 1; j <= ny-2; j++)             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;             dery[j * nx + i]      = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);             err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
}         }
}     }

// consider the edges      /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
ignore the edges for the error terms as those arrays have been initialized to zero.
this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/

i=0;     i=0;
for (j = 0; j <= ny-1; j++)     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;         derx[j * nx + i]      = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
}     }

i=nx-1;     i=nx-1;
for (j = 0; j <= ny-1; j++)     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;         derx[j * nx + i]      = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
}     }

j=0;     j=0;
for (i = 0; i <= nx-1; i++)     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;         dery[j * nx + i]      = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
}     }

j=ny-1;     j=ny-1;
for (i = 0; i <= nx-1; i++)     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;         dery[j * nx + i]      = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );

}     }

 Line 837  int computeHelicity(float *jz_err, float
 Line 833  int computeHelicity(float *jz_err, float
/* Example function 12:  Sum of Absolute Value per polarity  */ /* Example function 12:  Sum of Absolute Value per polarity  */

//  The Sum of the Absolute Value per polarity is defined as the following: //  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 per arcsecond.  //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
//  The units of jz are in G/pix. In this case, we would have the following: //  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), //  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) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
 Line 996  int computeShearAngle(float *bx_err, flo
 Line 992  int computeShearAngle(float *bx_err, flo
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]);             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_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]) );             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]) );

//printf("i,j=%d,%d\n",i,j);
//printf("dotproduct=%f\n",dotproduct);             //printf("dotproduct=%f\n",dotproduct);
//printf("magnitude_potential=%f\n",magnitude_potential);             //printf("magnitude_potential=%f\n",magnitude_potential);
//printf("magnitude_vector=%f\n",magnitude_vector);             //printf("magnitude_vector=%f\n",magnitude_vector);
//printf("dotproduct/(magnitude_potential*magnitude_vector)=%f\n",dotproduct/(magnitude_potential*magnitude_vector));
//printf("acos(dotproduct/(magnitude_potential*magnitude_vector))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector)));
//printf("acos(dotproduct/(magnitude_potential*magnitude_vector)*(180./PI))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI));

shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);             shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
sumsum                  += shear_angle;             sumsum                  += shear_angle;
 Line 1321  void greenpot(float *bx, float *by, floa
 Line 1322  void greenpot(float *bx, float *by, floa
i2e = inx + iwindow;             i2e = inx + iwindow;
if (i2s <   0){i2s =   0;}             if (i2s <   0){i2s =   0;}
if (i2e > nnx){i2e = nnx;}             if (i2e > nnx){i2e = nnx;}
//            for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)              for (j2=j2s;j2<2;j2++){for (i2=i2s;i2<2;i2++)
{             {
float val1 = bztmp[nnx*j2 + 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;                 int   di, dj;
di = abs(i2 - inx);                 di = abs(i2 - inx);
dj = abs(j2 - iny);                 dj = abs(j2 - iny);
sum = sum + val1 * rdist[nnx * dj + di] * dz;                 sum = sum + val1 * rdist[nnx * dj + di] * dz;
//        }                  printf("sum,val1,rdist[nnx * dj + di]=%f,%f,%f\n",sum,val1,rdist[nnx * dj + di]);
} }             } }
pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.             pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
}         }
} } // end of OpenMP parallelism     } } // end of OpenMP parallelism

printf("bztmp[0]=%f\n",bztmp[0]);
printf("bztmp[91746]=%f\n",bztmp[91746]);
printf("pfpot[0]=%f\n",pfpot[0]);
printf("pfpot[91746]=%f\n",pfpot[91746]);

// #pragma omp parallel for private(inx)     // #pragma omp parallel for private(inx)
for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)     for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
{     {
 Line 1348  void greenpot(float *bx, float *by, floa
 Line 1348  void greenpot(float *bx, float *by, floa
by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;         by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
} } // end of OpenMP parallelism     } } // end of OpenMP parallelism

printf("bx[0]=%f\n",bx[0]);
printf("by[0]=%f\n",by[0]);
printf("bx[91746]=%f\n",bx[91746]);
printf("by[91746]=%f\n",by[91746]);
free(rdist);     free(rdist);
free(pfpot);     free(pfpot);
free(bztmp);     free(bztmp);
 Line 1361  char *sw_functions_version() // Returns
 Line 1365  char *sw_functions_version() // Returns
return strdup("\$Id");     return strdup("\$Id");
} }

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

Legend:
 Removed from v.1.33 changed lines Added in v.1.36