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.28 and 1.39

version 1.28, 2014/03/13 18:40:43 version 1.39, 2021/03/18 01:53:40
 Line 1
 Line 1

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

The following 14 functions calculate the following spaceweather indices:  The following 14 functions calculate the following spaceweather indices:
 Line 17
 Line 18
MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter  MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
TOTPOT Total photospheric 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  MEANSHR Mean shear angle (measured using Btotal) in degrees
R_VALUE Karel Schrijver's R parameter

The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and  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  pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
 Line 245  int computeB_total(float *bx_err, float
 Line 247  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_termAt, float *err_termBt)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 265  int computeBtotalderivative(float *bt, i
 Line 267  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_termAt[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 274  int computeBtotalderivative(float *bt, i
 Line 277  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_termBt[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 303  int computeBtotalderivative(float *bt, i
 Line 308  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 319  int computeBtotalderivative(float *bt, i
 Line 324  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_termBt[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_termAt[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 337  int computeBtotalderivative(float *bt, i
 Line 342  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_termAh, float *err_termBh)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 357  int computeBhderivative(float *bh, float
 Line 362  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_termAh[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 366  int computeBhderivative(float *bh, float
 Line 372  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_termBh[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 411  int computeBhderivative(float *bh, float
 Line 419  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_termBh[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_termAh[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 428  int computeBhderivative(float *bh, float
 Line 436  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_termA, float *err_termB)
{ {

int nx = dims[0];     int nx = dims[0];
 Line 447  int computeBzderivative(float *bz, float
 Line 455  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_termA[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 457  int computeBzderivative(float *bz, float
 Line 465  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_termB[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 508  int computeBzderivative(float *bz, float
 Line 513  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_termB[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_termA[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 562  int computeBzderivative(float *bz, float
 Line 565  int computeBzderivative(float *bz, float
//              float *noiseby, float *noisebz) //              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 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 *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)

{ {
 Line 581  int computeJz(float *bx_err, float *by_e
 Line 584  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]);
}         }
}     }

 Line 590  int computeJz(float *bx_err, float *by_e
 Line 593  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]);
}         }
}     }

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

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

for (i = 1; i <= nx-2; i++)
for (i = 0; i <= nx-1; i++)
{     {
for (j = 1; j <= ny-2; j++)          for (j = 0; j <= ny-1; 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( err_term1[j * nx + i] + err_term2[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)]) ) ;
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 831  int computeHelicity(float *jz_err, float
 Line 831  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.  //  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 868  int computeSumAbsPerPolarity(float *jz_e
 Line 868  int computeSumAbsPerPolarity(float *jz_e
}         }
}     }

*totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */      *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are Amperes per arcsecond */
*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);
 Line 1022  int computeShearAngle(float *bx_err, flo
 Line 1022  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 = (sumsum)/(count);                 /* Units are degrees */     *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */

// For the error in the mean 3D shear angle:
// If count_mask is 0, then we run into a divide by zero error. In this case, set *meanshear_angle_err_ptr to NAN
// If count_mask is greater than zero, then compute the error.
if (count_mask == 0)
*meanshear_angle_err_ptr = NAN;
else

/* The area here is a fractional area -- the % of the total area. This has no error associated with it. */     /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */

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

return 0;         return 0;
} }

/*===========================================*/ /*===========================================*/
/* Example function 15: R parameter as defined in Schrijver, 2007 */
//
// Note that there is a restriction on the function fsample()
// If the following occurs:
//      nx_out > floor((ny_in-1)/scale + 1)
//      ny_out > floor((ny_in-1)/scale + 1),
// where n*_out are the dimensions of the output array and n*_in
// are the dimensions of the input array, fsample() will usually result
// in a segfault (though not always, depending on how the segfault was accessed.)

int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1, int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,              float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
float *pmap, int nx1, int ny1)               float *pmap, int nx1, int ny1,
{               int scale, float *p1pad, int nxp, int nyp, float *pmapn)

{
int nx = dims[0];     int nx = dims[0];
int ny = dims[1];     int ny = dims[1];
int i = 0;     int i = 0;
int j = 0;     int j = 0;
int index;      int index, index1;
double sum = 0.0;     double sum = 0.0;
double err = 0.0;     double err = 0.0;
*Rparam = 0.0;     *Rparam = 0.0;
struct fresize_struct fresboxcar, fresgauss;     struct fresize_struct fresboxcar, fresgauss;
int scale = round(2.0/cdelt1);      struct fint_struct fints;
float sigma = 10.0/2.3548;     float sigma = 10.0/2.3548;

// set up convolution kernels
init_fresize_boxcar(&fresboxcar,1,1);     init_fresize_boxcar(&fresboxcar,1,1);

// set up convolution kernel
init_fresize_gaussian(&fresgauss,sigma,20,1);     init_fresize_gaussian(&fresgauss,sigma,20,1);

// =============== [STEP 1] ===============
// bin the line-of-sight magnetogram down by a factor of scale
fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);     fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);

// =============== [STEP 2] ===============
// identify positive and negative pixels greater than +/- 150 gauss
// and label those pixels with a 1.0 in arrays p1p0 and p1n0
for (i = 0; i < nx1; i++)     for (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
 Line 1074  int computeR(float *bz_err, float *los,
 Line 1096  int computeR(float *bz_err, float *los,
}         }
}     }

// =============== [STEP 3] ===============
// smooth each of the negative and positive pixel bitmaps
fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);     fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);     fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);

// =============== [STEP 4] ===============
// find the pixels for which p1p and p1n are both equal to 1.
// this defines the polarity inversion line
for (i = 0; i < nx1; i++)     for (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
{         {
index = j * nx1 + i;             index = j * nx1 + i;
if (p1p[index] > 0 && p1n[index] > 0)              if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
p1[index]=1.0;                 p1[index]=1.0;
else             else
p1[index]=0.0;                 p1[index]=0.0;
}         }
}     }

fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);      // pad p1 with zeroes so that the gaussian colvolution in step 5
// does not cut off data within hwidth of the edge

// step i: zero p1pad
for (i = 0; i < nxp; i++)
{
for (j = 0; j < nyp; j++)
{
index = j * nxp + i;
}
}

// step ii: place p1 at the center of p1pad
for (i = 0; i < nx1; i++)
{
for (j = 0; j < ny1; j++)
{
index  = j * nx1 + i;
index1 = (j+20) * nxp + (i+20);
}
}

// =============== [STEP 5] ===============
// convolve the polarity inversion line map with a gaussian
// to identify the region near the plarity inversion line
// the resultant array is called pmap
fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);

// select out the nx1 x ny1 non-padded array  within pmap
for (i = 0; i < nx1; i++)
{
for (j = 0; j < ny1; j++)
{
index  = j * nx1 + i;
index1 = (j+20) * nxp + (i+20);
pmapn[index]=pmap[index1];
}
}

// =============== [STEP 6] ===============
// the R parameter is calculated
for (i = 0; i < nx1; i++)     for (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
{         {
index = j * nx1 + i;             index = j * nx1 + i;
sum += pmap[index]*abs(rim[index]);              if isnan(pmapn[index]) continue;
if isnan(rim[index]) continue;
sum += pmapn[index]*abs(rim[index]);
}         }
}     }

 Line 1105  int computeR(float *bz_err, float *los,
 Line 1176  int computeR(float *bz_err, float *los,
else     else
*Rparam = log10(sum);         *Rparam = log10(sum);

//printf("R_VALUE=%f\n",*Rparam);

free_fresize(&fresboxcar);     free_fresize(&fresboxcar);
free_fresize(&fresgauss);     free_fresize(&fresgauss);

return 0;     return 0;

}

/*===========================================*/
/* Example function 16: Lorentz force as defined in Fisher, 2012 */
//
// This calculation is adapted from Xudong's code
// at /proj/cgem/lorentz/apps/lorentz.c

int computeLorentz(float *bx,  float *by, float *bz, float *fx, float *fy, float *fz, int *dims,
float *totfx_ptr, float *totfy_ptr, float *totfz_ptr, float *totbsq_ptr,
float *epsx_ptr, float *epsy_ptr, float *epsz_ptr, int *mask, int *bitmask,
float cdelt1, double rsun_ref, double rsun_obs)

{

int nx = dims[0];
int ny = dims[1];
int nxny = nx*ny;
int j = 0;
int index;
double totfx = 0, totfy = 0, totfz = 0;
double bsq = 0, totbsq = 0;
double epsx = 0, epsy = 0, epsz = 0;
double area = cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
double k_h = -1.0 * area / (4. * PI) / 1.0e20;
double k_z = area / (8. * PI) / 1.0e20;

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

for (int i = 0; i < nxny; i++)
{
if ( mask[i] < 70 || bitmask[i] < 30 ) continue;
if isnan(bx[i]) continue;
if isnan(by[i]) continue;
if isnan(bz[i]) continue;
fx[i]  = bx[i] * bz[i] * k_h;
fy[i]  = by[i] * bz[i] * k_h;
fz[i]  = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];
totbsq += bsq;
}

*totfx_ptr  = totfx;
*totfy_ptr  = totfy;
*totfz_ptr  = totfz;
*totbsq_ptr = totbsq;
*epsx_ptr   = (totfx / k_h) / totbsq;
*epsy_ptr   = (totfy / k_h) / totbsq;
*epsz_ptr   = (totfz / k_z) / totbsq;

//printf("TOTBSQ=%f\n",*totbsq_ptr);

return 0;

}

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

/* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field */

//  To compute the unsigned flux, we simply calculate
//  flux = surface integral [(vector LOS) dot (normal vector)],
//       = surface integral [(magnitude LOS)*(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*cm^2

int computeAbsFlux_los(float *los, int *dims, float *absFlux,
float *mean_vf_ptr, float *count_mask_ptr,
int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{

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

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

for (i = 0; i < nx; i++)
{
for (j = 0; j < ny; j++)
{
if ( bitmask[j * nx + i] < 30 ) continue;
if isnan(los[j * nx + i]) continue;
sum += (fabs(los[j * nx + i]));
}
}

*mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;

return 0;
}

/*===========================================*/
/* Example function 18:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */

int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
{

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

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

/* 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_los[j * nx + i] = (los[j * nx + i+1] - los[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_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
}
}

/* 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;
for (j = 0; j <= ny-1; j++)
{
derx_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
}

i=nx-1;
for (j = 0; j <= ny-1; j++)
{
derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
}

j=0;
for (i = 0; i <= nx-1; i++)
{
dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5;
} }

j=ny-1;
for (i = 0; i <= nx-1; i++)
{
dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
}

for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{
if ( bitmask[j * nx + i] < 30 ) continue;
if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue;
if isnan(los[j * nx + i])      continue;
if isnan(los[(j+1) * nx + i])  continue;
if isnan(los[(j-1) * nx + i])  continue;
if isnan(los[j * nx + i-1])    continue;
if isnan(los[j * nx + i+1])    continue;
if isnan(derx_los[j * nx + i]) continue;
if isnan(dery_los[j * nx + i]) continue;
sum += sqrt( derx_los[j * nx + i]*derx_los[j * nx + i]  + dery_los[j * nx + i]*dery_los[j * nx + i] ); /* Units of Gauss */
}
}

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

return 0;
}

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

Legend:
 Removed from v.1.28 changed lines Added in v.1.39