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

version 1.29, 2014/05/16 21:55:54

version 1.31, 2014/06/05 21:27:04

Line 1048 int computeShearAngle(float *bx_err, flo

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)

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;

int index, index1;

double sum = 0.0;

double err = 0.0;

*Rparam = 0.0;

struct fresize_struct fresboxcar, fresgauss;

int scale = round(2.0/cdelt1);

struct fint_struct fints;

float sigma = 10.0/2.3548;

// set up convolution kernels

init_fresize_boxcar(&fresboxcar,1,1);

// set up convolution kernel

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

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, 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 && p1n[index] > 0)

if ((p1p[index] > 0.0) && (p1n[index] > 0.0))

p1[index]=1.0;

else

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

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

sum += pmap[index]*abs(rim[index]);

sum += pmapn[index]*abs(rim[index]);

}

Line 1120 int computeR(float *bz_err, float *los,

Line 1173 int computeR(float *bz_err, float *los,

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;

/* Multiplier */

float vectorMulti[] = {1.,-1.,1.};

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] * vectorMulti[0]) * (bz[i] * vectorMulti[2]) * k_h;

fy[i] = (by[i] * vectorMulti[1]) * (bz[i] * vectorMulti[2]) * 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;

return 0;

}

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

Legend:

Removed from v.1.29
changed lines
	Added in v.1.31

Karen Tian