(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.14 and 1.28

version 1.14, 2013/07/04 02:16:52 version 1.28, 2014/03/13 18:40:43
Line 20 
Line 20 
  
    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 41 
  
    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 63 
 //  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  
   
 //  The disambig mask value selects only the pixels with values of 5 or 7 -- that is,  
 //  5: pixels for which the radial acute disambiguation solution was chosen  
 //  7: pixels for which the radial acute and NRWA disambiguation agree  
  
 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 78  int computeAbsFlux(float *bz_err, float
Line 76  int computeAbsFlux(float *bz_err, float
     int i=0;     int i=0;
     int j=0;     int j=0;
     int count_mask=0;     int count_mask=0;
     float sum=0.0;      double sum = 0.0;
     float err=0.0;      double err = 0.0;
     *absFlux = 0.0;     *absFlux = 0.0;
     *mean_vf_ptr = 0.0;     *mean_vf_ptr = 0.0;
  
Line 93  int computeAbsFlux(float *bz_err, float
Line 91  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 102  int computeAbsFlux(float *bz_err, float
Line 99  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 126  int computeBh(float *bx_err, float *by_e
Line 116  int computeBh(float *bx_err, float *by_e
     int i=0;     int i=0;
     int j=0;     int j=0;
     int count_mask=0;     int count_mask=0;
     float sum=0.0;      double sum = 0.0;
     *mean_hf_ptr = 0.0;     *mean_hf_ptr = 0.0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
Line 135  int computeBh(float *bx_err, float *by_e
Line 125  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 152  int computeBh(float *bx_err, float *by_e
Line 152  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 162  int computeGamma(float *bz_err, float *b
Line 164  int computeGamma(float *bz_err, float *b
     int i=0;     int i=0;
     int j=0;     int j=0;
     int count_mask=0;     int count_mask=0;
     float sum=0.0;      double sum = 0.0;
     float err=0.0;      double err = 0.0;
     float err_value=0.0;  
     *mean_gamma_ptr=0.0;     *mean_gamma_ptr=0.0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
Line 179  int computeGamma(float *bz_err, float *b
Line 180  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 += (atan(fabs(bz[j * nx + i]/bh[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 213  int computeB_total(float *bx_err, float
Line 217  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 234  int computeBtotalderivative(float *bt, i
Line 253  int computeBtotalderivative(float *bt, i
     int i=0;     int i=0;
     int j=0;     int j=0;
     int count_mask=0;     int count_mask=0;
     float sum=0.0;      double sum = 0.0;
     float err = 0.0;      double err = 0.0;
     *mean_derivative_btotal_ptr = 0.0;     *mean_derivative_btotal_ptr = 0.0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
Line 285  int computeBtotalderivative(float *bt, i
Line 304  int computeBtotalderivative(float *bt, i
           }           }
  
  
         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 += (((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])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[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)])) / (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 317  int computeBhderivative(float *bh, float
Line 345  int computeBhderivative(float *bh, float
      int i=0;      int i=0;
      int j=0;      int j=0;
      int count_mask=0;      int count_mask=0;
      float sum=0.0;      double sum= 0.0;
      float err =0.0;      double err =0.0;
      *mean_derivative_bh_ptr = 0.0;      *mean_derivative_bh_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
Line 373  int computeBhderivative(float *bh, float
Line 401  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 += (((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])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[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)])) / (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++;
             }             }
           }           }
  
         *mean_derivative_bh_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram         *mean_derivative_bh_ptr     = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
         *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)         *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
         printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);      //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
         printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);      //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
  
         return 0;         return 0;
 } }
Line 400  int computeBzderivative(float *bz, float
Line 436  int computeBzderivative(float *bz, float
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float sum = 0.0;          double sum = 0.0;
         float err = 0.0;      double err = 0.0;
         *mean_derivative_bz_ptr = 0.0;         *mean_derivative_bz_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
Line 461  int computeBzderivative(float *bz, float
Line 497  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 += (((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])) /
               (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[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)])) /
               (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++;
             }             }
           }           }
  
         *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_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
         *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)         *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
         printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);      //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
         printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);      //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
  
         return 0;         return 0;
 } }
Line 528  int computeJz(float *bx_err, float *by_e
Line 571  int computeJz(float *bx_err, float *by_e
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float curl=0.0;  
         float us_i=0.0;  
         float test_perimeter=0.0;  
         float mean_curl=0.0;  
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
Line 585  int computeJz(float *bx_err, float *by_e
Line 624  int computeJz(float *bx_err, float *by_e
              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 = 0; j <= ny-1; j++)          for (j = 1; j <= ny-2; 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( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * 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)]) ) ;                                             (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;
 } }
  
 /*===========================================*/ /*===========================================*/
  
   
 /* 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 619  int computeJzsmooth(float *bx, float *by
Line 657  int computeJzsmooth(float *bx, float *by
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float curl=0.0;          double curl = 0.0;
         float us_i=0.0;      double us_i = 0.0;
         float test_perimeter=0.0;      double err = 0.0;
         float mean_curl=0.0;  
         float err=0.0;  
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
Line 643  int computeJzsmooth(float *bx, float *by
Line 679  int computeJzsmooth(float *bx, float *by
             }             }
           }           }
  
         /* Calculate mean vertical current density (mean_curl) 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);
  
         printf("TOTUSJZ=%g\n",*us_i_ptr);      //printf("TOTUSJZ=%g\n",*us_i_ptr);
         printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);      //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
  
         return 0;         return 0;
  
Line 665  int computeJzsmooth(float *bx, float *by
Line 701  int computeJzsmooth(float *bx, float *by
 /* Example function 10:  Twist Parameter, alpha */ /* Example function 10:  Twist Parameter, alpha */
  
 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
 // for alpha is calculated in the following way (different from Leka and Barnes' approach):  // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
  
        // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz  // numerator   = sum of all Jz*Bz
        // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz  // denominator = sum of Bz*Bz
        // avg alpha = avg Jz / avg Bz  // alpha       = numerator/denominator
   
 // The sign is assigned as follows:  
 // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.  
 // If this value is > 0, then alpha is > 0.  
 // If this value is < 0, then alpha is <0.  
 //  
 // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.  
 // If this value is > 0, then alpha is < 0.  
 // If this value is < 0, then alpha is > 0.  
  
 // The units of alpha are in 1/Mm // The units of alpha are in 1/Mm
 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
Line 694  int computeAlpha(float *jz_err, float *b
Line 721  int computeAlpha(float *jz_err, float *b
         int ny = dims[1];         int ny = dims[1];
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;          double alpha_total         = 0.0;
         float a=0.0;      double C                   = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
         float b=0.0;      double total               = 0.0;
         float c=0.0;      double A                   = 0.0;
         float d=0.0;      double B                   = 0.0;
         float bznew=0.0;  
         float alpha2=0.0;  
         float sum1=0.0;  
         float sum2=0.0;  
         float sum3=0.0;  
         float sum4=0.0;  
         float sum=0.0;  
         float sum5=0.0;  
         float sum6=0.0;  
         float sum_err=0.0;  
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
   
         for (i = 1; i < nx-1; i++)         for (i = 1; i < nx-1; i++)
           {           {
             for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
Line 720  int computeAlpha(float *jz_err, float *b
Line 736  int computeAlpha(float *jz_err, float *b
                 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 (jz[j * nx + i]     == 0.0) continue;                 if (jz[j * nx + i]     == 0.0) continue;
                 if (bz_err[j * nx + i] == 0.0) continue;  
                 if (bz[j * nx + i]     == 0.0) continue;                 if (bz[j * nx + i]     == 0.0) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i] ); a++;              A += jz[j*nx+i]*bz[j*nx+i];
                 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;              B += bz[j*nx+i]*bz[j*nx+i];
                 if (bz[j * nx + i] >  0) sum3 += ( jz[j * nx + i] ); c++;  
                 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;  
                 sum5    += bz[j * nx + i];  
                 /* sum_err is a fractional uncertainty */  
                 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)*(1000000.));  
                 count_mask++;  
               }               }
           }           }
  
         sum     = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */          for (i = 1; i < nx-1; i++)
       {
         /* Determine the sign of alpha */              for (j = 1; j < ny-1; j++)
         if ((sum5 > 0) && (sum3 >  0)) sum=sum;          {
         if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;              if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
         if ((sum5 < 0) && (sum4 <= 0)) sum=sum;              if isnan(jz[j * nx + i])   continue;
         if ((sum5 < 0) && (sum4 >  0)) sum=-sum;              if isnan(bz[j * nx + i])   continue;
               if (jz[j * nx + i] == 0.0) continue;
         *mean_alpha_ptr = sum; /* Units are 1/Mm */              if (bz[j * nx + i] == 0.0) continue;
         *mean_alpha_err_ptr    = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.0); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent              total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
           }
       }
  
         printf("MEANALP=%f\n",*mean_alpha_ptr);      /* Determine the absolute value of alpha. The units for alpha are 1/Mm */
         printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);      alpha_total              = ((A/B)*C);
       *mean_alpha_ptr          = alpha_total;
       *mean_alpha_err_ptr      = (C/B)*(sqrt(total));
  
         return 0;         return 0;
 } }
