JSOC/proj/sharp/apps/sw_functions.c - diff

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

version 1.18, 2013/10/01 01:57:44

version 1.38, 2020/06/30 22:38:17

Line 1

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

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

TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter

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

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

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:

50 : not all solutions agree (weak field method applied)

Line 39

Line 43

Written by Monica Bobra 15 August 2012

Potential Field code (appended) written by Keiji Hayashi

Error analysis modification 21 October 2013

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

#include <math.h>

Line 60

Line 65

// 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/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,

float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,

Line 93 int computeAbsFlux(float *bz_err, float

if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;

if isnan(bz[j * nx + i]) continue;

sum += (fabs(bz[j * nx + i]));

//printf("i,j,bz[j * nx + i]=%d,%d,%f\n",i,j,bz[j * nx + i]);

err += bz_err[j * nx + i]*bz_err[j * nx + i];

count_mask++;

}

Line 102 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_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux

*count_mask_ptr = count_mask;

//printf("cdelt1=%f\n",cdelt1);

//printf("rsun_obs=%f\n",rsun_obs);

//printf("rsun_ref=%f\n",rsun_ref);

//printf("CMASK=%g\n",*count_mask_ptr);

//printf("USFLUX=%g\n",*mean_vf_ptr);

//printf("sum=%f\n",sum);

//printf("USFLUX_err=%g\n",*mean_vf_err_ptr);

return 0;

}

Line 135 int computeBh(float *bx_err, float *by_e

Line 127 int computeBh(float *bx_err, float *by_e

{

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] );

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

Line 152 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 */

// 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,

float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)

Line 178 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_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;

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])) );

count_mask++;

}

*mean_gamma_ptr = sum/count_mask;

*mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.0); // error in the quantity (sum)/(count_mask)

*mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);

//printf("MEANGAM=%f\n",*mean_gamma_ptr);

//printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);

return 0;

