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.21 and 1.32

version 1.21, 2013/11/11 23:18:40 version 1.32, 2014/09/05 21:59:48
 Line 17
 Line 17
MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter     MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
TOTPOT Total photospheric 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     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    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    pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
 Line 831  int computeHelicity(float *jz_err, float
 Line 832  int computeHelicity(float *jz_err, float
/* Example function 12:  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: //  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: //  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), //  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) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
 Line 868  int computeSumAbsPerPolarity(float *jz_e
 Line 869  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));         *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);
 Line 1027  int computeShearAngle(float *bx_err, flo
 Line 1028  int computeShearAngle(float *bx_err, flo
/* The area here is a fractional area -- the % of the total area. This has no error associated with it. */         /* 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);         *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);      //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);

return 0;         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 =========================*/ /*==================KEIJI'S CODE =========================*/

 Line 1153  void greenpot(float *bx, float *by, floa
 Line 1353  void greenpot(float *bx, float *by, floa

char *sw_functions_version() // Returns CVS version of sw_functions.c char *sw_functions_version() // Returns CVS version of sw_functions.c
{ {
return strdup("\$Id\$");      return strdup("\$Id");
} }

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

Legend:
 Removed from v.1.21 changed lines Added in v.1.32

 Karen Tian Powered byViewCVS 0.9.4