Line 772  int computeHelicity(float *jz_err, float
Line 783  int computeHelicity(float *jz_err, float
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float sum=0.0;          double sum     = 0.0;
         float sum2=0.0;          double sum2    = 0.0;
         float sum_err=0.0;          double err     = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
Line 785  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 798  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);
  
         printf("TOTUSJH=%f\n",*total_us_ih_ptr);      //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
         printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);      //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
  
         printf("ABSNJZH=%f\n",*total_abs_ih_ptr);      //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
         printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);      //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
  
         return 0;         return 0;
 } }
Line 834  int computeSumAbsPerPolarity(float *jz_e
Line 847  int computeSumAbsPerPolarity(float *jz_e
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float sum1=0.0;          double sum1=0.0;
         float sum2=0.0;      double sum2=0.0;
         float err=0.0;      double err=0.0;
         *totaljzptr=0.0;         *totaljzptr=0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
Line 847  int computeSumAbsPerPolarity(float *jz_e
Line 860  int computeSumAbsPerPolarity(float *jz_e
               {               {
                 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;
               if isnan(jz[j * nx + i]) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);                 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
Line 856  int computeSumAbsPerPolarity(float *jz_e
Line 870  int computeSumAbsPerPolarity(float *jz_e
  
         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */
         *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);
  
         return 0;         return 0;
 } }
Line 885  int computeFreeEnergy(float *bx_err, flo
Line 899  int computeFreeEnergy(float *bx_err, flo
         int i=0;         int i=0;
         int j=0;         int j=0;
         int count_mask=0;         int count_mask=0;
         float sum=0.0;          double sum = 0.0;
         float sum1=0.0;      double sum1 = 0.0;
         float err=0.0;      double err = 0.0;
         *totpotptr=0.0;         *totpotptr=0.0;
         *meanpotptr=0.0;         *meanpotptr=0.0;
  
Line 900  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);
         *totpot_err_ptr  = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));         *totpot_err_ptr  = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
  
         printf("MEANPOT=%g\n",*meanpotptr);      //printf("MEANPOT=%g\n",*meanpotptr);
         printf("MEANPOT_err=%g\n",*meanpot_err_ptr);      //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
  
         printf("TOTPOT=%g\n",*totpotptr);      //printf("TOTPOT=%g\n",*totpotptr);
         printf("TOTPOT_err=%g\n",*totpot_err_ptr);      //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
  
         return 0;         return 0;
 } }
Line 926  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 dotproduct = 0.0;      float count = 0;
         float magnitude_potential = 0.0;      double dotproduct = 0.0;
         float magnitude_vector=0.0;      double magnitude_potential = 0.0;
         float shear_angle=0.0;      double magnitude_vector = 0.0;
         float err=0.0;      double shear_angle = 0.0;
         float sum = 0.0;      double denominator = 0.0;
         float count=0.0;      double term1 = 0.0;
       double term2 = 0.0;
       double term3 = 0.0;
       double sumsum = 0.0;
       double err = 0.0;
       double part1 = 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 957  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_err=%f\n",*meanshear_angle_err_ptr);
       //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
   
           return 0;
   }
   
   /*===========================================*/
   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 nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int index;
       double sum = 0.0;
       double err = 0.0;
       *Rparam = 0.0;
       struct fresize_struct fresboxcar, fresgauss;
       int scale = round(2.0/cdelt1);
       float sigma = 10.0/2.3548;
   
       init_fresize_boxcar(&fresboxcar,1,1);
   
       // set up convolution kernel
       init_fresize_gaussian(&fresgauss,sigma,20,1);
   
       fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
       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;
           }
       }
   
       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);
   
       for (i = 0; i < nx1; i++)
       {
           for (j = 0; j < ny1; j++)
           {
               index = j * nx1 + i;
               if (p1p[index] > 0 && p1n[index] > 0)
                   p1[index]=1.0;
               else
                   p1[index]=0.0;
           }
       }
   
       fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
   
       for (i = 0; i < nx1; i++)
       {
           for (j = 0; j < ny1; j++)
           {
               index = j * nx1 + i;
               sum += pmap[index]*abs(rim[index]);
           }
       }
   
       if (sum < 1.0)
           *Rparam = 0.0;
       else
           *Rparam = log10(sum);
  
         printf("MEANSHR=%f\n",*meanshear_angleptr);      free_fresize(&fresboxcar);
         printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);      free_fresize(&fresgauss);
  
         return 0;         return 0;
 } }
Line 1099  void greenpot(float *bx, float *by, floa
Line 1230  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.14  
changed lines
  Added in v.1.28

Karen Tian
Powered by
ViewCVS 0.9.4