 Return to smarp_functions.c CVS log  Up to [Development] / JSOC / proj / sharp / apps

### Diff for /JSOC/proj/sharp/apps/smarp_functions.c between version 1.2 and 1.3

version 1.2, 2020/05/04 21:11:07 version 1.3, 2020/06/29 22:19:51
 Line 1
 Line 1

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

The following 3 functions calculate the following spaceweather indices:   The following 3 functions calculate the following spaceweather indices on line-of-sight
magnetic field data:

USFLUX Total unsigned flux in Maxwells  USFLUX Total unsigned flux in Maxwells
MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm   MEANGBZ Mean value of the line-of-sight (approximately vertical in remapped and reprojected data) field gradient, in Gauss/Mm
R_VALUE Karel Schrijver's R parameter   R_VALUE R parameter (Schrijver, 2007)

The indices are calculated on the pixels in which the bitmap segment is greater than 36.  The indices are calculated on the pixels in which the bitmap segment is greater than 36.

The SMARP bitmap has 13 unique values because they describe three different characteristics:  The SMARP bitmap has 13 unique values because they describe three different characteristics:
the location of the pixel (whether the pixel is off disk or a member of the patch), the  the location of the pixel (whether the pixel is off disk or a member of the patch), the
character of the magnetic field (quiet or active), and the character of the continuum  character of the magnetic field (quiet or active), and the character of the continuum
intensity data (quiet, part of faculae, part of a sunspot).   intensity data (quiet, faculae, sunspot).

Here are all the possible values:  Here are all the possible values:

 Line 58
 Line 59
/* Example function 1: Compute total unsigned flux in units of G/cm^2 */ /* Example function 1: Compute total unsigned flux in units of G/cm^2 */

//  To compute the unsigned flux, we simply calculate //  To compute the unsigned flux, we simply calculate
//  flux = surface integral [(vector Bz) dot (normal vector)],  //  flux = surface integral [(vector LOS) dot (normal vector)],
//       = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].  //       = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
//  However, since the field is radial, we will assume cos theta = 1. //  However, since the field is radial, we will assume cos theta = 1.
//  Therefore the pixels only need to be corrected for the projection. //  Therefore the pixels only need to be corrected for the projection.

 Line 68
 Line 69
//  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2 //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
//  =Gauss*cm^2 //  =Gauss*cm^2

int computeAbsFlux(float *bz, int *dims, float *absFlux,  int computeAbsFlux_los(float *los, int *dims, float *absFlux,
int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                    int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

 Line 91  int computeAbsFlux(float *bz, int *dims,
 Line 92  int computeAbsFlux(float *bz, int *dims,
for (j = 0; j < ny; j++)            for (j = 0; j < ny; j++)
{            {
if ( bitmask[j * nx + i] < 36 ) continue;             if ( bitmask[j * nx + i] < 36 ) continue;
if isnan(bz[j * nx + i]) continue;              if isnan(los[j * nx + i]) continue;
sum += (fabs(bz[j * nx + i]));              sum += (fabs(los[j * nx + i]));
}            }
}         }
 Line 109  int computeAbsFlux(float *bz, int *dims,
 Line 110  int computeAbsFlux(float *bz, int *dims,
} }

/*===========================================*/ /*===========================================*/
/* Example function 2:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */  /* Example function 2:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */

int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *bitmask, float *derx_bz, float *dery_bz)  int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
{ {

int nx = dims;     int nx = dims;
 Line 120  int computeBzderivative(float *bz, int *
 Line 121  int computeBzderivative(float *bz, int *
int j = 0;     int j = 0;
double sum = 0.0;     double sum = 0.0;
*mean_derivative_bz_ptr = 0.0;      *mean_derivative_los_ptr = 0.0;

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

 Line 129  int computeBzderivative(float *bz, int *
 Line 130  int computeBzderivative(float *bz, int *
{     {
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{         {
derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;             derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
}         }
}     }

 Line 138  int computeBzderivative(float *bz, int *
 Line 139  int computeBzderivative(float *bz, int *
{     {
for (j = 1; j <= ny-2; j++)         for (j = 1; j <= ny-2; j++)
{         {
dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;             dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
}         }
}     }

 Line 148  int computeBzderivative(float *bz, int *
 Line 149  int computeBzderivative(float *bz, int *
i=0;     i=0;
for (j = 0; j <= ny-1; j++)     for (j = 0; j <= ny-1; j++)
{     {
derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;          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;     i=nx-1;
for (j = 0; j <= ny-1; j++)     for (j = 0; j <= ny-1; j++)
{     {
derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;          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;     j=0;
for (i = 0; i <= nx-1; i++)     for (i = 0; i <= nx-1; i++)
{     {
dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;          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;     j=ny-1;
for (i = 0; i <= nx-1; i++)     for (i = 0; i <= nx-1; i++)
{     {
dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;          dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
}     }

 Line 175  int computeBzderivative(float *bz, int *
 Line 176  int computeBzderivative(float *bz, int *
for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
{         {
if ( bitmask[j * nx + i] < 36 ) continue;             if ( bitmask[j * nx + i] < 36 ) continue;
if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;              if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue;
if isnan(bz[j * nx + i])      continue;              if isnan(los[j * nx + i])      continue;
if isnan(bz[(j+1) * nx + i])  continue;              if isnan(los[(j+1) * nx + i])  continue;
if isnan(bz[(j-1) * nx + i])  continue;              if isnan(los[(j-1) * nx + i])  continue;
if isnan(bz[j * nx + i-1])    continue;              if isnan(los[j * nx + i-1])    continue;
if isnan(bz[j * nx + i+1])    continue;              if isnan(los[j * nx + i+1])    continue;
if isnan(derx_bz[j * nx + i]) continue;              if isnan(derx_los[j * nx + i]) continue;
if isnan(dery_bz[j * nx + i]) continue;              if isnan(dery_los[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 */              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 */
}         }
}     }

*mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram      *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
//printf("mean_derivative_bz_ptr=%f\n",*mean_derivative_bz_ptr);      //printf("mean_derivative_los_ptr=%f\n",*mean_derivative_los_ptr);

return 0;         return 0;
} }
 Line 205  int computeBzderivative(float *bz, int *
 Line 206  int computeBzderivative(float *bz, int *
// are the dimensions of the input array, fsample() will usually result // are the dimensions of the input array, fsample() will usually result
// in a segfault (though not always, depending on how the segfault was accessed.) // in a segfault (though not always, depending on how the segfault was accessed.)

int computeR(float *los, int *dims, float *Rparam, float cdelt1,  int computeR_los(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 scale, float *p1pad, int nxp, int nyp, float *pmapn)

Legend:
 Removed from v.1.2 changed lines Added in v.1.3