(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.11 and 1.16

version 1.11, 2013/05/30 23:26:02 version 1.16, 2013/07/18 00:06:57
Line 73  int computeAbsFlux(float *bz_err, float
Line 73  int computeAbsFlux(float *bz_err, float
  
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
     double sum,err=0.0;      int i = 0;
       int j = 0;
     if (nx <= 0 || ny <= 0) return 1;      int count_mask = 0;
       double sum = 0.0;
       double err = 0.0;
     *absFlux = 0.0;     *absFlux = 0.0;
     *mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;
  
   
       if (nx <= 0 || ny <= 0) return 1;
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
         {         {
                 for (j = 0; j < ny; j++)                 for (j = 0; j < ny; j++)
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 97  int computeAbsFlux(float *bz_err, float
Line 102  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("CMASK=%g\n",*count_mask_ptr);       //printf("cdelt1=%f\n",cdelt1);
      printf("USFLUX=%g\n",*mean_vf_ptr);       //printf("rsun_obs=%f\n",rsun_obs);
      printf("sum=%f\n",sum);       //printf("rsun_ref=%f\n",rsun_ref);
      printf("USFLUX_err=%g\n",*mean_vf_err_ptr);       //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 113  int computeBh(float *bx_err, float *by_e
Line 121  int computeBh(float *bx_err, float *by_e
  
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
     float sum=0.0;      int i = 0;
       int j = 0;
       int count_mask = 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 146  int computeBh(float *bx_err, float *by_e
Line 157  int computeBh(float *bx_err, float *by_e
 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)
 { {
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
     if (nx <= 0 || ny <= 0) return 1;      int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       double err = 0.0;
     *mean_gamma_ptr=0.0;     *mean_gamma_ptr=0.0;
     float sum,err,err_value=0.0;  
  
       if (nx <= 0 || ny <= 0) return 1;
  
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
Line 174  int computeGamma(float *bz_err, float *b
Line 187  int computeGamma(float *bz_err, float *b
           }           }
  
      *mean_gamma_ptr = sum/count_mask;      *mean_gamma_ptr = sum/count_mask;
      *mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.); // error in the quantity (sum)/(count_mask)       *mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.0); // error in the quantity (sum)/(count_mask)
      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 187  int computeGamma(float *bz_err, float *b
Line 200  int computeGamma(float *bz_err, float *b
 int computeB_total(float *bx_err, float *by_err, float *bz_err, float *bt_err, float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask) int computeB_total(float *bx_err, float *by_err, float *bz_err, float *bt_err, float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask)
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
Line 212  int computeB_total(float *bx_err, float
Line 228  int computeB_total(float *bx_err, float
 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)
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
     if (nx <= 0 || ny <= 0) return 1;      int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       double err = 0.0;
     *mean_derivative_btotal_ptr = 0.0;     *mean_derivative_btotal_ptr = 0.0;
     float sum, err = 0.0;  
  
       if (nx <= 0 || ny <= 0) return 1;
  
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
Line 281  int computeBtotalderivative(float *bt, i
Line 299  int computeBtotalderivative(float *bt, i
  
         *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); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
         *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); // error in the quantity (sum)/(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 293  int computeBtotalderivative(float *bt, i
Line 311  int computeBtotalderivative(float *bt, i
 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)
 { {
  
         int nx = dims[0], ny = dims[1];       int nx = dims[0];
         int i, j, count_mask=0;       int ny = dims[1];
        int i = 0;
        int j = 0;
        int count_mask = 0;
        double sum= 0.0;
        double err =0.0;
        *mean_derivative_bh_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_derivative_bh_ptr = 0.0;  
         float sum,err = 0.0;  
   
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
Line 361  int computeBhderivative(float *bh, float
Line 382  int computeBhderivative(float *bh, float
  
         *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 373  int computeBhderivative(float *bh, float
Line 394  int computeBhderivative(float *bh, float
 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)
 { {
  
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double err = 0.0;
           *mean_derivative_bz_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_derivative_bz_ptr = 0.0;  
         float sum,err = 0.0;  
   
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
Line 450  int computeBzderivative(float *bz, float
Line 474  int computeBzderivative(float *bz, float
  
         *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 498  int computeJz(float *bx_err, float *by_e
Line 522  int computeJz(float *bx_err, float *by_e
  
  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;  
   
  
         /* Calculate the derivative*/         /* Calculate the derivative*/
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
  
   
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  //if isnan(by[j * nx + i]) continue;                   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;
               }               }
           }           }
Line 522  int computeJz(float *bx_err, float *by_e
Line 546  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;                   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;
               }               }
           }           }
Line 563  int computeJz(float *bx_err, float *by_e
Line 587  int computeJz(float *bx_err, float *by_e
             {             {
                // 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
   
                //printf("jz[j * nx + i]=%f,i=%d,j=%d\n,",jz[j * nx + i],i,j); // tmp.txt  
                //printf("i=%d,j=%d,jz[j * nx + i]=%f\n,",i,j,jz[j * nx + i]); // tmp1.txt  
   
   
                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]);
Line 590  int computeJzsmooth(float *bx, float *by
Line 609  int computeJzsmooth(float *bx, float *by
  
 { {
  
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double curl = 0.0;
           double us_i = 0.0;
           double err = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         float curl,us_i,test_perimeter,mean_curl,err=0.0;  
   
   
         /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions*/         /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions*/
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
Line 614  int computeJzsmooth(float *bx, float *by
Line 636  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))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // 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))*fabs((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 661  int computeJzsmooth(float *bx, float *by
Line 683  int computeJzsmooth(float *bx, float *by
 int computeAlpha(float *jz_err, float *bz_err, float *bz, int *dims, float *jz, float *jz_smooth, float *mean_alpha_ptr, float *mean_alpha_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs) int computeAlpha(float *jz_err, float *bz_err, float *bz, int *dims, float *jz, float *jz_smooth, float *mean_alpha_ptr, float *mean_alpha_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask, a,b,c,d=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double a = 0.0;
           double b = 0.0;
           double c = 0.0;
           double d = 0.0;
           double sum1 = 0.0;
           double sum2 = 0.0;
           double sum3 = 0.0;
           double sum4 = 0.0;
           double sum = 0.0;
           double sum5 = 0.0;
           double sum6 = 0.0;
           double sum_err = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6, sum_err=0.0;  
   
         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++)
               {               {
                 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_smooth[j * nx + i]) 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 (jz_smooth[j * nx + i] == 0) 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_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++;                 if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i] ); a++;
                 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;                 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;
                 //if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);  
                 //if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);  
                 if (bz[j * nx + i] >  0) sum3 += ( jz[j * nx + i] ); c++;                 if (bz[j * nx + i] >  0) sum3 += ( jz[j * nx + i] ); c++;
                 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;                 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;
                 sum5    += bz[j * nx + i];                 sum5    += bz[j * nx + i];
Line 702  int computeAlpha(float *jz_err, float *b
Line 733  int computeAlpha(float *jz_err, float *b
         if ((sum5 < 0) && (sum4 >  0)) sum=-sum;         if ((sum5 < 0) && (sum4 >  0)) sum=-sum;
  
         *mean_alpha_ptr = sum; /* Units are 1/Mm */         *mean_alpha_ptr = sum; /* Units are 1/Mm */
         *mean_alpha_err_ptr    = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent          *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
   
         printf("a=%d\n",a);  
         printf("b=%d\n",b);  
         printf("d=%d\n",d);  
         printf("c=%d\n",c);  
  
         printf("MEANALP=%f\n",*mean_alpha_ptr);          //printf("MEANALP=%f\n",*mean_alpha_ptr);
         printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);          //printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);
  
         return 0;         return 0;
 } }
Line 730  int computeHelicity(float *jz_err, float
Line 756  int computeHelicity(float *jz_err, float
  
 { {
  
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double sum2 = 0.0;
           double sum_err = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         float sum,sum2,sum_err=0.0;  
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
         {         {
                 for (j = 0; j < ny; j++)                 for (j = 0; j < ny; j++)
Line 757  int computeHelicity(float *jz_err, float
Line 787  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.)    ;  // error in the quantity MEANJZH          *mean_ih_err_ptr      = (sqrt(sum_err*sum_err)) / (count_mask*100.0)    ;  // error in the quantity MEANJZH
         *total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity TOTUSJH          *total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity TOTUSJH
         *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.)               ;  // error in the quantity ABSNJZH          *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0)               ;  // 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 788  int computeSumAbsPerPolarity(float *jz_e
Line 818  int computeSumAbsPerPolarity(float *jz_e
                                                          int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                                                          int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i=0;
           int j=0;
           int count_mask=0;
           double sum1=0.0;
           double sum2=0.0;
           double err=0.0;
           *totaljzptr=0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *totaljzptr=0.0;  
         float sum1,sum2,err=0.0;  
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
Line 811  int computeSumAbsPerPolarity(float *jz_e
Line 845  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 835  int computeFreeEnergy(float *bx_err, flo
Line 869  int computeFreeEnergy(float *bx_err, flo
                                           float cdelt1, double rsun_ref, double rsun_obs)                                           float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
         if (nx <= 0 || ny <= 0) return 1;          int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double sum1 = 0.0;
           double err = 0.0;
         *totpotptr=0.0;         *totpotptr=0.0;
         *meanpotptr=0.0;         *meanpotptr=0.0;
         float sum,err=0.0;  
           if (nx <= 0 || ny <= 0) return 1;
  
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
Line 852  int computeFreeEnergy(float *bx_err, flo
Line 890  int computeFreeEnergy(float *bx_err, flo
                  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 += ( ((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);
                    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])) );
                  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]*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]);
                  count_mask++;                  count_mask++;
               }               }
           }           }
  
         *meanpotptr      = (sum/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */          *meanpotptr      = (sum1/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */
         *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)          *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)
  
         /* 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 879  int computeFreeEnergy(float *bx_err, flo
Line 918  int computeFreeEnergy(float *bx_err, flo
 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 *bh_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], ny = dims[1];          int nx = dims[0];
         int i, j;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double dotproduct = 0.0;
           double magnitude_potential = 0.0;
           double magnitude_vector = 0.0;
           double shear_angle = 0.0;
           double err = 0.0;
           double sum = 0.0;
           double count = 0.0;
           *area_w_shear_gt_45ptr = 0.0;
           *meanshear_angleptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         //*area_w_shear_gt_45ptr=0.0;  
         //*meanshear_angleptr=0.0;  
         float dotproduct, magnitude_potential, magnitude_vector, shear_angle,err=0.0, sum = 0.0, count=0.0, count_mask=0.0;  
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
Line 914  int computeShearAngle(float *bx_err, flo
Line 961  int computeShearAngle(float *bx_err, flo
         /* For mean 3D shear angle, area with shear greater than 45*/         /* For mean 3D shear angle, area with shear greater than 45*/
         *meanshear_angleptr = (sum)/(count);                 /* Units are degrees */         *meanshear_angleptr = (sum)/(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*err))/(count);  // error in the quantity (sum)/(count_mask)
         *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.);/* The area here is a fractional area -- the % of the total area */          *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */
  
         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);
  
         return 0;         return 0;
 } }
Line 1038  void greenpot(float *bx, float *by, floa
Line 1085  void greenpot(float *bx, float *by, floa
  
  
 /*===========END OF KEIJI'S CODE =========================*/ /*===========END OF KEIJI'S CODE =========================*/
   
   char *sw_functions_version() // Returns CVS version of sw_functions.c
   {
     return strdup("$Id$");
   }
   
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/


Legend:
Removed from v.1.11  
changed lines
  Added in v.1.16

Karen Tian
Powered by
ViewCVS 0.9.4