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

Diff for /JSOC/proj/sharp/apps/sw_functions.c between version 1.35 and 1.43

version 1.35, 2015/03/02 21:41:31 version 1.43, 2021/05/26 04:45:48
Line 1 
Line 1 
   
 /*=========================================== /*===========================================
  
  The following 14 functions calculate the following spaceweather indices:   The following functions calculate these spaceweather indices from the vector magnetic field data:
  
  USFLUX Total unsigned flux in Maxwells  USFLUX Total unsigned flux in Maxwells
  MEANGAM Mean inclination angle, gamma, in degrees  MEANGAM Mean inclination angle, gamma, in degrees
Line 17 
Line 18 
  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
    CMASK The total number of pixels that contributed to the calculation of all the indices listed above
   
    And these spaceweather indices from the line-of-sight magnetic field data:
    USFLUXL Total unsigned flux in Maxwells
    MEANGBL Mean value of the line-of-sight field gradient, in Gauss/Mm
    CMASKL The total number of pixels that contributed to the calculation of USFLUXL and MEANGBL
  R_VALUE Karel Schrijver's R parameter  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
Line 830  int computeHelicity(float *jz_err, float
Line 837  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 per arcsecond.  //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square 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 1021  int computeShearAngle(float *bx_err, flo
Line 1028  int computeShearAngle(float *bx_err, flo
     }     }
     /* For mean 3D shear angle, area with shear greater than 45*/     /* For mean 3D shear angle, area with shear greater than 45*/
     *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */     *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */
   
       // For the error in the mean 3D shear angle:
       // If count_mask is 0, then we run into a divide by zero error. In this case, set *meanshear_angle_err_ptr to NAN
       // If count_mask is greater than zero, then compute the error.
       if (count_mask == 0)
           *meanshear_angle_err_ptr = NAN;
       else
     *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);     *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
  
     /* 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("ERRMSHA=%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;
 } }
  
Line 1233  int computeLorentz(float *bx, float *by
Line 1246  int computeLorentz(float *bx, float *by
  
 } }
  
   /*===========================================*/
   
   /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field */
   
   //  To compute the unsigned flux, we simply calculate
   //  flux = surface integral [(vector LOS) dot (normal vector)],
   //       = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
   //  However, since the field is radial, we will assume cos theta = 1.
   //  Therefore the pixels only need to be corrected for the projection.
   
   //  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
   //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
   //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
   //  =Gauss*cm^2
   
   int computeAbsFlux_los(float *los, int *dims, float *absFlux_los,
                          float *mean_vf_los_ptr, float *count_mask_los_ptr,
                          int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
   
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask_los = 0;
       double sum = 0.0;
       *absFlux_los = 0.0;
       *mean_vf_los_ptr = 0.0;
   
   
       if (nx <= 0 || ny <= 0) return 1;
   
           for (i = 0; i < nx; i++)
           {
              for (j = 0; j < ny; j++)
              {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if isnan(los[j * nx + i]) continue;
               sum += (fabs(los[j * nx + i]));
               count_mask_los++;
              }
           }
   
       *mean_vf_los_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       *count_mask_los_ptr  = count_mask_los;
   
       //printf("USFLUXL=%f\n",*mean_vf_los_ptr);
       //printf("CMASKL=%f\n",*count_mask_los_ptr);
   
       return 0;
   }
   
   /*===========================================*/
   /* Example function 18:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
   
   int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       *mean_derivative_los_ptr = 0.0;
   
       if (nx <= 0 || ny <= 0) return 1;
   
       /* brute force method of calculating the derivative (no consideration for edges) */
       for (i = 1; i <= nx-2; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
              derx_los[j * nx + i] = (los[j * nx + i+1] - los[j * nx + i-1])*0.5;
           }
       }
   
       /* brute force method of calculating the derivative (no consideration for edges) */
       for (i = 0; i <= nx-1; i++)
       {
           for (j = 1; j <= ny-2; j++)
           {
              dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
           }
       }
   
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
       ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
       i=0;
       for (j = 0; j <= ny-1; j++)
       {
           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;
       for (j = 0; j <= ny-1; j++)
       {
           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;
       for (i = 0; i <= nx-1; i++)
       {
           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;
       for (i = 0; i <= nx-1; i++)
       {
           dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
       }
   
   
       for (i = 0; i <= nx-1; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if isnan(los[j * nx + i])      continue;
               if isnan(los[(j+1) * nx + i])  continue;
               if isnan(los[(j-1) * nx + i])  continue;
               if isnan(los[j * nx + i-1])    continue;
               if isnan(los[j * nx + i+1])    continue;
               if isnan(derx_los[j * nx + i]) continue;
               if isnan(dery_los[j * nx + i]) continue;
               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++;
           }
       }
   
       *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
   
       //printf("MEANGBL=%f\n",*mean_derivative_los_ptr);
   
           return 0;
   }
   
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/
  
 // #include <omp.h> // #include <omp.h>

Legend:
Removed from v.1.35  
changed lines
  Added in v.1.43
Karen Tian
 Powered by
ViewCVS 0.9.4