Line 212 int computeB_total(float *bx_err, float

Line 219 int computeB_total(float *bx_err, float

{

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_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 225 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 ) */

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

Line 245 int computeBtotalderivative(float *bt, i

Line 267 int computeBtotalderivative(float *bt, i

for (j = 0; j <= ny-1; j++)

{

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 254 int computeBtotalderivative(float *bt, i

Line 277 int computeBtotalderivative(float *bt, i

for (j = 1; j <= ny-2; j++)

{

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;

for (j = 0; j <= ny-1; j++)

{

Line 283 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;

}

// Calculate the sum only

for (i = 0; i <= nx-1; i++)

for (i = 1; i <= nx-2; i++)

{

for (j = 0; j <= ny-1; j++)

for (j = 1; j <= ny-2; j++)

{

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(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 */

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] )) ;

count_mask++;

}

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

*mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)

*mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask);

//printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);

//printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);

return 0;

}

Line 308 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 ) */

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

Line 328 int computeBhderivative(float *bh, float

Line 362 int computeBhderivative(float *bh, float

for (j = 0; j <= ny-1; j++)

{

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 337 int computeBhderivative(float *bh, float

Line 372 int computeBhderivative(float *bh, float

for (j = 1; j <= ny-2; j++)

{

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;

for (j = 0; j <= ny-1; j++)

{

Line 372 int computeBhderivative(float *bh, float

Line 409 int computeBhderivative(float *bh, float

for (j = 0; j <= ny-1; j++)

{

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(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 */

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] )) ;

count_mask++;

}

Line 391 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 ) */

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

Line 410 int computeBzderivative(float *bz, float

Line 455 int computeBzderivative(float *bz, float

{

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;

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 420 int computeBzderivative(float *bz, float

Line 465 int computeBzderivative(float *bz, float

{

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;

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;

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;

}

i=nx-1;

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;

}

j=0;

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;

}

j=ny-1;

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;

}

Line 460 int computeBzderivative(float *bz, float

Line 502 int computeBzderivative(float *bz, float

{

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 ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) 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(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 */

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] )) ;

count_mask++;

}

Line 518 int computeBzderivative(float *bz, float

Line 565 int computeBzderivative(float *bz, float

// 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 *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 537 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++)

{

if isnan(by[j * nx + i]) continue;

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 546 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++)

{

if isnan(bx[j * nx + i]) continue;

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;

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;

}

i=nx-1;

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;

}

j=0;

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;

}

j=ny-1;

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;

}

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

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] ) ;

// the next 7 lines can be used with a for loop that goes from i=0;i<=nx-1 and j=0;j<=ny-1.

//int i1, j1,i2, j2;

//i1 = i + 1 ; if (i1 >nx-1){i1=nx-1;}

//j1 = j + 1 ; if (j1 >ny-1){j1=ny-1;}

//i2 = i - 1; if (i2 < 0){i2 = 0;}

//j2 = j - 1; if (j2 < 0){i2 = 0;}

//jz_err[j * nx + i] = 0.5*sqrt( (bx_err[j1 * nx + i]*bx_err[j1 * nx + i]) + (bx_err[j2 * nx + i]*bx_err[j2 * nx + i]) +

// (by_err[j * nx + i1]*by_err[j * nx + i1]) + (by_err[j * nx + i2]*by_err[j * nx + i2]) ) ;

jz_err[j * nx + i] = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +

(by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;

jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);

count_mask++;

}

return 0;

Line 609 int computeJz(float *bx_err, float *by_e

Line 643 int computeJz(float *bx_err, float *by_e

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

/* Example function 9: Compute quantities on Jz array */

// Compute mean and total current on Jz array.

Line 648 int computeJzsmooth(float *bx, float *by

Line 681 int computeJzsmooth(float *bx, float *by

/* 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_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_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_err=%f\n",*mean_jz_err_ptr);

Line 668 int computeJzsmooth(float *bx, float *by

Line 701 int computeJzsmooth(float *bx, float *by

/* Example function 10: Twist Parameter, alpha */

// The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation

// for alpha is calculated in the following way (different from Leka and Barnes' approach):

// 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 Jz are in Gauss/pix; the units of Bz are in Gauss.

Line 697 int computeAlpha(float *jz_err, float *b

Line 721 int computeAlpha(float *jz_err, float *b

int ny = dims[1];

int i = 0;

int j = 0;

int count_mask = 0;

double alpha_total = 0.0;

double a = 0.0;

double C = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));

double b = 0.0;

double total = 0.0;

double c = 0.0;

double A = 0.0;

double d = 0.0;

double B = 0.0;

double sum1 = 0.0;

double sum2 = 0.0;

double sum3 = 0.0;

double sum4 = 0.0;

double sum = 0.0;

double sum5 = 0.0;

double sum6 = 0.0;

double sum_err = 0.0;

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

for (i = 1; i < nx-1; i++)

{

for (j = 1; j < ny-1; j++)

Line 721 int computeAlpha(float *jz_err, float *b

Line 736 int computeAlpha(float *jz_err, float *b

if isnan(jz[j * nx + i]) continue;

if isnan(bz[j * nx + i]) 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) 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[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.));

count_mask++;

}

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.0); // 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("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;

}

Line 773 int computeHelicity(float *jz_err, float

Line 785 int computeHelicity(float *jz_err, float

int count_mask = 0;

double sum = 0.0;

double sum2 = 0.0;

double sum_err = 0.0;

double err = 0.0;

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

Line 784 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 isnan(jz[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 (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

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]);

count_mask++;

}

Line 797 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_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.0) ; // 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.0) ; // 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.0) ; // 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_err=%f\n",*mean_ih_err_ptr);

Line 817 int computeHelicity(float *jz_err, float

Line 831 int computeHelicity(float *jz_err, float

/* Example function 12: Sum of Absolute Value per polarity */

// 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:

// 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)

Line 854 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));

//printf("SAVNCPP=%g\n",*totaljzptr);

//printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);

Line 900 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 isnan(bx[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])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);

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

sum1 += ( ((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])) );

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 += (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]);

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]);

count_mask++;

}

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

/* Units of meanpotptr are ergs per centimeter */

*meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)

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

*meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)

/* 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)/(8.*PI);

Line 926 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 */

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,

float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)

{

int nx = dims[0];

int ny = dims[1];

int i = 0;

int j = 0;

int count_mask = 0;

float count_mask = 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 sum = 0.0;

double part1 = 0.0;

double count = 0.0;

double part2 = 0.0;

double part3 = 0.0;

*area_w_shear_gt_45ptr = 0.0;

*meanshear_angleptr = 0.0;

Line 957 int computeShearAngle(float *bx_err, flo

Line 982 int computeShearAngle(float *bx_err, flo

if isnan(bz[j * nx + i]) continue;

if isnan(bx[j * nx + i]) continue;

if isnan(by[j * nx + i]) continue;

/* For mean 3D shear angle, area with shear greater than 45*/

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]);

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]) );

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

sumsum += shear_angle;

//printf("shear_angle=%f\n",shear_angle);

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

// 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*/

*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.0);/* 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. */

*area_w_shear_gt_45ptr = (count_mask/(count))*(100.0);

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

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

//printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);

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

// step i: zero p1pad

for (i = 0; i < nxp; i++)

{

for (j = 0; j < nyp; j++)

{

index = j * nxp + i;

p1pad[index]=0.0;

}

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

p1pad[index1]=p1[index];

}

// =============== [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 *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]));

count_mask++;

}

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

*count_mask_ptr = count_mask;

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 */

count_mask++;

}

*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 =========================*/

Line 1099 void greenpot(float *bx, float *by, floa

Line 1488 void greenpot(float *bx, float *by, floa

char *sw_functions_version() // Returns CVS version of sw_functions.c

{

return strdup("$Id$");

return strdup("$Id");

}

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

Legend:

Removed from v.1.18
changed lines
	Added in v.1.38

Karen Tian