(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.19 and 1.37

version 1.19, 2013/10/18 23:36:02 version 1.37, 2015/10/30 20:21:28
Line 1 
Line 1 
   
 /*=========================================== /*===========================================
  
    The following 14 functions calculate the following spaceweather indices:    The following 14 functions calculate the following spaceweather indices:
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
    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
    pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD    pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
    coordinate bitmaps are interpolated.   coordinate bitmaps are interpolated for certain data (at the time of this CVS submit, all data
    prior to 2013.08.21_17:24:00_TAI contain interpolated bitmaps; data post-2013.08.21_17:24:00_TAI
    contain nearest-neighbor bitmaps).
  
    In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig    In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
    and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:   and the pixels that equal 33 or 34 in bitmap. Here are the definitions of the pixel values:
  
    For conf_disambig:    For conf_disambig:
    50 : not all solutions agree (weak field method applied)    50 : not all solutions agree (weak field method applied)
Line 39 
Line 43 
  
    Written by Monica Bobra 15 August 2012    Written by Monica Bobra 15 August 2012
    Potential Field code (appended) written by Keiji Hayashi    Potential Field code (appended) written by Keiji Hayashi
    Error analysis modification 21 October 2013
  
 ===========================================*/ ===========================================*/
 #include <math.h> #include <math.h>
Line 60 
Line 65 
 //  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel. //  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). //  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/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2  //  =Gauss*cm^2
 //  =(1.30501e15)Gauss*cm^2  
  
 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux, int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
                    float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,                    float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
Line 89  int computeAbsFlux(float *bz_err, float
Line 93  int computeAbsFlux(float *bz_err, float
                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(bz[j * nx + i]) continue;                   if isnan(bz[j * nx + i]) continue;
                   sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
                   //printf("i,j,bz[j * nx + i]=%d,%d,%f\n",i,j,bz[j * nx + i]);  
                   err += bz_err[j * nx + i]*bz_err[j * nx + i];                   err += bz_err[j * nx + i]*bz_err[j * nx + i];
                   count_mask++;                   count_mask++;
                 }                 }
Line 98  int computeAbsFlux(float *bz_err, float
Line 101  int computeAbsFlux(float *bz_err, float
      *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;      *mean_vf_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
      *mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux      *mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux
      *count_mask_ptr  = count_mask;      *count_mask_ptr  = count_mask;
      //printf("cdelt1=%f\n",cdelt1);  
      //printf("rsun_obs=%f\n",rsun_obs);  
      //printf("rsun_ref=%f\n",rsun_ref);  
      //printf("CMASK=%g\n",*count_mask_ptr);  
      //printf("USFLUX=%g\n",*mean_vf_ptr);  
      //printf("sum=%f\n",sum);  
      //printf("USFLUX_err=%g\n",*mean_vf_err_ptr);  
      return 0;      return 0;
 } }
  
Line 131  int computeBh(float *bx_err, float *by_e
Line 127  int computeBh(float *bx_err, float *by_e
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 if isnan(bx[j * nx + i]) continue;              if isnan(bx[j * nx + i])
                 if isnan(by[j * nx + i]) continue;              {
                   bh[j * nx + i] = NAN;
                   bh_err[j * nx + i] = NAN;
                   continue;
               }
               if isnan(by[j * nx + i])
               {
                   bh[j * nx + i] = NAN;
                   bh_err[j * nx + i] = NAN;
                   continue;
               }
                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
                 sum += bh[j * nx + i];                 sum += bh[j * nx + i];
                 bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];                 bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];
Line 148  int computeBh(float *bx_err, float *by_e
Line 154  int computeBh(float *bx_err, float *by_e
 /*===========================================*/ /*===========================================*/
 /* Example function 3: Calculate Gamma in units of degrees */ /* Example function 3: Calculate Gamma in units of degrees */
 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI) // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
 // Redo calculation in radians for error analysis (since derivatives are only true in units of radians).  //
   // Error analysis calculations are done in radians (since derivatives are only true in units of radians),
   // and multiplied by (180./PI) at the end for consistency in units.
  
 int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims, int computeGamma(float *bz_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
                  float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)                  float *mean_gamma_ptr, float *mean_gamma_err_ptr, int *mask, int *bitmask)
Line 174  int computeGamma(float *bz_err, float *b
Line 182  int computeGamma(float *bz_err, float *b
                     if isnan(bz[j * nx + i]) continue;                     if isnan(bz[j * nx + i]) continue;
                     if isnan(bz_err[j * nx + i]) continue;                     if isnan(bz_err[j * nx + i]) continue;
                     if isnan(bh_err[j * nx + i]) continue;                     if isnan(bh_err[j * nx + i]) continue;
                   if isnan(bh[j * nx + i]) continue;
                     if (bz[j * nx + i] == 0) continue;                     if (bz[j * nx + i] == 0) continue;
                     sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);                     sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
                     err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i])))  * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);                  err += (1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))))*(1/(1+((bh[j * nx + i]*bh[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])))) *
                   ( ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) +
                    ((bh[j * nx + i]*bh[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]*bz[j * nx + i])) );
                     count_mask++;                     count_mask++;
                   }                   }
               }               }
           }           }
  
      *mean_gamma_ptr = sum/count_mask;      *mean_gamma_ptr = sum/count_mask;
      *mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.0); // error in the quantity (sum)/(count_mask)      *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
      //printf("MEANGAM=%f\n",*mean_gamma_ptr);      //printf("MEANGAM=%f\n",*mean_gamma_ptr);
      //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);      //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
      return 0;      return 0;
