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

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

version 1.1, 2018/03/28 01:27:06 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 conf_disambig segment is greater than 70 and   The indices are calculated on the pixels in which the bitmap segment is greater than 36.
  pixels in which the bitmap segment is greater than 42. These ranges are selected because the CCD  
  coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data   The SMARP bitmap has 13 unique values because they describe three different characteristics:
  prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI   the location of the pixel (whether the pixel is off disk or a member of the patch), the
  contain nearest-neighbor bitmaps).   character of the magnetic field (quiet or active), and the character of the continuum
    intensity data (quiet, faculae, sunspot).
  In the CCD coordinates, this means that we are selecting the pixels that that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:  
    Here are all the possible values:
   
  For bitmap:   Value   Keyword         Definition
  1  : weak field outside smooth bounding curve   0       OFFDISK         The pixel is located somewhere off the solar disk.
  2  : strong field outside smooth bounding curve   1       QUIET           The pixel is associated with weak line-of-sight magnetic field.
  33 : weak field inside smooth bounding curve   2       ACTIVE          The pixel is associated with strong line-of-sight magnetic field.
  34 : strong field inside smooth bounding curve   4       SPTQUIET        The pixel is associated with no features in the continuum intensity data.
    8       SPTFACUL        The pixel is associated with faculae in the continuum intensity data.
  Written by Monica Bobra 15 August 2012   16      SPTSPOT         The pixel is associated with sunspots in the continuum intensity data.
  Potential Field code (appended) written by Keiji Hayashi   32      ON_PATCH        The pixel is inside the patch.
  Error analysis modification 21 October 2013  
    Here are all the possible permutations:
   
    Value   Definition
    0       The pixel is located somewhere off the solar disk.
    5       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
    9       The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
    17      The pixel is located outside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
    6       The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data.
    10      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
    18      The pixel is located outside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.
    37      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and no features in the continuum intensity data.
    41      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and faculae in the continuum intensity data.
    49      The pixel is located inside the patch, associated with weak line-of-sight magnetic field, and sunspots in the continuum intensity data.
    38      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and no features in the continuum intensity data.
    42      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and faculae in the continuum intensity data.
    50      The pixel is located inside the patch, associated with strong line-of-sight magnetic field, and sunspots in the continuum intensity data.
   
    Thus pixels with a value greater than 36 are located inside the patch.
   
    Written by Monica Bobra
  
  ===========================================*/  ===========================================*/
 #include <math.h> #include <math.h>
Line 38 
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 48 
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,
                    float *mean_vf_ptr, float *count_mask_ptr,                    float *mean_vf_ptr, float *count_mask_ptr,
                    int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                    int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
Line 70  int computeAbsFlux(float *bz, int *dims,
Line 91  int computeAbsFlux(float *bz, int *dims,
         {         {
            for (j = 0; j < ny; j++)            for (j = 0; j < ny; j++)
            {            {
             if ( bitmask[j * nx + i] < 42 ) 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]));
             count_mask++;             count_mask++;
            }            }
         }         }
Line 89  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[0];     int nx = dims[0];
Line 100  int computeBzderivative(float *bz, int *
Line 121  int computeBzderivative(float *bz, int *
     int j = 0;     int j = 0;
     int count_mask = 0;     int count_mask = 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 109  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 118  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 128  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 154  int computeBzderivative(float *bz, int *
Line 175  int computeBzderivative(float *bz, int *
     {     {
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
         {         {
             if ( bitmask[j * nx + i] < 42 ) 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 */
             count_mask++;             count_mask++;
         }         }
     }     }
  
     *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 185  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.1  
changed lines
  Added in v.1.3
Karen Tian
 Powered by
ViewCVS 0.9.4