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.3 and 1.8

version 1.3, 2012/10/23 18:42:56 version 1.8, 2013/02/09 02:39:20
 Line 1
 Line 1
/*=========================================== /*===========================================

The following 13 functions calculate the following spaceweather indices:     The following 14 functions calculate the following spaceweather indices:

USFLUX Total unsigned flux in Maxwells     USFLUX Total unsigned flux in Maxwells
MEANGAM Mean inclination angle, gamma, in degrees     MEANGAM Mean inclination angle, gamma, in degrees
 Line 81  int computeAbsFlux(float *bz, int *dims,
 Line 81  int computeAbsFlux(float *bz, int *dims,
*absFlux = 0.0;     *absFlux = 0.0;
*mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)                 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 ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[j * nx + i]) continue;
sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
}                 }
 Line 113  int computeBh(float *bx, float *by, floa
 Line 114  int computeBh(float *bx, float *by, floa

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

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)             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] );                 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];                 sum += bh[j * nx + i];
 Line 152  int computeGamma(float *bx, float *by, f
 Line 155  int computeGamma(float *bx, float *by, f
if (bh[j * nx + i] > 100)                 if (bh[j * nx + i] > 100)
{                   {
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;
sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
}                   }
 Line 179  int computeB_total(float *bx, float *by,
 Line 183  int computeB_total(float *bx, float *by,
{           {
for (j = 0; j < ny; j++)             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[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]);
}               }
}           }
 Line 245  int computeBtotalderivative(float *bt, i
 Line 252  int computeBtotalderivative(float *bt, i
}           }

/* Just some print statements
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
printf("j=%d\n",j);
printf("i=%d\n",i);
printf("dery_bt[j*nx+i]=%f\n",dery_bt[j*nx+i]);
printf("derx_bt[j*nx+i]=%f\n",derx_bt[j*nx+i]);
printf("bt[j*nx+i]=%f\n",bt[j*nx+i]);
}
}
*/

for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{             {
// if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;
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(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 */                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 */
}             }
 Line 335  int computeBhderivative(float *bh, int *
 Line 329  int computeBhderivative(float *bh, int *
}           }

/*Just some print statements
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
printf("j=%d\n",j);
printf("i=%d\n",i);
printf("dery_bh[j*nx+i]=%f\n",dery_bh[j*nx+i]);
printf("derx_bh[j*nx+i]=%f\n",derx_bh[j*nx+i]);
printf("bh[j*nx+i]=%f\n",bh[j*nx+i]);
}
}
*/

for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{             {
// if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;
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(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 */                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 */
}             }
 Line 383  int computeBzderivative(float *bz, int *
 Line 364  int computeBzderivative(float *bz, int *
{           {
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;
}               }
}           }
 Line 392  int computeBzderivative(float *bz, int *
 Line 374  int computeBzderivative(float *bz, int *
{           {
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;
}               }
}           }
 Line 401  int computeBzderivative(float *bz, int *
 Line 384  int computeBzderivative(float *bz, int *
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;
}           }

/*Just some print statements
for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
printf("j=%d\n",j);
printf("i=%d\n",i);
printf("dery_bz[j*nx+i]=%f\n",dery_bz[j*nx+i]);
printf("derx_bz[j*nx+i]=%f\n",derx_bz[j*nx+i]);
printf("bz[j*nx+i]=%f\n",bz[j*nx+i]);
}
}
*/

for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{             {
// if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;                // 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 ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(bz[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 */                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 */
}             }
 Line 453  int computeBzderivative(float *bz, int *
 Line 429  int computeBzderivative(float *bz, int *
} }

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

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

//  In discretized space like data pixels, //  In discretized space like data pixels,
 Line 470  int computeBzderivative(float *bz, int *
 Line 445  int computeBzderivative(float *bz, int *
// //
//  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), //  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: //  one must perform the following unit conversions:
//  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or  //  (Gauss)(1/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),  //  (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. //  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). //  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: //  In that case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or //  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
//  jz * (35.0) //  jz * (35.0)
// //
//  The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following: //  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)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(1000.)  //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)
//  =(Gauss/pix)(1/CDELT1)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)  //  = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
//  =(Gauss/pix)(1/0.5)(10^-4)(4*PI*10^7)(722500)(1000.)
//  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)

int computeJz(float *bx, float *by, int *dims, float *jz, int computeJz(float *bx, float *by, int *dims, float *jz,
float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)

 Line 496  int computeJz(float *bx, float *by, int
 Line 471  int computeJz(float *bx, float *by, int

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

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

 Line 506  int computeJz(float *bx, float *by, int
 Line 479  int computeJz(float *bx, float *by, int
{           {
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;
}               }
}           }
 Line 515  int computeJz(float *bx, float *by, int
 Line 489  int computeJz(float *bx, float *by, int
{           {
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;
}               }
}           }
 Line 524  int computeJz(float *bx, float *by, int
 Line 499  int computeJz(float *bx, float *by, int
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;
}           }

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;
}           }

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;
}           }

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;
}           }

/* Just some print statements
for (i = 0; i < nx; i++)          for (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j < ny; j++)              for (j = 0; j <= ny-1; j++)
{               {
printf("j=%d\n",j);                 /* calculate jz at all points */
printf("i=%d\n",i);                 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                              /* jz is in units of Gauss/pix */
printf("derx[j*nx+i]=%f\n",derx[j*nx+i]);              }
printf("bx[j*nx+i]=%f\n",bx[j*nx+i]);
printf("by[j*nx+i]=%f\n",by[j*nx+i]);
}               }

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_smooth,
float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)

