(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.39

version 1.35, 2015/03/02 21:41:31 version 1.39, 2021/03/18 01:53:40
Line 1 
Line 1 
   
 /*=========================================== /*===========================================
  
  The following 14 functions calculate the following spaceweather indices:  The following 14 functions calculate the following spaceweather indices:
Line 830  int computeHelicity(float *jz_err, float
Line 831  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 1022  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 1240  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,
                          float *mean_vf_ptr, float *count_mask_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 = 0;
       double sum = 0.0;
       *absFlux = 0.0;
       *mean_vf_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++;
              }
           }
   
       *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       *count_mask_ptr  = count_mask;
   
       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 ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) 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("mean_derivative_los_ptr=%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.39
Karen Tian
 Powered by
ViewCVS 0.9.4