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.9 and 1.38

version 1.9, 2013/05/20 23:55:32 version 1.38, 2020/06/30 22:38:17
 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
coordinate bitmaps are interpolated.   coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
contain nearest-neighbor bitmaps).

In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig    In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:   and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:

For conf_disambig:    For conf_disambig:
50 : not all solutions agree (weak field method applied)    50 : not all solutions agree (weak field method applied)
 Line 39
 Line 43

Written by Monica Bobra 15 August 2012    Written by Monica Bobra 15 August 2012
Potential Field code (appended) written by Keiji Hayashi    Potential Field code (appended) written by Keiji Hayashi
Error analysis modification 21 October 2013

===========================================*/ ===========================================*/
#include <math.h> #include <math.h>
 Line 60
 Line 65
//  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel. //  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). //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2 //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
//  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2  //  =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,
 Line 73  int computeAbsFlux(float *bz_err, float
 Line 73  int computeAbsFlux(float *bz_err, float

{ {

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 97  int computeAbsFlux(float *bz_err, float
 Line 101  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("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 113  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 124  int computeBh(float *bx_err, float *by_e
 Line 127  int computeBh(float *bx_err, float *by_e
{           {
for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
{               {
if isnan(bx[j * nx + i]) continue;              if isnan(bx[j * nx + i])
if isnan(by[j * nx + i]) continue;              {
bh[j * nx + i] = NAN;
bh_err[j * nx + i] = NAN;
continue;
}
if isnan(by[j * nx + i])
{
bh[j * nx + i] = NAN;
bh_err[j * nx + i] = NAN;
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];
bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];                 bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];
 Line 141  int computeBh(float *bx_err, float *by_e
 Line 154  int computeBh(float *bx_err, float *by_e
/*===========================================*/ /*===========================================*/
/* Example function 3: Calculate Gamma in units of degrees */ /* Example function 3: Calculate Gamma in units of degrees */
// Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI) // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
// Redo calculation in radians for error analysis (since derivatives are only true in units of radians).  //
// Error analysis calculations are done in radians (since derivatives are only true in units of radians),
// and multiplied by (180./PI) at the end for consistency in units.

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 165  int computeGamma(float *bz_err, float *b
 Line 182  int computeGamma(float *bz_err, float *b
if isnan(bz[j * nx + i]) continue;                     if isnan(bz[j * nx + i]) continue;
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 isnan(bh[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 += (1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))))*(1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])))) *
( ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) +
((bh[j * nx + i]*bh[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i])) );
}                   }
}               }
}           }

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 207  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 196  int computeB_total(float *bx_err, float
 Line 219  int computeB_total(float *bx_err, float
{           {
for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
{               {
if isnan(bx[j * nx + i]) continue;              if isnan(bx[j * nx + i])
if isnan(by[j * nx + i]) continue;              {
if isnan(bz[j * nx + i]) continue;                  bt[j * nx + i] = NAN;
bt_err[j * nx + i] = NAN;
continue;
}
if isnan(by[j * nx + i])
{
bt[j * nx + i] = NAN;
bt_err[j * nx + i] = NAN;
continue;
}
if isnan(bz[j * nx + i])
{
bt[j * nx + i] = NAN;
bt_err[j * nx + i] = NAN;
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]);
bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];                 bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];
}               }
 Line 209  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], 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 = 3; i <= nx-4; i++)      for (i = 1; i <= nx-2; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{               {
derx_bt[j * nx + i] = (-bt[j * nx + (i-3)] + 9.0*bt[j * nx + (i-2)] - 45.0*bt[j * nx + (i-1)] + 45*bt[j * nx + (i+1)] - 9.0*bt[j * nx + (i+2)] + bt[j * nx + (i+3)])*(1.0/60.0);             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)])) ;
}               }
}           }

/* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 3; j <= ny-4; j++)          for (j = 1; j <= ny-2; j++)
{               {
dery_bt[j * nx + i] = (-bt[(j-3) * nx + i] + 9.0*bt[(j-2) * nx + i] - 45.0*bt[(j-1) * nx + i] + 45*bt[(j+1) * nx + i] - 9.0*bt[(j+2) * nx + i] + bt[(j+3) * nx + i])*(1.0/60.0);             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: 3 pixels on each edge, for a total of 12 edge formulae 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++)
{           {
derx_bt[j * nx + i] = (-147.0*bt[j * nx + i] + 360.0*bt[j * nx + (i+1)] - 450.0*bt[j * nx + (i+2)] + 400.0*bt[j * nx + (i+3)] - 225.0*bt[j * nx + (i+4)] + 72.0*bt[j * nx + (i+5)] - 10.0*bt[j * nx + (i+6)])*(1.0/60.0);          derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
}           }

i=nx-1;         i=nx-1;
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{           {
derx_bt[j * nx + i] = ( 147.0*bt[j * nx + i] - 360.0*bt[j * nx + (i-1)] + 450.0*bt[j * nx + (i-2)] - 400.0*bt[j * nx + (i-3)] + 225.0*bt[j * nx + (i-4)] - 72.0*bt[j * nx + (i-5)] + 10.0*bt[j * nx + (i-6)])*(1.0/60.0);          derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
}

i=1;
for (j = 0; j <= ny-2; j++)
{
derx_bt[j * nx + i] = (-10.0*bt[j * nx + i] - 77.0*bt[j * nx + (i+1)] + 150.0*bt[j * nx + (i+2)] - 100.0*bt[j * nx + (i+3)] + 50.0*bt[j * nx + (i+4)] - 15.0*bt[j * nx + (i+5)] + 2.0*bt[j * nx + (i+6)])*(1.0/60.0);
}           }

i=nx-2;
for (j = 0; j <= ny-2; j++)
{
derx_bt[j * nx + i] = ( 10.0*bt[j * nx + i] + 77.0*bt[j * nx + (i-1)] - 150.0*bt[j * nx + (i-2)] + 100.0*bt[j * nx + (i-3)] - 50.0*bt[j * nx + (i-4)] + 15.0*bt[j * nx + (i-5)] - 2.0*bt[j * nx + (i-6)])*(1.0/60.0);
}

i=2;
for (j = 0; j <= ny-2; j++)
{
derx_bt[j * nx + i] = ( 2.0*bt[j * nx + i] - 24.0*bt[j * nx + (i+1)] - 35.0*bt[j * nx + (i+2)] + 80.0*bt[j * nx + (i+3)] - 30.0*bt[j * nx + (i+4)] + 8.0*bt[j * nx + (i+5)] - bt[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-3;
for (j = 0; j <= ny-2; j++)
{
derx_bt[j * nx + i] = (-2.0*bt[j * nx + i] + 24.0*bt[j * nx + (i-1)] + 35.0*bt[j * nx + (i-2)] - 80.0*bt[j * nx + (i-3)] + 30.0*bt[j * nx + (i-4)] - 8.0*bt[j * nx + (i-5)] + bt[j * nx + (i-6)])*(1.0/60.0);
}

j=0;         j=0;
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
dery_bt[j * nx + i] = (-147.0*bt[j * nx + i] + 360.0*bt[(j+1) * nx + i] - 450.0*bt[(j+2) * nx + i] + 400.0*bt[(j+3) * nx + i] - 225.0*bt[(j+4) * nx + i] + 72.0*bt[(j+5) * nx + i] - 10.0*bt[(j+6) * nx + i])*(1.0/60.0);          dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5;
}           }

j=ny-1;         j=ny-1;
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
dery_bt[j * nx + i] = ( 147.0*bt[j * nx + i] - 360.0*bt[(j-1) * nx + i] + 450.0*bt[(j-2) * nx + i] - 400.0*bt[(j-3) * nx + i] + 225.0*bt[(j-4) * nx + i] - 72.0*bt[(j-5) * nx + i] + 10.0*bt[(j-6) * nx + i])*(1.0/60.0);          dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
}           }

j=1;      // Calculate the sum only
for (i = 0; i <= nx-2; i++)      for (i = 1; i <= nx-2; i++)
{           {
dery_bt[j * nx + i] = (-10.0*bt[j * nx + i] - 77.0*bt[(j+1) * nx + i] + 150.0*bt[(j+2) * nx + i] - 100.0*bt[(j+3) * nx + i] + 50.0*bt[(j+4) * nx + i] - 15.0*bt[(j+5) * nx + i] + 2.0*bt[(j+6) * nx + i])*(1.0/60.0);          for (j = 1; j <= ny-2; j++)
}

j=ny-2;
for (i = 0; i <= nx-2; i++)
{
dery_bt[j * nx + i] = ( 10.0*bt[j * nx + i] + 77.0*bt[(j-1) * nx + i] - 150.0*bt[(j-2) * nx + i] + 100.0*bt[(j-3) * nx + i] - 50.0*bt[(j-4) * nx + i] + 15.0*bt[(j-5) * nx + i] - 2.0*bt[(j-6) * nx + i])*(1.0/60.0);
}

j=2;
for (i = 0; i <= nx-3; i++)
{
dery_bt[j * nx + i] = ( 2.0*bt[j * nx + i] - 24.0*bt[(j+1) * nx + i] - 35.0*bt[(j+2) * nx + i] + 80.0*bt[(j+3) * nx + i] - 30.0*bt[(j+4) * nx + i] + 8.0*bt[(j+5) * nx + i] - bt[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-3;
for (i = 0; i <= nx-3; i++)
{
dery_bt[j * nx + i] = (-2.0*bt[j * nx + i] + 24.0*bt[(j-1) * nx + i] + 35.0*bt[(j-2) * nx + i] - 80.0*bt[(j-3) * nx + i] + 30.0*bt[(j-4) * nx + i] - 8.0*bt[(j-5) * nx + i] + bt[(j-6) * nx + i])*(1.0/60.0);
}

for (i = 0; i <= nx-1; i++)
{
for (j = 0; j <= ny-1; j++)
{             {
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
if isnan(bt[j * nx + i])      continue;
if isnan(bt[(j+1) * nx + i])  continue;
if isnan(bt[(j-1) * nx + i])  continue;
if isnan(bt[j * nx + i-1])    continue;
if isnan(bt[j * nx + i+1])    continue;
if isnan(bt_err[j * nx + i])  continue;
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 += (2.0)*bt_err[j * nx + i]*bt_err[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]  ))+
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]  )) ;
}             }
}           }

*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);
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 341  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], 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 = 3; i <= nx-4; i++)      for (i = 1; i <= nx-2; i++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{               {
derx_bh[j * nx + i] = (-bh[j * nx + (i-3)] + 9.0*bh[j * nx + (i-2)] - 45.0*bh[j * nx + (i-1)] + 45*bh[j * nx + (i+1)] - 9.0*bh[j * nx + (i+2)] + bh[j * nx + (i+3)])*(1.0/60.0);             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)]));
}               }
}           }

/* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 3; j <= ny-4; j++)         for (j = 1; j <= ny-2; j++)
{               {
dery_bh[j * nx + i] = (-bh[(j-3) * nx + i] + 9.0*bh[(j-2) * nx + i] - 45.0*bh[(j-1) * nx + i] + 45*bh[(j+1) * nx + i] - 9.0*bh[(j+2) * nx + i] + bh[(j+3) * nx + i])*(1.0/60.0);            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: 3 pixels on each edge, for a total of 12 edge formulae 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++)
{           {
derx_bh[j * nx + i] = (-147.0*bh[j * nx + i] + 360.0*bh[j * nx + (i+1)] - 450.0*bh[j * nx + (i+2)] + 400.0*bh[j * nx + (i+3)] - 225.0*bh[j * nx + (i+4)] + 72.0*bh[j * nx + (i+5)] - 10.0*bh[j * nx + (i+6)])*(1.0/60.0);          derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
}           }

i=nx-1;         i=nx-1;
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{           {
derx_bh[j * nx + i] = ( 147.0*bh[j * nx + i] - 360.0*bh[j * nx + (i-1)] + 450.0*bh[j * nx + (i-2)] - 400.0*bh[j * nx + (i-3)] + 225.0*bh[j * nx + (i-4)] - 72.0*bh[j * nx + (i-5)] + 10.0*bh[j * nx + (i-6)])*(1.0/60.0);          derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
}

i=1;
for (j = 0; j <= ny-2; j++)
{
derx_bh[j * nx + i] = (-10.0*bh[j * nx + i] - 77.0*bh[j * nx + (i+1)] + 150.0*bh[j * nx + (i+2)] - 100.0*bh[j * nx + (i+3)] + 50.0*bh[j * nx + (i+4)] - 15.0*bh[j * nx + (i+5)] + 2.0*bh[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-2;
for (j = 0; j <= ny-2; j++)
{
derx_bh[j * nx + i] = ( 10.0*bh[j * nx + i] + 77.0*bh[j * nx + (i-1)] - 150.0*bh[j * nx + (i-2)] + 100.0*bh[j * nx + (i-3)] - 50.0*bh[j * nx + (i-4)] + 15.0*bh[j * nx + (i-5)] - 2.0*bh[j * nx + (i-6)])*(1.0/60.0);
}           }

i=2;
for (j = 0; j <= ny-2; j++)
{
derx_bh[j * nx + i] = ( 2.0*bh[j * nx + i] - 24.0*bh[j * nx + (i+1)] - 35.0*bh[j * nx + (i+2)] + 80.0*bh[j * nx + (i+3)] - 30.0*bh[j * nx + (i+4)] + 8.0*bh[j * nx + (i+5)] - bh[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-3;
for (j = 0; j <= ny-2; j++)
{
derx_bh[j * nx + i] = (-2.0*bh[j * nx + i] + 24.0*bh[j * nx + (i-1)] + 35.0*bh[j * nx + (i-2)] - 80.0*bh[j * nx + (i-3)] + 30.0*bh[j * nx + (i-4)] - 8.0*bh[j * nx + (i-5)] + bh[j * nx + (i-6)])*(1.0/60.0);
}

j=0;         j=0;
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
dery_bh[j * nx + i] = (-147.0*bh[j * nx + i] + 360.0*bh[(j+1) * nx + i] - 450.0*bh[(j+2) * nx + i] + 400.0*bh[(j+3) * nx + i] - 225.0*bh[(j+4) * nx + i] + 72.0*bh[(j+5) * nx + i] - 10.0*bh[(j+6) * nx + i])*(1.0/60.0);          dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5;
}           }

j=ny-1;         j=ny-1;
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
dery_bh[j * nx + i] = ( 147.0*bh[j * nx + i] - 360.0*bh[(j-1) * nx + i] + 450.0*bh[(j-2) * nx + i] - 400.0*bh[(j-3) * nx + i] + 225.0*bh[(j-4) * nx + i] - 72.0*bh[(j-5) * nx + i] + 10.0*bh[(j-6) * nx + i])*(1.0/60.0);          dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
}

j=1;
for (i = 0; i <= nx-2; i++)
{
dery_bh[j * nx + i] = (-10.0*bh[j * nx + i] - 77.0*bh[(j+1) * nx + i] + 150.0*bh[(j+2) * nx + i] - 100.0*bh[(j+3) * nx + i] + 50.0*bh[(j+4) * nx + i] - 15.0*bh[(j+5) * nx + i] + 2.0*bh[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-2;
for (i = 0; i <= nx-2; i++)
{
dery_bh[j * nx + i] = ( 10.0*bh[j * nx + i] + 77.0*bh[(j-1) * nx + i] - 150.0*bh[(j-2) * nx + i] + 100.0*bh[(j-3) * nx + i] - 50.0*bh[(j-4) * nx + i] + 15.0*bh[(j-5) * nx + i] - 2.0*bh[(j-6) * nx + i])*(1.0/60.0);
}

j=2;
for (i = 0; i <= nx-3; i++)
{
dery_bh[j * nx + i] = ( 2.0*bh[j * nx + i] - 24.0*bh[(j+1) * nx + i] - 35.0*bh[(j+2) * nx + i] + 80.0*bh[(j+3) * nx + i] - 30.0*bh[(j+4) * nx + i] + 8.0*bh[(j+5) * nx + i] - bh[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-3;
for (i = 0; i <= nx-3; i++)
{
dery_bh[j * nx + i] = (-2.0*bh[j * nx + i] + 24.0*bh[(j-1) * nx + i] + 35.0*bh[(j-2) * nx + i] - 80.0*bh[(j-3) * nx + i] + 30.0*bh[(j-4) * nx + i] - 8.0*bh[(j-5) * nx + i] + bh[(j-6) * nx + i])*(1.0/60.0);
}           }

 Line 453  int computeBhderivative(float *bh, float
 Line 409  int computeBhderivative(float *bh, float
for (j = 0; j <= ny-1; j++)             for (j = 0; 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 ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
if isnan(bh[j * nx + i])      continue;
if isnan(bh[(j+1) * nx + i])  continue;
if isnan(bh[(j-1) * nx + i])  continue;
if isnan(bh[j * nx + i-1])    continue;
if isnan(bh[j * nx + i+1])    continue;
if isnan(bh_err[j * nx + i])  continue;
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 += (2.0)*bh_err[j * nx + i]*bh_err[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]  ))+
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]  )) ;
}             }
}           }

*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 472  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], 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_bz_ptr = 0.0;         *mean_derivative_bz_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 = 3; i <= nx-4; i++)      for (i = 1; i <= nx-2; i++)
{           {
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-3)] + 9.0*bz[j * nx + (i-2)] - 45.0*bz[j * nx + (i-1)] + 45*bz[j * nx + (i+1)] - 9.0*bz[j * nx + (i+2)] + bz[j * nx + (i+3)])*(1.0/60.0);             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)]));
}               }
}           }

/* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 3; j <= ny-4; 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-3) * nx + i] + 9.0*bz[(j-2) * nx + i] - 45.0*bz[(j-1) * nx + i] + 45*bz[(j+1) * nx + i] - 9.0*bz[(j+2) * nx + i] + bz[(j+3) * nx + i])*(1.0/60.0);             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: 3 pixels on each edge, for a total of 12 edge formulae 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(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] = (-147.0*bz[j * nx + i] + 360.0*bz[j * nx + (i+1)] - 450.0*bz[j * nx + (i+2)] + 400.0*bz[j * nx + (i+3)] - 225.0*bz[j * nx + (i+4)] + 72.0*bz[j * nx + (i+5)] - 10.0*bz[j * nx + (i+6)])*(1.0/60.0);
}           }

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] = ( 147.0*bz[j * nx + i] - 360.0*bz[j * nx + (i-1)] + 450.0*bz[j * nx + (i-2)] - 400.0*bz[j * nx + (i-3)] + 225.0*bz[j * nx + (i-4)] - 72.0*bz[j * nx + (i-5)] + 10.0*bz[j * nx + (i-6)])*(1.0/60.0);
}           }

i=1;
for (j = 0; j <= ny-2; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = (-10.0*bz[j * nx + i] - 77.0*bz[j * nx + (i+1)] + 150.0*bz[j * nx + (i+2)] - 100.0*bz[j * nx + (i+3)] + 50.0*bz[j * nx + (i+4)] - 15.0*bz[j * nx + (i+5)] + 2.0*bz[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-2;
for (j = 0; j <= ny-2; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = ( 10.0*bz[j * nx + i] + 77.0*bz[j * nx + (i-1)] - 150.0*bz[j * nx + (i-2)] + 100.0*bz[j * nx + (i-3)] - 50.0*bz[j * nx + (i-4)] + 15.0*bz[j * nx + (i-5)] - 2.0*bz[j * nx + (i-6)])*(1.0/60.0);
}

i=2;
for (j = 0; j <= ny-2; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = ( 2.0*bz[j * nx + i] - 24.0*bz[j * nx + (i+1)] - 35.0*bz[j * nx + (i+2)] + 80.0*bz[j * nx + (i+3)] - 30.0*bz[j * nx + (i+4)] + 8.0*bz[j * nx + (i+5)] - bz[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-3;
for (j = 0; j <= ny-2; j++)
{
if isnan(bz[j * nx + i]) continue;
derx_bz[j * nx + i] = (-2.0*bz[j * nx + i] + 24.0*bz[j * nx + (i-1)] + 35.0*bz[j * nx + (i-2)] - 80.0*bz[j * nx + (i-3)] + 30.0*bz[j * nx + (i-4)] - 8.0*bz[j * nx + (i-5)] + bz[j * nx + (i-6)])*(1.0/60.0);
}

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] = (-147.0*bz[j * nx + i] + 360.0*bz[(j+1) * nx + i] - 450.0*bz[(j+2) * nx + i] + 400.0*bz[(j+3) * nx + i] - 225.0*bz[(j+4) * nx + i] + 72.0*bz[(j+5) * nx + i] - 10.0*bz[(j+6) * nx + i])*(1.0/60.0);
}           }

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] = ( 147.0*bz[j * nx + i] - 360.0*bz[(j-1) * nx + i] + 450.0*bz[(j-2) * nx + i] - 400.0*bz[(j-3) * nx + i] + 225.0*bz[(j-4) * nx + i] - 72.0*bz[(j-5) * nx + i] + 10.0*bz[(j-6) * nx + i])*(1.0/60.0);
}

j=1;
for (i = 0; i <= nx-2; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = (-10.0*bz[j * nx + i] - 77.0*bz[(j+1) * nx + i] + 150.0*bz[(j+2) * nx + i] - 100.0*bz[(j+3) * nx + i] + 50.0*bz[(j+4) * nx + i] - 15.0*bz[(j+5) * nx + i] + 2.0*bz[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-2;
for (i = 0; i <= nx-2; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = ( 10.0*bz[j * nx + i] + 77.0*bz[(j-1) * nx + i] - 150.0*bz[(j-2) * nx + i] + 100.0*bz[(j-3) * nx + i] - 50.0*bz[(j-4) * nx + i] + 15.0*bz[(j-5) * nx + i] - 2.0*bz[(j-6) * nx + i])*(1.0/60.0);
}

j=2;
for (i = 0; i <= nx-3; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = ( 2.0*bz[j * nx + i] - 24.0*bz[(j+1) * nx + i] - 35.0*bz[(j+2) * nx + i] + 80.0*bz[(j+3) * nx + i] - 30.0*bz[(j+4) * nx + i] + 8.0*bz[(j+5) * nx + i] - bz[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-3;
for (i = 0; i <= nx-3; i++)
{
if isnan(bz[j * nx + i]) continue;
dery_bz[j * nx + i] = (-2.0*bz[j * nx + i] + 24.0*bz[(j-1) * nx + i] + 35.0*bz[(j-2) * nx + i] - 80.0*bz[(j-3) * nx + i] + 30.0*bz[(j-4) * nx + i] - 8.0*bz[(j-5) * nx + i] + bz[(j-6) * nx + i])*(1.0/60.0);
}           }

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

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

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

/* 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 627  int computeBzderivative(float *bz, float
 Line 535  int computeBzderivative(float *bz, float
//  the integration of the field Bx and By along //  the integration of the field Bx and By along
//  the circumference of the data pixel divided by the area of the pixel. //  the circumference of the data pixel divided by the area of the pixel.
//  One form of differencing (a word for the differential operator //  One form of differencing (a word for the differential operator
//  in the discretized space) of the curl is expressed as the following,  //  in the discretized space) of the curl is expressed as
//  which utilizes a second-order finite difference method:

//  (dx * (Bx(i,j-1)+Bx(i,j)) / 2 //  (dx * (Bx(i,j-1)+Bx(i,j)) / 2
//  +dy * (By(i+1,j)+By(i,j)) / 2 //  +dy * (By(i+1,j)+By(i,j)) / 2
//  -dx * (Bx(i,j+1)+Bx(i,j)) / 2 //  -dx * (Bx(i,j+1)+Bx(i,j)) / 2
//  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) //  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
// //
// //
//  However, for the purposes of this calculation, we will use a sixth-order finite difference
//  method taken from the pencil code:
//
//  dBy/dx = (    -By*(i-3,j) + 9By(i-2,j) - 45By(i-1,j) + 45By(i+1,j) - 9By(i+2,j) + By(i+3,j)    )/ 60
//  and similarly for dBx/dy.
//
//  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)(1/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
 Line 664  int computeBzderivative(float *bz, float
 Line 564  int computeBzderivative(float *bz, float
//            int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx, //            int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
//              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)

{ {
int nx = dims[0], ny = dims[1];      int nx = dims[0];
int i, j, count_mask=0;      int ny = dims[1];
printf("nx=%d\n",nx);      int i = 0;
printf("ny=%d\n",ny);      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 = 3; i <= nx-4; i++)
for (i = 1; i <= nx-2; i++)
{           {
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-3)] + 9.0*by[j * nx + (i-2)] - 45.0*by[j * nx + (i-1)] + 45*by[j * nx + (i+1)] - 9.0*by[j * nx + (i+2)] + by[j * nx + (i+3)])*(1.0/60.0);             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]);
}               }
}           }

/* brute force method of calculating the derivative (no consideration for edges) */
for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
{           {
for (j = 3; j <= ny-4; 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-3) * nx + i] + 9.0*bx[(j-2) * nx + i] - 45.0*bx[(j-1) * nx + i] + 45*bx[(j+1) * nx + i] - 9.0*bx[(j+2) * nx + i] + bx[(j+3) * nx + i])*(1.0/60.0);             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: 3 pixels on each edge, for a total of 12 edge formulae below */      /* 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] = (-147.0*by[j * nx + i] + 360.0*by[j * nx + (i+1)] - 450.0*by[j * nx + (i+2)] + 400.0*by[j * nx + (i+3)] - 225.0*by[j * nx + (i+4)] + 72.0*by[j * nx + (i+5)] - 10.0*by[j * nx + (i+6)])*(1.0/60.0);
}           }

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] = ( 147.0*by[j * nx + i] - 360.0*by[j * nx + (i-1)] + 450.0*by[j * nx + (i-2)] - 400.0*by[j * nx + (i-3)] + 225.0*by[j * nx + (i-4)] - 72.0*by[j * nx + (i-5)] + 10.0*by[j * nx + (i-6)])*(1.0/60.0);
}

i=1;
for (j = 0; j <= ny-2; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = (-10.0*by[j * nx + i] - 77.0*by[j * nx + (i+1)] + 150.0*by[j * nx + (i+2)] - 100.0*by[j * nx + (i+3)] + 50.0*by[j * nx + (i+4)] - 15.0*by[j * nx + (i+5)] + 2.0*by[j * nx + (i+6)])*(1.0/60.0);
}           }

i=nx-2;
for (j = 0; j <= ny-2; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = ( 10.0*by[j * nx + i] + 77.0*by[j * nx + (i-1)] - 150.0*by[j * nx + (i-2)] + 100.0*by[j * nx + (i-3)] - 50.0*by[j * nx + (i-4)] + 15.0*by[j * nx + (i-5)] - 2.0*by[j * nx + (i-6)])*(1.0/60.0);
}

i=2;
for (j = 0; j <= ny-2; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = ( 2.0*by[j * nx + i] - 24.0*by[j * nx + (i+1)] - 35.0*by[j * nx + (i+2)] + 80.0*by[j * nx + (i+3)] - 30.0*by[j * nx + (i+4)] + 8.0*by[j * nx + (i+5)] - by[j * nx + (i+6)])*(1.0/60.0);
}

i=nx-3;
for (j = 0; j <= ny-2; j++)
{
if isnan(by[j * nx + i]) continue;
derx[j * nx + i] = (-2.0*by[j * nx + i] + 24.0*by[j * nx + (i-1)] + 35.0*by[j * nx + (i-2)] - 80.0*by[j * nx + (i-3)] + 30.0*by[j * nx + (i-4)] - 8.0*by[j * nx + (i-5)] + by[j * nx + (i-6)])*(1.0/60.0);
}

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] = (-147.0*bx[j * nx + i] + 360.0*bx[(j+1) * nx + i] - 450.0*bx[(j+2) * nx + i] + 400.0*bx[(j+3) * nx + i] - 225.0*bx[(j+4) * nx + i] + 72.0*bx[(j+5) * nx + i] - 10.0*bx[(j+6) * nx + i])*(1.0/60.0);
}           }

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] = ( 147.0*bx[j * nx + i] - 360.0*bx[(j-1) * nx + i] + 450.0*bx[(j-2) * nx + i] - 400.0*bx[(j-3) * nx + i] + 225.0*bx[(j-4) * nx + i] - 72.0*bx[(j-5) * nx + i] + 10.0*bx[(j-6) * nx + i])*(1.0/60.0);
}

j=1;
for (i = 0; i <= nx-2; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = (-10.0*bx[j * nx + i] - 77.0*bx[(j+1) * nx + i] + 150.0*bx[(j+2) * nx + i] - 100.0*bx[(j+3) * nx + i] + 50.0*bx[(j+4) * nx + i] - 15.0*bx[(j+5) * nx + i] + 2.0*bx[(j+6) * nx + i])*(1.0/60.0);
}

j=ny-2;
for (i = 0; i <= nx-2; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = ( 10.0*bx[j * nx + i] + 77.0*bx[(j-1) * nx + i] - 150.0*bx[(j-2) * nx + i] + 100.0*bx[(j-3) * nx + i] - 50.0*bx[(j-4) * nx + i] + 15.0*bx[(j-5) * nx + i] - 2.0*bx[(j-6) * nx + i])*(1.0/60.0);
}

j=2;
for (i = 0; i <= nx-3; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = ( 2.0*bx[j * nx + i] - 24.0*bx[(j+1) * nx + i] - 35.0*bx[(j+2) * nx + i] + 80.0*bx[(j+3) * nx + i] - 30.0*bx[(j+4) * nx + i] + 8.0*bx[(j+5) * nx + i] - bx[(j+6) * nx + i])*(1.0/60.0);
}           }

j=ny-3;
for (i = 0; i <= nx-3; i++)
{
if isnan(bx[j * nx + i]) continue;
dery[j * nx + i] = (-2.0*bx[j * nx + i] + 24.0*bx[(j-1) * nx + i] + 35.0*bx[(j-2) * nx + i] - 80.0*bx[(j-3) * nx + i] + 30.0*bx[(j-4) * nx + i] - 8.0*bx[(j-5) * nx + i] + bx[(j-6) * nx + i])*(1.0/60.0);
}

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++)
{             {
/* 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( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;

/* calculate the error in jz at all points */
jz_err[j * nx + i]=sqrt( (1.0/60.0)*bx_err[(j-3) * nx + i]*(1.0/60.0)*bx_err[(j-3) * nx + i] + (9.0/60.0)*bx_err[(j-2) * nx + i]*(9.0/60.0)*bx_err[(j-2) * nx + i] + (45.0/60.0)*bx_err[(j-1) * nx + i]*(45.0/60.0)*bx_err[(j-1) * nx + i] +
(1.0/60.0)*bx_err[(j+3) * nx + i]*(1.0/60.0)*bx_err[(j+3) * nx + i] + (9.0/60.0)*bx_err[(j+2) * nx + i]*(9.0/60.0)*bx_err[(j+2) * nx + i] + (45.0/60.0)*bx_err[(j+1) * nx + i]*(45.0/60.0)*bx_err[(j+1) * nx + i] +
(1.0/60.0)*by_err[j * nx + (i-3)]*(1.0/60.0)*by_err[j * nx + (i-3)] + (9.0/60.0)*by_err[j * nx + (i-2)]*(9.0/60.0)*by_err[j * nx + (i-2)] + (45.0/60.0)*by_err[j * nx + (i-1)]*(45.0/60.0)*by_err[j * nx + (i-1)] +
(1.0/60.0)*by_err[j * nx + (i+3)]*(1.0/60.0)*by_err[j * nx + (i+3)] + (9.0/60.0)*by_err[j * nx + (i+2)]*(9.0/60.0)*by_err[j * nx + (i+2)] + (45.0/60.0)*by_err[j * nx + (i+1)]*(45.0/60.0)*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]);
}             }
//printf("\n");
}           }

