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.15 and 1.19

version 1.15, 2013/07/08 23:02:22

version 1.19, 2013/10/18 23:36:02

Line 63

// =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^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,

int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

Line 102 int computeAbsFlux(float *bz_err, float

Line 98 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("cdelt1=%f\n",cdelt1);

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

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

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

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

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

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

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

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

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

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

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

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

return 0;

}

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

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

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

if isnan(bh_err[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);

count_mask++;

}

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

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

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

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

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

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

Line 576 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

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 643 int computeJzsmooth(float *bx, float *by

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

*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 658 int computeJzsmooth(float *bx, float *by

Line 654 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 687 int computeAlpha(float *jz_err, float *b

Line 674 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 711 int computeAlpha(float *jz_err, float *b

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

//if (jz_err[j * nx + i] > abs(jz[j * nx + i]) ) continue;

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

//if (bz_err[j * nx + i] > abs(bz[j * nx + i]) ) continue;

if (bz[j * nx + i] > 0) sum3 += ( jz[j * nx + i] ); c++;

A += jz[j*nx+i]*bz[j*nx+i];

if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;

B += bz[j*nx+i]*bz[j*nx+i];

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 791 int computeHelicity(float *jz_err, float

Line 768 int computeHelicity(float *jz_err, float

*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.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 836 int computeSumAbsPerPolarity(float *jz_e

Line 813 int computeSumAbsPerPolarity(float *jz_e

{

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

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

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

if (bz[j * nx + i] > 0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);

if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);

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

Line 845 int computeSumAbsPerPolarity(float *jz_e

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

Line 881 int computeFreeEnergy(float *bx_err, flo

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

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 963 int computeShearAngle(float *bx_err, flo

Line 941 int computeShearAngle(float *bx_err, flo

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

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

}

Legend:

Removed from v.1.15
changed lines
	Added in v.1.19

Karen Tian