{

int nx = dims[0], ny = dims[1];

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

*mean_jz_ptr = 0.0;
float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=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 (i = 0; i <= nx-1; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{             {
// if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                 if isnan(derx[j * nx + i]) continue;
us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */                 if isnan(dery[j * nx + i]) continue;
jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                 if isnan(jz_smooth[j * nx + i]) continue;
//printf("%d,%d,%f\n",i,j,jz_smooth[j * nx + i]);
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 */
}             }
}           }

/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
printf("mean_curl=%f\n",mean_curl);         printf("mean_curl=%f\n",mean_curl);
printf("cdelt1, what is it?=%f\n",cdelt1);         printf("cdelt1, what is it?=%f\n",cdelt1);
*mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */
*us_i_ptr        = (us_i);                   /* us_i gets populated as MEANJZD */          *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
return 0;         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  /* 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. // 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 // 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) //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
//                               = 1/Mm //                               = 1/Mm

int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask, int *bitmask,  int computeAlpha(float *bz, int *dims, float *jz_smooth, float *mean_alpha_ptr, int *mask, int *bitmask,
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], ny = dims[1];
int i, j, count_mask=0;          int i, j=0;

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

*mean_alpha_ptr = 0.0;         *mean_alpha_ptr = 0.0;
float aa, bb, cc, bznew, alpha2, sum=0.0;          float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=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[j * nx + i]) continue;                  if isnan(jz_smooth[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
sum += (jz[j * nx + i] / bz[j * nx + i])*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */                  if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i]);
count_mask++;                  if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i]);
if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);
if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
sum5 += bz[j * nx + i];
}               }
}           }

printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);          sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
printf("sum=%f\n",sum);          /* Determine the sign of alpha */
*mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */          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 */
return 0;         return 0;
} }

/*===========================================*/ /*===========================================*/
/* Example function 10:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */  /* Example function 11:  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 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. //  The units of Jz are in G/pix; the units of Bz are in G.
 Line 637  int computeAlpha(float *bz, int *dims, f
 Line 667  int computeAlpha(float *bz, int *dims, f
//                                                      = G^2 / m. //                                                      = G^2 / m.

int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,  int computeHelicity(float *bz, int *dims, float *jz_smooth, float *mean_ih_ptr, float *total_us_ih_ptr,
float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                                         float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{ {
 Line 650  int computeHelicity(float *bz, int *dims
 Line 680  int computeHelicity(float *bz, int *dims
*mean_ih_ptr = 0.0;         *mean_ih_ptr = 0.0;
float sum=0.0, sum2=0.0;         float sum=0.0, sum2=0.0;

for (j = 0; j < ny; j++)
{
for (i = 0; i < nx; i++)                 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 ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(jz[j * nx + i]) continue;                  if isnan(jz_smooth[j * nx + i]) continue;
if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
if (jz[j * nx + i] == 0.0) continue;                  if (jz_smooth[j * nx + i] == 0.0) continue;
sum  +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum  +=     (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
}                 }
}          }
 Line 675  int computeHelicity(float *bz, int *dims
 Line 705  int computeHelicity(float *bz, int *dims
} }

/*===========================================*/ /*===========================================*/
/* Example function 11:  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. //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
 Line 683  int computeHelicity(float *bz, int *dims
 Line 713  int computeHelicity(float *bz, int *dims
//  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)

int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,  int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr,

{ {
 Line 700  int computeSumAbsPerPolarity(float *bz,
 Line 730  int computeSumAbsPerPolarity(float *bz,
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 (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                  if isnan(bz[j * nx + i]) continue;
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) sum1 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
if (bz[j * nx + i] <= 0) sum2 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
}               }
}           }

 Line 710  int computeSumAbsPerPolarity(float *bz,
 Line 741  int computeSumAbsPerPolarity(float *bz,
} }

/*===========================================*/ /*===========================================*/
/* Example function 12:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */  /* 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 // 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. // automatically yields erg per cubic centimeter for an input B in Gauss.
// //
 Line 740  int computeFreeEnergy(float *bx, float *
 Line 771  int computeFreeEnergy(float *bx, float *
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(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) 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);                  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);
}               }
 Line 751  int computeFreeEnergy(float *bx, float *
 Line 784  int computeFreeEnergy(float *bx, float *
} }

/*===========================================*/ /*===========================================*/
/* 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 */  /* 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, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, 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_angleptr, float *area_w_shear_gt_45ptr,
 Line 777  int computeShearAngle(float *bx, float *
 Line 810  int computeShearAngle(float *bx, float *
if isnan(bpy[j * nx + i]) continue;                  if isnan(bpy[j * nx + i]) continue;
if isnan(bpz[j * nx + i]) continue;                  if isnan(bpz[j * nx + i]) continue;
if isnan(bz[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*/                  /* 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]);                  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]));

Legend:
 Removed from v.1.3 changed lines Added in v.1.8