return 0;         return 0;
} }

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

/* Example function 9:  Compute quantities on Jz array */
/* Example function 9:  Compute quantities on smoothed Jz array */  // Compute mean and total current on Jz array.

// All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's
// fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels
// and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values
// of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course,
// give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels.

int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth, int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth,
float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,                     float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
 Line 826  int computeJzsmooth(float *bx, float *by
 Line 652  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++)
{           {
for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
{             {
//printf("%f ",jz_smooth[j * nx + i]);
if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
if isnan(derx[j * nx + i]) continue;                if isnan(derx[j * nx + i]) continue;
if isnan(dery[j * nx + i]) continue;                if isnan(dery[j * nx + i]) continue;
//if isnan(jz_smooth[j * nx + i]) continue;
//curl +=     (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
//us_i += fabs(jz_smooth[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
//jz_rms_err[j * nx + i] = sqrt(jz_err_squared_smooth[j * nx + i]);
//err += (jz_rms_err[j * nx + i]*jz_rms_err[j * nx + i]);
if isnan(jz[j * nx + i]) continue;                if isnan(jz[j * nx + i]) continue;
curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */                us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);                err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);
}             }
//printf("\n");
}           }

/* 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)/count_mask)*((1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.)); // 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))*((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 879  int computeJzsmooth(float *bx, float *by
 Line 701  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 904  int computeJzsmooth(float *bx, float *by
 Line 717  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++;              A += jz[j*nx+i]*bz[j*nx+i];
if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;              B += bz[j*nx+i]*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]);
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 973  int computeHelicity(float *jz_err, float
 Line 778  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 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 987  int computeHelicity(float *jz_err, float
 Line 796  int computeHelicity(float *jz_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(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 isnan(jz_err[j * nx + i]) continue;
if isnan(bz_err[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[j * nx + i] == 0.0) continue;
sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH                   sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
sum2    += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH                   sum2    += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));              err     += (jz_err[j * nx + i]*jz_err[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]) + (bz_err[j * nx + i]*bz_err[j * nx + i]*jz[j * nx + i]*jz[j * nx + i]);
}                 }
}          }
 Line 1000  int computeHelicity(float *jz_err, float
 Line 811  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(err)/count_mask)*(1/cdelt1)*(rsun_obs/rsun_ref) ; // 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(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // 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(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // 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 1020  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 1031  int computeSumAbsPerPolarity(float *jz_e
 Line 842  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 1052  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);

return 0;         return 0;
} }
 Line 1063  int computeSumAbsPerPolarity(float *jz_e
 Line 879  int computeSumAbsPerPolarity(float *jz_e
/*===========================================*/ /*===========================================*/
/* Example function 13:  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. Note that the 8*PI can come out of the integral; thus,
// the integral is over B^2 dV and the 8*PI is divided at the end.
// //
// Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
// ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm: // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
 Line 1077  int computeFreeEnergy(float *bx_err, flo
 Line 894  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,err=0.0;
if (nx <= 0 || ny <= 0) return 1;

for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
{           {
 Line 1093  int computeFreeEnergy(float *bx_err, flo
 Line 914  int computeFreeEnergy(float *bx_err, flo
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(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) continue;                  if isnan(by[j * nx + i]) continue;
sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) ) / 8.*PI;              sum  += ( ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);              sum1 += (  ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) );
//err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];              err  += 4.0*(bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])*(bx_err[j * nx + i]*bx_err[j * nx + i]) +
4.0*(by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])*(by_err[j * nx + i]*by_err[j * nx + i]);
}               }
}           }

*meanpotptr    = (sum) / (count_mask);                    /* Units are ergs per cubic centimeter */      /* Units of meanpotptr are ergs per centimeter */
*meanpotptr      = (sum1) / (count_mask*8.*PI) ;     /* Units are ergs per cubic centimeter */

/* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */         /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
*totpotptr     = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/8.*PI)) ;      *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 1120  int computeFreeEnergy(float *bx_err, flo
 Line 942  int computeFreeEnergy(float *bx_err, flo
/*===========================================*/ /*===========================================*/
/* 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 */ /* Example function 14:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */

int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,  int computeShearAngle(float *bx_err, float *by_err, float *bz_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;
float count = 0;
double dotproduct = 0.0;
double magnitude_potential = 0.0;
double magnitude_vector = 0.0;
double shear_angle = 0.0;
double denominator = 0.0;
double term1 = 0.0;
double term2 = 0.0;
double term3 = 0.0;
double sumsum = 0.0;
double err = 0.0;
double part1 = 0.0;
double part2 = 0.0;
double part3 = 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 1143  int computeShearAngle(float *bx_err, flo
 Line 982  int computeShearAngle(float *bx_err, flo
if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
if isnan(bx[j * nx + i]) continue;                  if isnan(bx[j * nx + i]) continue;
if isnan(by[j * nx + i]) continue;                  if isnan(by[j * nx + i]) continue;
/* For mean 3D shear angle, area with shear greater than 45*/              if isnan(bx_err[j * nx + i]) continue;
if isnan(by_err[j * nx + i]) continue;
if isnan(bz_err[j * nx + i]) continue;

/* For mean 3D shear angle, percentage 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]));
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("dotproduct=%f\n",dotproduct);
//printf("magnitude_potential=%f\n",magnitude_potential);
//printf("magnitude_vector=%f\n",magnitude_vector);

shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
sumsum                  += shear_angle;
//printf("shear_angle=%f\n",shear_angle);
count ++;                  count ++;
sum += shear_angle ;
err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
if (shear_angle > 45) count_mask ++;                  if (shear_angle > 45) count_mask ++;

// For the error analysis

term1 = bx[j * nx + i]*by[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bz[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i];
term2 = bx[j * nx + i]*bx[j * nx + i]*bpy[j * nx + i] - bx[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bx[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i] - bz[j * nx + i]*by[j * nx + i]*bpz[j * nx + i];
term3 = bx[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bx[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*by[j * nx + i]*bpz[j * nx + i] - by[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i];

part1 = bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i];
part2 = bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i];
part3 = bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i];

// denominator is squared
denominator = part1*part1*part1*part2*(1.0-((part3*part3)/(part1*part2)));

err = (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
(term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
(term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) ;

}               }
}           }

/* 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 = (sumsum)/(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)/count_mask)*(180./PI);
*area_w_shear_gt_45ptr   = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */
/* 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("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
//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,
float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
float *pmap, int nx1, int ny1,
int scale, float *p1pad, int nxp, int nyp, float *pmapn)

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

// set up convolution kernels
init_fresize_boxcar(&fresboxcar,1,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);

// =============== [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 (j = 0; j < ny1; j++)
{
index = j * nx1 + i;
if (rim[index] > 150)
p1p0[index]=1.0;
else
p1p0[index]=0.0;
if (rim[index] < -150)
p1n0[index]=1.0;
else
p1n0[index]=0.0;
}
}

// =============== [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, 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 (j = 0; j < ny1; j++)
{
index = j * nx1 + i;
if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
p1[index]=1.0;
else
p1[index]=0.0;
}
}

// pad p1 with zeroes so that the gaussian colvolution in step 5
// does not cut off data within hwidth of the edge

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 (j = 0; j < ny1; j++)
{
index = j * nx1 + i;
if isnan(pmapn[index]) continue;
if isnan(rim[index]) continue;
sum += pmapn[index]*abs(rim[index]);
}
}

if (sum < 1.0)
*Rparam = 0.0;
else
*Rparam = log10(sum);

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

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

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 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 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,
int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{

int nx = dims[0];
int ny = dims[1];
int i = 0;
int j = 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;
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 =========================*/

 Line 1282  void greenpot(float *bx, float *by, floa
 Line 1485  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.9 changed lines Added in v.1.38