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

version 1.29, 2014/05/16 21:55:54 version 1.32, 2014/09/05 21:59:48
 Line 832  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 869  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 1048  int computeShearAngle(float *bx_err, flo
 Line 1048  int computeShearAngle(float *bx_err, flo

int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1, 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 *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
float *pmap, int nx1, int ny1)               float *pmap, int nx1, int ny1,
int scale, float *p1pad, int nxp, int nyp, float *pmapn)

{ {
int nx = dims[0];     int nx = dims[0];
int ny = dims[1];     int ny = dims[1];
int i = 0;     int i = 0;
int j = 0;     int j = 0;
int index;      int index, index1;
double sum = 0.0;     double sum = 0.0;
double err = 0.0;     double err = 0.0;
*Rparam = 0.0;     *Rparam = 0.0;
struct fresize_struct fresboxcar, fresgauss;     struct fresize_struct fresboxcar, fresgauss;
int scale = round(2.0/cdelt1);      struct fint_struct fints;
float sigma = 10.0/2.3548;     float sigma = 10.0/2.3548;

// set up convolution kernels
init_fresize_boxcar(&fresboxcar,1,1);     init_fresize_boxcar(&fresboxcar,1,1);

// set up convolution kernel
init_fresize_gaussian(&fresgauss,sigma,20,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);     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 (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
 Line 1085  int computeR(float *bz_err, float *los,
 Line 1091  int computeR(float *bz_err, float *los,
}         }
}     }

// =============== [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, 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);     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 (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
{         {
index = j * nx1 + i;             index = j * nx1 + i;
if (p1p[index] > 0 && p1n[index] > 0)              if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
p1[index]=1.0;                 p1[index]=1.0;
else             else
p1[index]=0.0;                 p1[index]=0.0;
}         }
}     }

fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);      // pad p1 with zeroes so that the gaussian colvolution in step 5
// does not cut off data within hwidth of the edge

for (i = 0; i < nxp; i++)
{
for (j = 0; j < nyp; j++)
{
index = j * nxp + i;
}
}

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

// =============== [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 (i = 0; i < nx1; i++)
{     {
for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
{         {
index = j * nx1 + i;             index = j * nx1 + i;
sum += pmap[index]*abs(rim[index]);              if isnan(pmapn[index]) continue;
if isnan(rim[index]) continue;
sum += pmapn[index]*abs(rim[index]);
}         }
}     }

 Line 1116  int computeR(float *bz_err, float *los,
 Line 1171  int computeR(float *bz_err, float *los,
else     else
*Rparam = log10(sum);         *Rparam = log10(sum);

//printf("R_VALUE=%f\n",*Rparam);

free_fresize(&fresboxcar);     free_fresize(&fresboxcar);
free_fresize(&fresgauss);     free_fresize(&fresgauss);

return 0;     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 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 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 =========================*/

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