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.36 and 1.40

version 1.36, 2015/10/30 18:53:02

version 1.40, 2021/05/24 22:17:06

Line 1

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

The following 14 functions calculate the following spaceweather indices:

The following functions calculate these spaceweather indices from the vector magnetic field data:

USFLUX Total unsigned flux in Maxwells

MEANGAM Mean inclination angle, gamma, in degrees

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

CMASK The total number of pixels that contributed to the calculation of all the indices listed above

And these spaceweather indices from the line-of-sight magnetic field data:

USFLUXL Total unsigned flux in Maxwells

MEANGBL Mean value of the line-of-sight field gradient, in Gauss/Mm

CMASKL The total number of pixels that contributed to the calculation of USFLUXL and MEANGBL

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

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

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

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

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

printf("bh[200,200],bh_err[200,200]=%f,%f\n",bh[200,200],bh_err[200,200]);

return 0;

}

Line 249 int computeB_total(float *bx_err, float

Line 253 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, float *err_term1, float *err_term2)

{

int nx = dims[0];

Line 269 int computeBtotalderivative(float *bt, i

Line 273 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_term1[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)])) ;

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

Line 283 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_term2[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])) ;

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

}

Line 326 int computeBtotalderivative(float *bt, i

Line 330 int computeBtotalderivative(float *bt, i

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 += err_term2[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 += 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_term1[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++;

}

Line 344 int computeBtotalderivative(float *bt, i

Line 348 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, float *err_term1, float *err_term2)

{

int nx = dims[0];

Line 364 int computeBhderivative(float *bh, float

Line 368 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_term1[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)]));

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

Line 378 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_term2[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]));

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

}

Line 421 int computeBhderivative(float *bh, float

Line 425 int computeBhderivative(float *bh, float

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 += err_term2[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 += 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_term1[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 438 int computeBhderivative(float *bh, float

Line 442 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, float *err_term1, float *err_term2)

{

int nx = dims[0];

Line 458 int computeBzderivative(float *bz, float

Line 462 int computeBzderivative(float *bz, float

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

{

derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;

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

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

Line 472 int computeBzderivative(float *bz, float

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

{

dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;

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

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

}

Line 515 int computeBzderivative(float *bz, float

Line 519 int computeBzderivative(float *bz, float

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 += err_term2[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 += 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_term1[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 992 int computeShearAngle(float *bx_err, flo

Line 996 int computeShearAngle(float *bx_err, flo

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("i,j=%d,%d\n",i,j);

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

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

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

//printf("dotproduct/(magnitude_potential*magnitude_vector)=%f\n",dotproduct/(magnitude_potential*magnitude_vector));

//printf("acos(dotproduct/(magnitude_potential*magnitude_vector))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector)));

//printf("acos(dotproduct/(magnitude_potential*magnitude_vector)*(180./PI))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI));

shear_angle = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);

sumsum += shear_angle;

Line 1029 int computeShearAngle(float *bx_err, flo

Line 1028 int computeShearAngle(float *bx_err, flo

}

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

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

// For the error in the mean 3D shear angle:

// If count_mask is 0, then we run into a divide by zero error. In this case, set *meanshear_angle_err_ptr to NAN

// If count_mask is greater than zero, then compute the error.

if (count_mask == 0)

*meanshear_angle_err_ptr = NAN;

else

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

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

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

return 0;

}

Line 1241 int computeLorentz(float *bx, float *by

Line 1246 int computeLorentz(float *bx, float *by

}

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

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

float *mean_vf_los_ptr, float *count_mask_los_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_los = 0;

double sum = 0.0;

*absFlux_los = 0.0;

*mean_vf_los_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_los++;

}

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

*count_mask_los_ptr = count_mask_los;

printf("USFLUXL=%f\n",*mean_vf_los_ptr);

printf("CMASKL=%f\n",*count_mask_los_ptr);

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("MEANGBL=%f\n",*mean_derivative_los_ptr);

return 0;

}

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

// #include <omp.h>

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

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

i2e = inx + iwindow;

if (i2s < 0){i2s = 0;}

if (i2e > nnx){i2e = nnx;}

// for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)

for (j2=j2s;j2<2;j2++){for (i2=i2s;i2<2;i2++)

for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)

{

float val1 = bztmp[nnx*j2 + i2];

float rr, r1, r2;

// r1 = (float)(i2 - inx);

// r2 = (float)(j2 - iny);

// rr = r1*r1 + r2*r2;

// if (rr < rwindow)

// {

int di, dj;

di = abs(i2 - inx);

dj = abs(j2 - iny);

sum = sum + val1 * rdist[nnx * dj + di] * dz;

printf("sum,val1,rdist[nnx * dj + di]=%f,%f,%f\n",sum,val1,rdist[nnx * dj + di]);

// }

} }

pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.

}

} } // end of OpenMP parallelism

printf("bztmp[0]=%f\n",bztmp[0]);

printf("bztmp[91746]=%f\n",bztmp[91746]);

printf("pfpot[0]=%f\n",pfpot[0]);

printf("pfpot[91746]=%f\n",pfpot[91746]);

// #pragma omp parallel for private(inx)

for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)

{

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

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

by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;

} } // end of OpenMP parallelism

printf("bx[0]=%f\n",bx[0]);

printf("by[0]=%f\n",by[0]);

printf("bx[91746]=%f\n",bx[91746]);

printf("by[91746]=%f\n",by[91746]);

free(rdist);

free(pfpot);

free(bztmp);

Line 1365 char *sw_functions_version() // Returns

Line 1507 char *sw_functions_version() // Returns

return strdup("$Id");

}

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

Legend:

Removed from v.1.36
changed lines
	Added in v.1.40

Karen Tian