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.4 and 1.8

version 1.4, 2012/12/18 23:05:10

version 1.8, 2013/02/09 02:39:20

Line 1

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

The following 13 functions calculate the following spaceweather indices:

The following 14 functions calculate the following spaceweather indices:

USFLUX Total unsigned flux in Maxwells

MEANGAM Mean inclination angle, gamma, in degrees

Line 81 int computeAbsFlux(float *bz, int *dims,

*absFlux = 0.0;

*mean_vf_ptr =0.0;

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

{

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

{

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

{

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

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

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

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

count_mask++;

}

Line 115 int computeBh(float *bx, float *by, floa

Line 114 int computeBh(float *bx, float *by, floa

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

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

{

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

{

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

{

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

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

bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );

Line 253 int computeBtotalderivative(float *bt, i

Line 252 int computeBtotalderivative(float *bt, i

}

/* Just some print statements

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

{

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

{

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

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

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

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

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

}

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

{

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

{

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

if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */

count_mask++;

}

Line 343 int computeBhderivative(float *bh, int *

Line 329 int computeBhderivative(float *bh, int *

}

/*Just some print statements

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

{

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

{

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

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

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

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

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

}

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

{

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

{

// if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;

if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */

count_mask++;

}

Line 437 int computeBzderivative(float *bz, int *

Line 410 int computeBzderivative(float *bz, int *

}

/*Just some print statements

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

{

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

{

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

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

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

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

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

}

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

{

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

Line 470 int computeBzderivative(float *bz, int *

Line 429 int computeBzderivative(float *bz, int *

}

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

/* Example function 8: Current Jz = (dBy/dx) - (dBx/dy) */

// In discretized space like data pixels,

Line 487 int computeBzderivative(float *bz, int *

Line 445 int computeBzderivative(float *bz, int *

// To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),

// one must perform the following unit conversions:

// (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or

// (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or

// (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),

// (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or

// (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),

// where a Tesla is represented as a Newton/Ampere*meter.

// As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).

Line 497 int computeBzderivative(float *bz, int *

Line 456 int computeBzderivative(float *bz, int *

// jz * (35.0)

// The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:

// (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(1000.)

// (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)

// =(Gauss/pix)(1/CDELT1)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)

// = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)

// =(Gauss/pix)(1/0.5)(10^-4)(4*PI*10^7)(722500)(1000.)

// =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)

int computeJz(float *bx, float *by, int *dims, float *jz,

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

int *mask, int *bitmask,

float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)

Line 514 int computeJz(float *bx, float *by, int

Line 471 int computeJz(float *bx, float *by, int

int i, j, count_mask=0;

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

*mean_jz_ptr = 0.0;

float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;

Line 574 int computeJz(float *bx, float *by, int

Line 529 int computeJz(float *bx, float *by, int

{

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

count_mask++;

}

return 0;

}

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

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

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

float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)

{

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

int i, j, count_mask=0;

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

*mean_jz_ptr = 0.0;

float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=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++)

{

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;

curl += (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */

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

us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT); /* us_i is in units of A / m^2 */

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

jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]); /* jz is in units of Gauss/pix */

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

count_mask++;

}

/* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */

mean_curl = (curl/count_mask);

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

printf("cdelt1, what is it?=%f\n",cdelt1);

Line 594 int computeJz(float *bx, float *by, int

Line 590 int computeJz(float *bx, float *by, int

}

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

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

// The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm

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

// (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz

// (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz

// avg alpha = avg Jz / avg Bz

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

// Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or

// = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)

// = 1/Mm

int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask, int *bitmask,

int computeAlpha(float *bz, int *dims, float *jz_smooth, float *mean_alpha_ptr, int *mask, int *bitmask,

float cdelt1, double rsun_ref, double rsun_obs)

{

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

int i, j, count_mask=0;

int i, j=0;

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

*mean_alpha_ptr = 0.0;

float aa, bb, cc, bznew, alpha2, sum=0.0;

float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=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[j * nx + i]) continue;

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

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

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

sum += (jz[j * nx + i] / bz[j * nx + i])*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */

if (bz[j * nx + i] > 0) sum1 += ( bz[j * nx + i]);

count_mask++;

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

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

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

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

}

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

sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */

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

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

/* Determine the sign of alpha */

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

if ((sum5 > 0) && (sum3 > 0)) sum=sum;

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

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

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

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

return 0;

}

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

/* Example function 10: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */

/* Example function 11: Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */

// The current helicity is defined as Bz*Jz and the units are G^2 / m

// The units of Jz are in G/pix; the units of Bz are in G.

Line 646 int computeAlpha(float *bz, int *dims, f

Line 667 int computeAlpha(float *bz, int *dims, f

// = G^2 / m.

int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,

int computeHelicity(float *bz, int *dims, float *jz_smooth, float *mean_ih_ptr, float *total_us_ih_ptr,

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

{

Line 659 int computeHelicity(float *bz, int *dims

Line 680 int computeHelicity(float *bz, int *dims

*mean_ih_ptr = 0.0;

float sum=0.0, sum2=0.0;

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

{

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

{

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

{

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

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

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

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

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

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

if (jz_smooth[j * nx + i] == 0.0) continue;

sum += (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);

sum += (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);

sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);

sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);

count_mask++;

}

Line 684 int computeHelicity(float *bz, int *dims

Line 705 int computeHelicity(float *bz, int *dims

}

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

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

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

// The Sum of the Absolute Value per polarity is defined as the following:

// fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.

Line 692 int computeHelicity(float *bz, int *dims

Line 713 int computeHelicity(float *bz, int *dims

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

int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,

int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr,

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

{

Line 710 int computeSumAbsPerPolarity(float *bz,

Line 731 int computeSumAbsPerPolarity(float *bz,

{

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

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

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

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

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

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

}

Line 720 int computeSumAbsPerPolarity(float *bz,

Line 741 int computeSumAbsPerPolarity(float *bz,

}

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

/* Example function 12: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */

/* Example function 13: Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */

// The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV

// automatically yields erg per cubic centimeter for an input B in Gauss.

Line 763 int computeFreeEnergy(float *bx, float *

Line 784 int computeFreeEnergy(float *bx, float *

}

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

/* Example function 13: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */

/* Example function 14: Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */

int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,

float *meanshear_angleptr, float *area_w_shear_gt_45ptr,

Legend:

Removed from v.1.4
changed lines
	Added in v.1.8

Karen Tian