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.20 and 1.35

version 1.20, 2013/11/02 19:53:05

version 1.35, 2015/03/02 21:41:31

Line 17

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

Line 245 int computeB_total(float *bx_err, float

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

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

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

Line 307 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 = 1; i <= nx-2; i++)

{

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

{

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

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

if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;

if isnan(dery_bt[j * nx + i]) 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 += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ))+

err += err_termBt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] ))+

(((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] )) ;

err_termAt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i] + dery_bt[j * nx + i]*dery_bt[j * nx + i] )) ;

count_mask++;

}

Line 332 int computeBtotalderivative(float *bt, i

Line 341 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 352 int computeBhderivative(float *bh, float

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

Line 371 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 396 int computeBhderivative(float *bh, float

Line 408 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 += (((bh[(j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ))+

err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] ))+

(((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] )) ;

err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i] + dery_bh[j * nx + i]*dery_bh[j * nx + i] )) ;

count_mask++;

}

Line 416 int computeBhderivative(float *bh, float

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

Line 454 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 445 int computeBzderivative(float *bz, float

Line 464 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 485 int computeBzderivative(float *bz, float

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

(((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)])) / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i] + dery_bz[j * nx + i]*dery_bz[j * nx + i] )) ;

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

Line 564 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 563 int computeJz(float *bx_err, float *by_e

Line 583 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 572 int computeJz(float *bx_err, float *by_e

Line 592 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( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +

jz_err[j * nx + i] = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;

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

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

count_mask++;

}

return 0;

Line 797 int computeHelicity(float *jz_err, float

Line 814 int computeHelicity(float *jz_err, float

*total_us_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ; // error in the quantity TOTUSJH

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

}

Line 813 int computeHelicity(float *jz_err, float

Line 830 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 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 850 int computeSumAbsPerPolarity(float *jz_e

Line 867 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 926 int computeFreeEnergy(float *bx_err, flo

Line 943 int computeFreeEnergy(float *bx_err, flo

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 err = 0.0;

double sum = 0.0;

double count = 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;

Line 959 int computeShearAngle(float *bx_err, flo

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

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

sum += shear_angle ;

sumsum += shear_angle;

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

count ++;

if (shear_angle > 45) count_mask ++;

/* For the error analysis*/

// For the error analysis

// terms 1,2, and 3 are not squared

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

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] + bz[j * nx + i]*(bz[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*bpz[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] - bz[j * nx + i]*bpy[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 = (bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) *

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

(bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) *

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

(bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] * bpz[j * nx + i]*bpz[j * nx + i]) *

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

( 1-( ((bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i]) *

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

(bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i])) /

((bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]) * (bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i])) ));

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

((term2*term2*by_err[j * nx + i]*by_err[j * nx + i])/(denominator)) +

((term3*term3*bz_err[j * nx + i]*bz_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))/(count))*(180./PI);

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

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

}

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

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

}

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

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

Line 1351 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.20
changed lines
	Added in v.1.35

Karen Tian