Line 208  int computeB_total(float *bx_err, float
Line 219  int computeB_total(float *bx_err, float
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 if isnan(bx[j * nx + i]) continue;              if isnan(bx[j * nx + i])
                 if isnan(by[j * nx + i]) continue;              {
                 if isnan(bz[j * nx + i]) continue;                  bt[j * nx + i] = NAN;
                   bt_err[j * nx + i] = NAN;
                   continue;
               }
               if isnan(by[j * nx + i])
               {
                   bt[j * nx + i] = NAN;
                   bt_err[j * nx + i] = NAN;
                   continue;
               }
               if isnan(bz[j * nx + i])
               {
                   bt[j * nx + i] = NAN;
                   bt_err[j * nx + i] = NAN;
                   continue;
               }
                 bt[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]);                 bt[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]);
                 bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];                 bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];
               }               }
Line 221  int computeB_total(float *bx_err, float
Line 247  int computeB_total(float *bx_err, float
 /*===========================================*/ /*===========================================*/
 /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
  
 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_termAt, float *err_termBt)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 241  int computeBtotalderivative(float *bt, i
Line 267  int computeBtotalderivative(float *bt, i
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;                 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
              err_termAt[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
               }               }
           }           }
  
Line 250  int computeBtotalderivative(float *bt, i
Line 277  int computeBtotalderivative(float *bt, i
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;                 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
              err_termBt[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      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;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
Line 279  int computeBtotalderivative(float *bt, i
Line 308  int computeBtotalderivative(float *bt, i
              dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;              dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
           }           }
  
       // Calculate the sum only
         for (i = 0; i <= nx-1; i++)      for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 0; j <= ny-1; j++)          for (j = 1; j <= ny-2; j++)
             {             {
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
               if ( (derx_bt[j * nx + i] + dery_bt[j * nx + i]) == 0) continue;
               if isnan(bt[j * nx + i])      continue;
               if isnan(bt[(j+1) * nx + i])  continue;
               if isnan(bt[(j-1) * nx + i])  continue;
               if isnan(bt[j * nx + i-1])    continue;
               if isnan(bt[j * nx + i+1])    continue;
               if isnan(bt_err[j * nx + i])  continue;
                if isnan(derx_bt[j * nx + i]) continue;                if isnan(derx_bt[j * nx + i]) continue;
                if isnan(dery_bt[j * nx + i]) continue;                if isnan(dery_bt[j * nx + i]) continue;
                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */
                err += (2.0)*bt_err[j * nx + i]*bt_err[j * nx + i];              err += err_termBt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+
                      err_termAt[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
  
         *mean_derivative_btotal_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram      *mean_derivative_btotal_ptr     = (sum)/(count_mask);
         *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)      *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask);
         //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);         //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
         //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);         //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
   
         return 0;         return 0;
 } }
  
Line 304  int computeBtotalderivative(float *bt, i
Line 342  int computeBtotalderivative(float *bt, i
 /*===========================================*/ /*===========================================*/
 /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
  
 int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)  int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_termAh, float *err_termBh)
 { {
  
      int nx = dims[0];      int nx = dims[0];
Line 324  int computeBhderivative(float *bh, float
Line 362  int computeBhderivative(float *bh, float
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;                 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
              err_termAh[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
               }               }
           }           }
  
Line 333  int computeBhderivative(float *bh, float
Line 372  int computeBhderivative(float *bh, float
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;                 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
             err_termBh[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      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;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
Line 368  int computeBhderivative(float *bh, float
Line 409  int computeBhderivative(float *bh, float
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
               if ( (derx_bh[j * nx + i] + dery_bh[j * nx + i]) == 0) continue;
               if isnan(bh[j * nx + i])      continue;
               if isnan(bh[(j+1) * nx + i])  continue;
               if isnan(bh[(j-1) * nx + i])  continue;
               if isnan(bh[j * nx + i-1])    continue;
               if isnan(bh[j * nx + i+1])    continue;
               if isnan(bh_err[j * nx + i])  continue;
                if isnan(derx_bh[j * nx + i]) continue;                if isnan(derx_bh[j * nx + i]) continue;
                if isnan(dery_bh[j * nx + i]) continue;                if isnan(dery_bh[j * nx + i]) continue;
                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */
                err += (2.0)*bh_err[j * nx + i]*bh_err[j * nx + i];              err += err_termBh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+
                      err_termAh[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
Line 387  int computeBhderivative(float *bh, float
Line 436  int computeBhderivative(float *bh, float
 /*===========================================*/ /*===========================================*/
 /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
  
 int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)  int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_termA, float *err_termB)
 { {
  
         int nx = dims[0];         int nx = dims[0];
Line 406  int computeBzderivative(float *bz, float
Line 455  int computeBzderivative(float *bz, float
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 if isnan(bz[j * nx + i]) continue;  
                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
              err_termA[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
               }               }
           }           }
  
Line 416  int computeBzderivative(float *bz, float
Line 465  int computeBzderivative(float *bz, float
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 if isnan(bz[j * nx + i]) continue;  
                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
              err_termB[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      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;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[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++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(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++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
           }           }
  
Line 456  int computeBzderivative(float *bz, float
Line 502  int computeBzderivative(float *bz, float
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;  
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
               if ( (derx_bz[j * nx + i] + dery_bz[j * nx + i]) == 0) continue;
                if isnan(bz[j * nx + i]) continue;                if isnan(bz[j * nx + i]) continue;
                //if isnan(bz_err[j * nx + i]) continue;              if isnan(bz[(j+1) * nx + i])  continue;
               if isnan(bz[(j-1) * nx + i])  continue;
               if isnan(bz[j * nx + i-1])    continue;
               if isnan(bz[j * nx + i+1])    continue;
               if isnan(bz_err[j * nx + i])  continue;
                if isnan(derx_bz[j * nx + i]) continue;                if isnan(derx_bz[j * nx + i]) continue;
                if isnan(dery_bz[j * nx + i]) continue;                if isnan(dery_bz[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_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */
                err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];              err += err_termB[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +
                      err_termA[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
Line 514  int computeBzderivative(float *bz, float
Line 565  int computeBzderivative(float *bz, float
 //              float *noiseby, float *noisebz) //              float *noiseby, float *noisebz)
  
 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared, int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
               int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)
  
  
 { {
Line 533  int computeJz(float *bx_err, float *by_e
Line 584  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  if isnan(by[j * nx + i]) continue;  
                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
              err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);
               }               }
           }           }
  
Line 542  int computeJz(float *bx_err, float *by_e
Line 593  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                  if isnan(bx[j * nx + i]) continue;  
                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
              err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
               }               }
           }           }
  
         // consider the edges      /* 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;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[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++)
           {           {
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(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++)
           {           {
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
           }           }
  
         for (i = 1; i <= nx-2; i++)  
       for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 1; j <= ny-2; j++)          for (j = 0; j <= ny-1; j++)
             {             {
                // calculate jz at all points                // calculate jz at all points
   
                jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix                jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix
                jz_err[j * nx + i]        = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +              jz_err[j * nx + i]        = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;
                                                      (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;  
                jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);                jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
                count_mask++;                count_mask++;
   
             }             }
           }           }
         return 0;         return 0;
Line 595  int computeJz(float *bx_err, float *by_e
Line 643  int computeJz(float *bx_err, float *by_e
  
 /*===========================================*/ /*===========================================*/
  
   
 /* Example function 9:  Compute quantities on Jz array */ /* Example function 9:  Compute quantities on Jz array */
 // Compute mean and total current on Jz array. // Compute mean and total current on Jz array.
  
Line 634  int computeJzsmooth(float *bx, float *by
Line 681  int computeJzsmooth(float *bx, float *by
  
         /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */         /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */
         *mean_jz_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // error in the quantity MEANJZD      *mean_jz_err_ptr = (sqrt(err)/count_mask)*((1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.)); // error in the quantity MEANJZD
  
         *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */         *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
         *us_i_err_ptr    = (sqrt(err))*fabs((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ      *us_i_err_ptr    = (sqrt(err))*((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
  
         //printf("MEANJZD=%f\n",*mean_jz_ptr);         //printf("MEANJZD=%f\n",*mean_jz_ptr);
         //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);         //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
Line 690  int computeAlpha(float *jz_err, float *b
Line 737  int computeAlpha(float *jz_err, float *b
                 if isnan(bz[j * nx + i])   continue;                 if isnan(bz[j * nx + i])   continue;
                 if (jz[j * nx + i] == 0.0) continue;                 if (jz[j * nx + i] == 0.0) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
                 //if (jz_err[j * nx + i] > abs(jz[j * nx + i]) ) continue;  
                 //if (bz_err[j * nx + i] > abs(bz[j * nx + i]) ) continue;  
                 A += jz[j*nx+i]*bz[j*nx+i];                 A += jz[j*nx+i]*bz[j*nx+i];
                 B += bz[j*nx+i]*bz[j*nx+i];                 B += bz[j*nx+i]*bz[j*nx+i];
               }               }
Line 740  int computeHelicity(float *jz_err, float
Line 785  int computeHelicity(float *jz_err, float
         int count_mask = 0;         int count_mask = 0;
         double sum = 0.0;         double sum = 0.0;
         double sum2 = 0.0;         double sum2 = 0.0;
         double sum_err = 0.0;          double err     = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
Line 751  int computeHelicity(float *jz_err, float
Line 796  int computeHelicity(float *jz_err, float
                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(jz[j * nx + i]) continue;                   if isnan(jz[j * nx + i]) continue;
                   if isnan(bz[j * nx + i]) continue;                   if isnan(bz[j * nx + i]) continue;
               if isnan(jz_err[j * nx + i]) continue;
               if isnan(bz_err[j * nx + i]) continue;
                   if (bz[j * nx + i] == 0.0) continue;                   if (bz[j * nx + i] == 0.0) continue;
                   if (jz[j * nx + i] == 0.0) continue;                   if (jz[j * nx + i] == 0.0) continue;
                   sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH                   sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
                   sum2    += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH                   sum2    += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to TOTUSJH
                   sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));              err     += (jz_err[j * nx + i]*jz_err[j * nx + i]*bz[j * nx + i]*bz[j * nx + i]) + (bz_err[j * nx + i]*bz_err[j * nx + i]*jz[j * nx + i]*jz[j * nx + i]);
                   count_mask++;                   count_mask++;
                 }                 }
          }          }
Line 764  int computeHelicity(float *jz_err, float
Line 811  int computeHelicity(float *jz_err, float
         *total_us_ih_ptr      = sum2           ; /* Units are G^2 / m ; keyword is TOTUSJH */         *total_us_ih_ptr      = sum2           ; /* Units are G^2 / m ; keyword is TOTUSJH */
         *total_abs_ih_ptr     = fabs(sum)      ; /* Units are G^2 / m ; keyword is ABSNJZH */         *total_abs_ih_ptr     = fabs(sum)      ; /* Units are G^2 / m ; keyword is ABSNJZH */
  
         *mean_ih_err_ptr      = (sqrt(sum_err*sum_err)) / (count_mask*100.0)    ;  // error in the quantity MEANJZH      *mean_ih_err_ptr      = (sqrt(err)/count_mask)*(1/cdelt1)*(rsun_obs/rsun_ref) ; // error in the quantity MEANJZH
         *total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity TOTUSJH      *total_us_ih_err_ptr  = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity TOTUSJH
         *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity ABSNJZH      *total_abs_ih_err_ptr = (sqrt(err))*(1/cdelt1)*(rsun_obs/rsun_ref) ;            // error in the quantity ABSNJZH
  
         //printf("MEANJZH=%f\n",*mean_ih_ptr);         //printf("MEANJZH=%f\n",*mean_ih_ptr);
         //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);         //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
Line 784  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.  //  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 821  int computeSumAbsPerPolarity(float *jz_e
Line 868  int computeSumAbsPerPolarity(float *jz_e
               }               }
           }           }
  
         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */      *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are Amperes per arcsecond */
         *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));         *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
         //printf("SAVNCPP=%g\n",*totaljzptr);         //printf("SAVNCPP=%g\n",*totaljzptr);
         //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);         //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
Line 867  int computeFreeEnergy(float *bx_err, flo
Line 914  int computeFreeEnergy(float *bx_err, flo
                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if isnan(bx[j * nx + i]) continue;                  if isnan(bx[j * nx + i]) continue;
                  if isnan(by[j * nx + i]) continue;                  if isnan(by[j * nx + i]) continue;
                  sum  += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);              sum  += ( ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
                  sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );              sum1 += (  ((bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])) + ((by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])) );
                  err  += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);              err  += 4.0*(bx[j * nx + i] - bpx[j * nx + i])*(bx[j * nx + i] - bpx[j * nx + i])*(bx_err[j * nx + i]*bx_err[j * nx + i]) +
               4.0*(by[j * nx + i] - bpy[j * nx + i])*(by[j * nx + i] - bpy[j * nx + i])*(by_err[j * nx + i]*by_err[j * nx + i]);
                  count_mask++;                  count_mask++;
               }               }
           }           }
  
         *meanpotptr      = (sum1/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */      /* Units of meanpotptr are ergs per centimeter */
         *meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)          *meanpotptr      = (sum1) / (count_mask*8.*PI) ;     /* Units are ergs per cubic centimeter */
       *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
  
         /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */         /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
         *totpotptr       = (sum)/(8.*PI);         *totpotptr       = (sum)/(8.*PI);
Line 893  int computeFreeEnergy(float *bx_err, flo
Line 942  int computeFreeEnergy(float *bx_err, flo
 /*===========================================*/ /*===========================================*/
 /* Example function 14:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */ /* Example function 14:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
  
 int computeShearAngle(float *bx_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,  int computeShearAngle(float *bx_err, float *by_err, float *bz_err, float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
                       float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)                       float *meanshear_angleptr, float *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
   
   
 { {
         int nx = dims[0];         int nx = dims[0];
         int ny = dims[1];         int ny = dims[1];
         int i = 0;         int i = 0;
         int j = 0;         int j = 0;
         int count_mask = 0;      float count_mask = 0;
       float count = 0;
         double dotproduct = 0.0;         double dotproduct = 0.0;
         double magnitude_potential = 0.0;         double magnitude_potential = 0.0;
         double magnitude_vector = 0.0;         double magnitude_vector = 0.0;
         double shear_angle = 0.0;         double shear_angle = 0.0;
       double denominator = 0.0;
       double term1 = 0.0;
       double term2 = 0.0;
       double term3 = 0.0;
       double sumsum = 0.0;
         double err = 0.0;         double err = 0.0;
         double sum = 0.0;      double part1 = 0.0;
         double count = 0.0;      double part2 = 0.0;
       double part3 = 0.0;
         *area_w_shear_gt_45ptr = 0.0;         *area_w_shear_gt_45ptr = 0.0;
         *meanshear_angleptr = 0.0;         *meanshear_angleptr = 0.0;
  
Line 924  int computeShearAngle(float *bx_err, flo
Line 982  int computeShearAngle(float *bx_err, flo
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                  if isnan(bx[j * nx + i]) continue;                  if isnan(bx[j * nx + i]) continue;
                  if isnan(by[j * nx + i]) continue;                  if isnan(by[j * nx + i]) continue;
                  /* For mean 3D shear angle, area with shear greater than 45*/              if isnan(bx_err[j * nx + i]) continue;
               if isnan(by_err[j * nx + i]) continue;
               if isnan(bz_err[j * nx + i]) continue;
   
               /* For mean 3D shear angle, percentage with shear greater than 45*/
                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );
               //printf("dotproduct=%f\n",dotproduct);
               //printf("magnitude_potential=%f\n",magnitude_potential);
               //printf("magnitude_vector=%f\n",magnitude_vector);
   
                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
               sumsum                  += shear_angle;
               //printf("shear_angle=%f\n",shear_angle);
                  count ++;                  count ++;
                  sum += shear_angle ;  
                  err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));  
                  if (shear_angle > 45) count_mask ++;                  if (shear_angle > 45) count_mask ++;
   
               // For the error analysis
   
               term1 = bx[j * nx + i]*by[j * nx + i]*bpy[j * nx + i] - by[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bz[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bz[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i];
               term2 = bx[j * nx + i]*bx[j * nx + i]*bpy[j * nx + i] - bx[j * nx + i]*by[j * nx + i]*bpx[j * nx + i] + bx[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i] - bz[j * nx + i]*by[j * nx + i]*bpz[j * nx + i];
               term3 = bx[j * nx + i]*bx[j * nx + i]*bpz[j * nx + i] - bx[j * nx + i]*bz[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*by[j * nx + i]*bpz[j * nx + i] - by[j * nx + i]*bz[j * nx + i]*bpy[j * nx + i];
   
               part1 = bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i];
               part2 = bpx[j * nx + i]*bpx[j * nx + i] + bpy[j * nx + i]*bpy[j * nx + i] + bpz[j * nx + i]*bpz[j * nx + i];
               part3 = bx[j * nx + i]*bpx[j * nx + i] + by[j * nx + i]*bpy[j * nx + i] + bz[j * nx + i]*bpz[j * nx + i];
   
               // denominator is squared
               denominator = part1*part1*part1*part2*(1.0-((part3*part3)/(part1*part2)));
   
               err = (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
               (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) +
               (term1*term1*bx_err[j * nx + i]*bx_err[j * nx + i])/(denominator) ;
   
               }               }
           }           }
   
         /* For mean 3D shear angle, area with shear greater than 45*/         /* For mean 3D shear angle, area with shear greater than 45*/
         *meanshear_angleptr = (sum)/(count);                 /* Units are degrees */      *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */
         *meanshear_angle_err_ptr = (sqrt(err*err))/(count);  // error in the quantity (sum)/(count_mask)      *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
         *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */  
       /* 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);
  
         //printf("MEANSHR=%f\n",*meanshear_angleptr);         //printf("MEANSHR=%f\n",*meanshear_angleptr);
         //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);         //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
       //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
  
         return 0;         return 0;
 } }
  
   /*===========================================*/
   /* Example function 15: R parameter as defined in Schrijver, 2007 */
   //
   // Note that there is a restriction on the function fsample()
   // If the following occurs:
   //      nx_out > floor((ny_in-1)/scale + 1)
   //      ny_out > floor((ny_in-1)/scale + 1),
   // where n*_out are the dimensions of the output array and n*_in
   // are the dimensions of the input array, fsample() will usually result
   // in a segfault (though not always, depending on how the segfault was accessed.)
   
   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,
                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, index1;
       double sum = 0.0;
       double err = 0.0;
       *Rparam = 0.0;
       struct fresize_struct fresboxcar, fresgauss;
       struct fint_struct fints;
       float sigma = 10.0/2.3548;
   
       // set up convolution kernels
       init_fresize_boxcar(&fresboxcar,1,1);
       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++)
           {
               index = j * nx1 + i;
               if (rim[index] > 150)
                   p1p0[index]=1.0;
               else
                   p1p0[index]=0.0;
               if (rim[index] < -150)
                   p1n0[index]=1.0;
               else
                   p1n0[index]=0.0;
           }
       }
   
       // =============== [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.0) && (p1n[index] > 0.0))
                   p1[index]=1.0;
               else
                   p1[index]=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;
               if isnan(pmapn[index]) continue;
               if isnan(rim[index]) continue;
               sum += pmapn[index]*abs(rim[index]);
           }
       }
   
       if (sum < 1.0)
           *Rparam = 0.0;
       else
           *Rparam = log10(sum);
   
       //printf("R_VALUE=%f\n",*Rparam);
   
       free_fresize(&fresboxcar);
       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;
   
       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] * bz[i] * k_h;
          fy[i]  = by[i] * bz[i] * 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;
   
       //printf("TOTBSQ=%f\n",*totbsq_ptr);
   
       return 0;
   
   }
  
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/
  
Line 1066  void greenpot(float *bx, float *by, floa
Line 1352  void greenpot(float *bx, float *by, floa
  
 char *sw_functions_version() // Returns CVS version of sw_functions.c char *sw_functions_version() // Returns CVS version of sw_functions.c
 { {
   return strdup("$Id$");      return strdup("$Id");
 } }
  
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

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