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.10 and 1.17

version 1.10, 2013/05/21 00:01:34

version 1.17, 2013/07/24 02:35:08

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;

*mean_vf_ptr =0.0;

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

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

{

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

Line 89 int computeAbsFlux(float *bz_err, float

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 97 int computeAbsFlux(float *bz_err, float

Line 102 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("CMASK=%g\n",*count_mask_ptr);

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

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

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

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

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

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

//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 113 int computeBh(float *bx_err, float *by_e

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

int count_mask = 0;

double sum = 0.0;

*mean_hf_ptr = 0.0;

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

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

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

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)

{

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;

int count_mask = 0;

double sum = 0.0;

double err = 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++)

{

Line 174 int computeGamma(float *bz_err, float *b

Line 187 int computeGamma(float *bz_err, float *b

}

*mean_gamma_ptr = sum/count_mask;

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

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

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;

}

Line 187 int computeGamma(float *bz_err, float *b

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

int count_mask = 0;

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

Line 212 int computeB_total(float *bx_err, float

Line 228 int computeB_total(float *bx_err, float

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

int count_mask = 0;

double sum = 0.0;

double err = 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) */

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

Line 281 int computeBtotalderivative(float *bt, i

Line 299 int computeBtotalderivative(float *bt, 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_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (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;

}

Line 293 int computeBtotalderivative(float *bt, i

Line 311 int computeBtotalderivative(float *bt, i

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

int count_mask = 0;

double sum= 0.0;

double err =0.0;

*mean_derivative_bh_ptr = 0.0;

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

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

{

Line 361 int computeBhderivative(float *bh, float

Line 382 int computeBhderivative(float *bh, float

*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_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)

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;

}

Line 373 int computeBhderivative(float *bh, float

Line 394 int computeBhderivative(float *bh, float

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

int count_mask = 0;

double sum = 0.0;

double err = 0.0;

*mean_derivative_bz_ptr = 0.0;

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

*mean_derivative_bz_ptr = 0.0;

float sum,err = 0.0;

/* brute force method of calculating the derivative (no consideration for edges) */

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

{

Line 450 int computeBzderivative(float *bz, float

Line 474 int computeBzderivative(float *bz, float

*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_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)

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;

}

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

Line 522 int computeJz(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];

int i = 0;

int j = 0;

int count_mask = 0;

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

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

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

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

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 = 0; j <= ny-1; j++)

for (j = 1; j <= ny-2; 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

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

}

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

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

/* Example function 9: Compute quantities on 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,

float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,

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

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

int count_mask = 0;

double curl = 0.0;

double us_i = 0.0;

double err = 0.0;

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

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

{

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 isnan(derx[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;

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

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

count_mask++;

}

//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_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // 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

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;

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

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

int count_mask = 0;

double a = 0.0;

double b = 0.0;

double c = 0.0;

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

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 (j = 1; j < ny-1; j++)

{

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

//if (jz_smooth[j * nx + i] == 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) sum1 += ( bz[j * nx + i] ); a++;

if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;

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

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

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

if ((sum5 < 0) && (sum4 > 0)) sum=-sum;

*mean_alpha_ptr = sum; /* Units are 1/Mm */

*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

*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

printf("a=%d\n",a);

//printf("MEANALP=%f\n",*mean_alpha_ptr);

printf("b=%d\n",b);

//printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);

printf("d=%d\n",d);

printf("c=%d\n",c);

printf("MEANALP=%f\n",*mean_alpha_ptr);

printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);

return 0;

}

Line 738 int computeHelicity(float *jz_err, float

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

int count_mask = 0;

double sum = 0.0;

double sum2 = 0.0;

double sum_err = 0.0;

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

float sum,sum2,sum_err=0.0;

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

{

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

Line 765 int computeHelicity(float *jz_err, float

Line 797 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.) ; // error in the quantity MEANJZH

*mean_ih_err_ptr = (sqrt(sum_err*sum_err)) / (count_mask*100.0) ; // 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(sum_err*sum_err)) / (100.0) ; // 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(sum_err*sum_err)) / (100.0) ; // 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 796 int computeSumAbsPerPolarity(float *jz_e

Line 828 int computeSumAbsPerPolarity(float *jz_e

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=0;

int ny = dims[1];

int i=0;

int j=0;

int count_mask=0;

double sum1=0.0;

double sum2=0.0;

double err=0.0;

*totaljzptr=0.0;

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

*totaljzptr=0.0;

float sum1,sum2,err=0.0;

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

{

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

Line 819 int computeSumAbsPerPolarity(float *jz_e

Line 855 int computeSumAbsPerPolarity(float *jz_e

*totaljzptr = fabs(sum1) + fabs(sum2); /* Units are A */

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

}

Line 828 int computeSumAbsPerPolarity(float *jz_e

Line 864 int computeSumAbsPerPolarity(float *jz_e

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

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

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

// ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:

Line 842 int computeFreeEnergy(float *bx_err, flo

Line 879 int computeFreeEnergy(float *bx_err, flo

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;

int count_mask = 0;

double sum = 0.0;

double sum1 = 0.0;

double err = 0.0;

*totpotptr=0.0;

*meanpotptr=0.0;

float sum,err=0.0;

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

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

{

Line 858 int computeFreeEnergy(float *bx_err, flo

Line 899 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])) ) / 8.*PI;

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

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

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 += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];

count_mask++;

}

*meanpotptr = (sum) / (count_mask); /* Units are ergs per cubic centimeter */

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

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

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

//*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // 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*(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;

}

Line 888 int computeFreeEnergy(float *bx_err, flo

Line 928 int computeFreeEnergy(float *bx_err, flo

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,

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

{

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

int nx = dims[0];

int i, j;

int ny = dims[1];

int i = 0;

int j = 0;

int count_mask = 0;

double dotproduct = 0.0;

double magnitude_potential = 0.0;

double magnitude_vector = 0.0;

double shear_angle = 0.0;

double err = 0.0;

double sum = 0.0;

double count = 0.0;

*area_w_shear_gt_45ptr = 0.0;

*meanshear_angleptr = 0.0;

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 (j = 0; j < ny; j++)

Line 923 int computeShearAngle(float *bx_err, flo

Line 971 int computeShearAngle(float *bx_err, flo

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

*meanshear_angleptr = (sum)/(count); /* Units are degrees */

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

*area_w_shear_gt_45ptr = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */

*area_w_shear_gt_45ptr = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */

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

return 0;

}

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

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

/*===========END OF KEIJI'S CODE =========================*/

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

{

return strdup("$Id$");

}

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

Legend:

Removed from v.1.10
changed lines
	Added in v.1.17

Karen Tian