(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.9 and 1.14

version 1.9, 2013/05/20 23:55:32 version 1.14, 2013/07/04 02:16:52
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;
       float sum=0.0;
       float 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("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("CMASK=%g\n",*count_mask_ptr);
      printf("USFLUX=%g\n",*mean_vf_ptr);      printf("USFLUX=%g\n",*mean_vf_ptr);
      printf("sum=%f\n",sum);      printf("sum=%f\n",sum);
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];
       int i=0;
       int j=0;
       int count_mask=0;
     float sum=0.0;     float sum=0.0;
     *mean_hf_ptr = 0.0;     *mean_hf_ptr = 0.0;
  
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;
       float sum=0.0;
       float err=0.0;
       float err_value=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 188  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 201  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 229  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;
       float sum=0.0;
       float 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 = 3; i <= nx-4; i++)          for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  derx_bt[j * nx + i] = (-bt[j * nx + (i-3)] + 9.0*bt[j * nx + (i-2)] - 45.0*bt[j * nx + (i-1)] + 45*bt[j * nx + (i+1)] - 9.0*bt[j * nx + (i+2)] + bt[j * nx + (i+3)])*(1.0/60.0);                  derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
               }               }
           }           }
  
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 3; j <= ny-4; j++)              for (j = 1; j <= ny-2; j++)
               {               {
                  dery_bt[j * nx + i] = (-bt[(j-3) * nx + i] + 9.0*bt[(j-2) * nx + i] - 45.0*bt[(j-1) * nx + i] + 45*bt[(j+1) * nx + i] - 9.0*bt[(j+2) * nx + i] + bt[(j+3) * nx + i])*(1.0/60.0);                  dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
               }               }
           }           }
  
  
         /* consider the edges: 3 pixels on each edge, for a total of 12 edge formulae below */          /* consider the edges */
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              derx_bt[j * nx + i] = (-147.0*bt[j * nx + i] + 360.0*bt[j * nx + (i+1)] - 450.0*bt[j * nx + (i+2)] + 400.0*bt[j * nx + (i+3)] - 225.0*bt[j * nx + (i+4)] + 72.0*bt[j * nx + (i+5)] - 10.0*bt[j * nx + (i+6)])*(1.0/60.0);               derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[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++)
           {           {
              derx_bt[j * nx + i] = ( 147.0*bt[j * nx + i] - 360.0*bt[j * nx + (i-1)] + 450.0*bt[j * nx + (i-2)] - 400.0*bt[j * nx + (i-3)] + 225.0*bt[j * nx + (i-4)] - 72.0*bt[j * nx + (i-5)] + 10.0*bt[j * nx + (i-6)])*(1.0/60.0);               derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5;
           }  
   
   
         i=1;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bt[j * nx + i] = (-10.0*bt[j * nx + i] - 77.0*bt[j * nx + (i+1)] + 150.0*bt[j * nx + (i+2)] - 100.0*bt[j * nx + (i+3)] + 50.0*bt[j * nx + (i+4)] - 15.0*bt[j * nx + (i+5)] + 2.0*bt[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bt[j * nx + i] = ( 10.0*bt[j * nx + i] + 77.0*bt[j * nx + (i-1)] - 150.0*bt[j * nx + (i-2)] + 100.0*bt[j * nx + (i-3)] - 50.0*bt[j * nx + (i-4)] + 15.0*bt[j * nx + (i-5)] - 2.0*bt[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         i=2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bt[j * nx + i] = ( 2.0*bt[j * nx + i] - 24.0*bt[j * nx + (i+1)] - 35.0*bt[j * nx + (i+2)] + 80.0*bt[j * nx + (i+3)] - 30.0*bt[j * nx + (i+4)] + 8.0*bt[j * nx + (i+5)] - bt[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-3;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bt[j * nx + i] = (-2.0*bt[j * nx + i] + 24.0*bt[j * nx + (i-1)] + 35.0*bt[j * nx + (i-2)] - 80.0*bt[j * nx + (i-3)] + 30.0*bt[j * nx + (i-4)] - 8.0*bt[j * nx + (i-5)] + bt[j * nx + (i-6)])*(1.0/60.0);  
           }           }
  
   
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              dery_bt[j * nx + i] = (-147.0*bt[j * nx + i] + 360.0*bt[(j+1) * nx + i] - 450.0*bt[(j+2) * nx + i] + 400.0*bt[(j+3) * nx + i] - 225.0*bt[(j+4) * nx + i] + 72.0*bt[(j+5) * nx + i] - 10.0*bt[(j+6) * nx + i])*(1.0/60.0);               dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(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++)
           {           {
              dery_bt[j * nx + i] = ( 147.0*bt[j * nx + i] - 360.0*bt[(j-1) * nx + i] + 450.0*bt[(j-2) * nx + i] - 400.0*bt[(j-3) * nx + i] + 225.0*bt[(j-4) * nx + i] - 72.0*bt[(j-5) * nx + i] + 10.0*bt[(j-6) * nx + i])*(1.0/60.0);               dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
           }  
   
         j=1;  
         for (i = 0; i <= nx-2; i++)  
           {  
              dery_bt[j * nx + i] = (-10.0*bt[j * nx + i] - 77.0*bt[(j+1) * nx + i] + 150.0*bt[(j+2) * nx + i] - 100.0*bt[(j+3) * nx + i] + 50.0*bt[(j+4) * nx + i] - 15.0*bt[(j+5) * nx + i] + 2.0*bt[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-2;  
         for (i = 0; i <= nx-2; i++)  
           {  
              dery_bt[j * nx + i] = ( 10.0*bt[j * nx + i] + 77.0*bt[(j-1) * nx + i] - 150.0*bt[(j-2) * nx + i] + 100.0*bt[(j-3) * nx + i] - 50.0*bt[(j-4) * nx + i] + 15.0*bt[(j-5) * nx + i] - 2.0*bt[(j-6) * nx + i])*(1.0/60.0);  
           }  
   
         j=2;  
         for (i = 0; i <= nx-3; i++)  
           {  
              dery_bt[j * nx + i] = ( 2.0*bt[j * nx + i] - 24.0*bt[(j+1) * nx + i] - 35.0*bt[(j+2) * nx + i] + 80.0*bt[(j+3) * nx + i] - 30.0*bt[(j+4) * nx + i] + 8.0*bt[(j+5) * nx + i] - bt[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-3;  
         for (i = 0; i <= nx-3; i++)  
           {  
              dery_bt[j * nx + i] = (-2.0*bt[j * nx + i] + 24.0*bt[(j-1) * nx + i] + 35.0*bt[(j-2) * nx + i] - 80.0*bt[(j-3) * nx + i] + 30.0*bt[(j-4) * nx + i] - 8.0*bt[(j-5) * nx + i] + bt[(j-6) * nx + i])*(1.0/60.0);  
           }           }
  
  
Line 344  int computeBtotalderivative(float *bt, i
Line 312  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;
        float sum=0.0;
        float 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 = 3; i <= nx-4; i++)          for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  derx_bh[j * nx + i] = (-bh[j * nx + (i-3)] + 9.0*bh[j * nx + (i-2)] - 45.0*bh[j * nx + (i-1)] + 45*bh[j * nx + (i+1)] - 9.0*bh[j * nx + (i+2)] + bh[j * nx + (i+3)])*(1.0/60.0);                  derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
               }               }
           }           }
  
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 3; j <= ny-4; j++)              for (j = 1; j <= ny-2; j++)
               {               {
                  dery_bh[j * nx + i] = (-bh[(j-3) * nx + i] + 9.0*bh[(j-2) * nx + i] - 45.0*bh[(j-1) * nx + i] + 45*bh[(j+1) * nx + i] - 9.0*bh[(j+2) * nx + i] + bh[(j+3) * nx + i])*(1.0/60.0);                  dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
               }               }
           }           }
  
  
         /* consider the edges: 3 pixels on each edge, for a total of 12 edge formulae below */          /* consider the edges */
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              derx_bh[j * nx + i] = (-147.0*bh[j * nx + i] + 360.0*bh[j * nx + (i+1)] - 450.0*bh[j * nx + (i+2)] + 400.0*bh[j * nx + (i+3)] - 225.0*bh[j * nx + (i+4)] + 72.0*bh[j * nx + (i+5)] - 10.0*bh[j * nx + (i+6)])*(1.0/60.0);               derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[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++)
           {           {
              derx_bh[j * nx + i] = ( 147.0*bh[j * nx + i] - 360.0*bh[j * nx + (i-1)] + 450.0*bh[j * nx + (i-2)] - 400.0*bh[j * nx + (i-3)] + 225.0*bh[j * nx + (i-4)] - 72.0*bh[j * nx + (i-5)] + 10.0*bh[j * nx + (i-6)])*(1.0/60.0);               derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5;
           }  
   
   
         i=1;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bh[j * nx + i] = (-10.0*bh[j * nx + i] - 77.0*bh[j * nx + (i+1)] + 150.0*bh[j * nx + (i+2)] - 100.0*bh[j * nx + (i+3)] + 50.0*bh[j * nx + (i+4)] - 15.0*bh[j * nx + (i+5)] + 2.0*bh[j * nx + (i+6)])*(1.0/60.0);  
           }           }
  
         i=nx-2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bh[j * nx + i] = ( 10.0*bh[j * nx + i] + 77.0*bh[j * nx + (i-1)] - 150.0*bh[j * nx + (i-2)] + 100.0*bh[j * nx + (i-3)] - 50.0*bh[j * nx + (i-4)] + 15.0*bh[j * nx + (i-5)] - 2.0*bh[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         i=2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bh[j * nx + i] = ( 2.0*bh[j * nx + i] - 24.0*bh[j * nx + (i+1)] - 35.0*bh[j * nx + (i+2)] + 80.0*bh[j * nx + (i+3)] - 30.0*bh[j * nx + (i+4)] + 8.0*bh[j * nx + (i+5)] - bh[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-3;  
         for (j = 0; j <= ny-2; j++)  
           {  
              derx_bh[j * nx + i] = (-2.0*bh[j * nx + i] + 24.0*bh[j * nx + (i-1)] + 35.0*bh[j * nx + (i-2)] - 80.0*bh[j * nx + (i-3)] + 30.0*bh[j * nx + (i-4)] - 8.0*bh[j * nx + (i-5)] + bh[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              dery_bh[j * nx + i] = (-147.0*bh[j * nx + i] + 360.0*bh[(j+1) * nx + i] - 450.0*bh[(j+2) * nx + i] + 400.0*bh[(j+3) * nx + i] - 225.0*bh[(j+4) * nx + i] + 72.0*bh[(j+5) * nx + i] - 10.0*bh[(j+6) * nx + i])*(1.0/60.0);               dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(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++)
           {           {
              dery_bh[j * nx + i] = ( 147.0*bh[j * nx + i] - 360.0*bh[(j-1) * nx + i] + 450.0*bh[(j-2) * nx + i] - 400.0*bh[(j-3) * nx + i] + 225.0*bh[(j-4) * nx + i] - 72.0*bh[(j-5) * nx + i] + 10.0*bh[(j-6) * nx + i])*(1.0/60.0);               dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
           }  
   
         j=1;  
         for (i = 0; i <= nx-2; i++)  
           {  
              dery_bh[j * nx + i] = (-10.0*bh[j * nx + i] - 77.0*bh[(j+1) * nx + i] + 150.0*bh[(j+2) * nx + i] - 100.0*bh[(j+3) * nx + i] + 50.0*bh[(j+4) * nx + i] - 15.0*bh[(j+5) * nx + i] + 2.0*bh[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-2;  
         for (i = 0; i <= nx-2; i++)  
           {  
              dery_bh[j * nx + i] = ( 10.0*bh[j * nx + i] + 77.0*bh[(j-1) * nx + i] - 150.0*bh[(j-2) * nx + i] + 100.0*bh[(j-3) * nx + i] - 50.0*bh[(j-4) * nx + i] + 15.0*bh[(j-5) * nx + i] - 2.0*bh[(j-6) * nx + i])*(1.0/60.0);  
           }  
   
         j=2;  
         for (i = 0; i <= nx-3; i++)  
           {  
              dery_bh[j * nx + i] = ( 2.0*bh[j * nx + i] - 24.0*bh[(j+1) * nx + i] - 35.0*bh[(j+2) * nx + i] + 80.0*bh[(j+3) * nx + i] - 30.0*bh[(j+4) * nx + i] + 8.0*bh[(j+5) * nx + i] - bh[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-3;  
         for (i = 0; i <= nx-3; i++)  
           {  
              dery_bh[j * nx + i] = (-2.0*bh[j * nx + i] + 24.0*bh[(j-1) * nx + i] + 35.0*bh[(j-2) * nx + i] - 80.0*bh[(j-3) * nx + i] + 30.0*bh[(j-4) * nx + i] - 8.0*bh[(j-5) * nx + i] + bh[(j-6) * nx + i])*(1.0/60.0);  
           }           }
  
  
Line 475  int computeBhderivative(float *bh, float
Line 395  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;
         if (nx <= 0 || ny <= 0) return 1;          int j=0;
           int count_mask=0;
           float sum = 0.0;
           float err = 0.0;
         *mean_derivative_bz_ptr = 0.0;         *mean_derivative_bz_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 = 3; i <= nx-4; i++)          for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                  derx_bz[j * nx + i] = (-bz[j * nx + (i-3)] + 9.0*bz[j * nx + (i-2)] - 45.0*bz[j * nx + (i-1)] + 45*bz[j * nx + (i+1)] - 9.0*bz[j * nx + (i+2)] + bz[j * nx + (i+3)])*(1.0/60.0);                  derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
               }               }
           }           }
  
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 3; j <= ny-4; j++)              for (j = 1; j <= ny-2; j++)
               {               {
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                  dery_bz[j * nx + i] = (-bz[(j-3) * nx + i] + 9.0*bz[(j-2) * nx + i] - 45.0*bz[(j-1) * nx + i] + 45*bz[(j+1) * nx + i] - 9.0*bz[(j+2) * nx + i] + bz[(j+3) * nx + i])*(1.0/60.0);                  dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
               }               }
           }           }
  
  
         /* consider the edges: 3 pixels on each edge, for a total of 12 edge formulae below */          /* consider the edges */
         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;              if isnan(bz[j * nx + i]) continue;
              derx_bz[j * nx + i] = (-147.0*bz[j * nx + i] + 360.0*bz[j * nx + (i+1)] - 450.0*bz[j * nx + (i+2)] + 400.0*bz[j * nx + (i+3)] - 225.0*bz[j * nx + (i+4)] + 72.0*bz[j * nx + (i+5)] - 10.0*bz[j * nx + (i+6)])*(1.0/60.0);               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;              if isnan(bz[j * nx + i]) continue;
              derx_bz[j * nx + i] = ( 147.0*bz[j * nx + i] - 360.0*bz[j * nx + (i-1)] + 450.0*bz[j * nx + (i-2)] - 400.0*bz[j * nx + (i-3)] + 225.0*bz[j * nx + (i-4)] - 72.0*bz[j * nx + (i-5)] + 10.0*bz[j * nx + (i-6)])*(1.0/60.0);               derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
           }           }
  
   
         i=1;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = (-10.0*bz[j * nx + i] - 77.0*bz[j * nx + (i+1)] + 150.0*bz[j * nx + (i+2)] - 100.0*bz[j * nx + (i+3)] + 50.0*bz[j * nx + (i+4)] - 15.0*bz[j * nx + (i+5)] + 2.0*bz[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( 10.0*bz[j * nx + i] + 77.0*bz[j * nx + (i-1)] - 150.0*bz[j * nx + (i-2)] + 100.0*bz[j * nx + (i-3)] - 50.0*bz[j * nx + (i-4)] + 15.0*bz[j * nx + (i-5)] - 2.0*bz[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         i=2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( 2.0*bz[j * nx + i] - 24.0*bz[j * nx + (i+1)] - 35.0*bz[j * nx + (i+2)] + 80.0*bz[j * nx + (i+3)] - 30.0*bz[j * nx + (i+4)] + 8.0*bz[j * nx + (i+5)] - bz[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-3;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = (-2.0*bz[j * nx + i] + 24.0*bz[j * nx + (i-1)] + 35.0*bz[j * nx + (i-2)] - 80.0*bz[j * nx + (i-3)] + 30.0*bz[j * nx + (i-4)] - 8.0*bz[j * nx + (i-5)] + bz[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         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;              if isnan(bz[j * nx + i]) continue;
              dery_bz[j * nx + i] = (-147.0*bz[j * nx + i] + 360.0*bz[(j+1) * nx + i] - 450.0*bz[(j+2) * nx + i] + 400.0*bz[(j+3) * nx + i] - 225.0*bz[(j+4) * nx + i] + 72.0*bz[(j+5) * nx + i] - 10.0*bz[(j+6) * nx + i])*(1.0/60.0);               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;              if isnan(bz[j * nx + i]) continue;
              dery_bz[j * nx + i] = ( 147.0*bz[j * nx + i] - 360.0*bz[(j-1) * nx + i] + 450.0*bz[(j-2) * nx + i] - 400.0*bz[(j-3) * nx + i] + 225.0*bz[(j-4) * nx + i] - 72.0*bz[(j-5) * nx + i] + 10.0*bz[(j-6) * nx + i])*(1.0/60.0);               dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
           }  
   
         j=1;  
         for (i = 0; i <= nx-2; i++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = (-10.0*bz[j * nx + i] - 77.0*bz[(j+1) * nx + i] + 150.0*bz[(j+2) * nx + i] - 100.0*bz[(j+3) * nx + i] + 50.0*bz[(j+4) * nx + i] - 15.0*bz[(j+5) * nx + i] + 2.0*bz[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-2;  
         for (i = 0; i <= nx-2; i++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( 10.0*bz[j * nx + i] + 77.0*bz[(j-1) * nx + i] - 150.0*bz[(j-2) * nx + i] + 100.0*bz[(j-3) * nx + i] - 50.0*bz[(j-4) * nx + i] + 15.0*bz[(j-5) * nx + i] - 2.0*bz[(j-6) * nx + i])*(1.0/60.0);  
           }  
   
         j=2;  
         for (i = 0; i <= nx-3; i++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( 2.0*bz[j * nx + i] - 24.0*bz[(j+1) * nx + i] - 35.0*bz[(j+2) * nx + i] + 80.0*bz[(j+3) * nx + i] - 30.0*bz[(j+4) * nx + i] + 8.0*bz[(j+5) * nx + i] - bz[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-3;  
         for (i = 0; i <= nx-3; i++)  
           {  
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = (-2.0*bz[j * nx + i] + 24.0*bz[(j-1) * nx + i] + 35.0*bz[(j-2) * nx + i] - 80.0*bz[(j-3) * nx + i] + 30.0*bz[(j-4) * nx + i] - 8.0*bz[(j-5) * nx + i] + bz[(j-6) * nx + i])*(1.0/60.0);  
           }           }
  
  
Line 619  int computeBzderivative(float *bz, float
Line 482  int computeBzderivative(float *bz, float
 } }
  
 /*===========================================*/ /*===========================================*/
   
 /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */ /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */
  
 //  In discretized space like data pixels, //  In discretized space like data pixels,
Line 627  int computeBzderivative(float *bz, float
Line 489  int computeBzderivative(float *bz, float
 //  the integration of the field Bx and By along //  the integration of the field Bx and By along
 //  the circumference of the data pixel divided by the area of the pixel. //  the circumference of the data pixel divided by the area of the pixel.
 //  One form of differencing (a word for the differential operator //  One form of differencing (a word for the differential operator
 //  in the discretized space) of the curl is expressed as the following,  //  in the discretized space) of the curl is expressed as
 //  which utilizes a second-order finite difference method:  
   
 //  (dx * (Bx(i,j-1)+Bx(i,j)) / 2 //  (dx * (Bx(i,j-1)+Bx(i,j)) / 2
 //  +dy * (By(i+1,j)+By(i,j)) / 2 //  +dy * (By(i+1,j)+By(i,j)) / 2
 //  -dx * (Bx(i,j+1)+Bx(i,j)) / 2 //  -dx * (Bx(i,j+1)+Bx(i,j)) / 2
 //  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) //  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy)
 // //
 // //
 //  However, for the purposes of this calculation, we will use a sixth-order finite difference  
 //  method taken from the pencil code:  
 //  
 //  dBy/dx = (    -By*(i-3,j) + 9By(i-2,j) - 45By(i-1,j) + 45By(i+1,j) - 9By(i+2,j) + By(i+3,j)    )/ 60  
 //  and similarly for dBx/dy.  
 //  
 //  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), //  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
 //  one must perform the following unit conversions: //  one must perform the following unit conversions:
 //  (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or //  (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
Line 664  int computeBzderivative(float *bz, float
Line 518  int computeBzderivative(float *bz, float
 //            int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx, //            int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *noisebx,
 //              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)
  
  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
         printf("nx=%d\n",nx);          int i=0;
         printf("ny=%d\n",ny);          int j=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;
         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 = 3; i <= nx-4; 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-3)] + 9.0*by[j * nx + (i-2)] - 45.0*by[j * nx + (i-1)] + 45*by[j * nx + (i+1)] - 9.0*by[j * nx + (i+2)] + by[j * nx + (i+3)])*(1.0/60.0);                   derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
               }               }
           }           }
  
         /* brute force method of calculating the derivative (no consideration for edges) */  
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 3; j <= ny-4; 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-3) * nx + i] + 9.0*bx[(j-2) * nx + i] - 45.0*bx[(j-1) * nx + i] + 45*bx[(j+1) * nx + i] - 9.0*bx[(j+2) * nx + i] + bx[(j+3) * nx + i])*(1.0/60.0);                   dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
               }               }
           }           }
  
         /* consider the edges: 3 pixels on each edge, for a total of 12 edge formulae below */          // consider the edges
         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;              if isnan(by[j * nx + i]) continue;
              derx[j * nx + i] = (-147.0*by[j * nx + i] + 360.0*by[j * nx + (i+1)] - 450.0*by[j * nx + (i+2)] + 400.0*by[j * nx + (i+3)] - 225.0*by[j * nx + (i+4)] + 72.0*by[j * nx + (i+5)] - 10.0*by[j * nx + (i+6)])*(1.0/60.0);               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;              if isnan(by[j * nx + i]) continue;
              derx[j * nx + i] = ( 147.0*by[j * nx + i] - 360.0*by[j * nx + (i-1)] + 450.0*by[j * nx + (i-2)] - 400.0*by[j * nx + (i-3)] + 225.0*by[j * nx + (i-4)] - 72.0*by[j * nx + (i-5)] + 10.0*by[j * nx + (i-6)])*(1.0/60.0);               derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
           }           }
  
   
         i=1;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = (-10.0*by[j * nx + i] - 77.0*by[j * nx + (i+1)] + 150.0*by[j * nx + (i+2)] - 100.0*by[j * nx + (i+3)] + 50.0*by[j * nx + (i+4)] - 15.0*by[j * nx + (i+5)] + 2.0*by[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( 10.0*by[j * nx + i] + 77.0*by[j * nx + (i-1)] - 150.0*by[j * nx + (i-2)] + 100.0*by[j * nx + (i-3)] - 50.0*by[j * nx + (i-4)] + 15.0*by[j * nx + (i-5)] - 2.0*by[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         i=2;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( 2.0*by[j * nx + i] - 24.0*by[j * nx + (i+1)] - 35.0*by[j * nx + (i+2)] + 80.0*by[j * nx + (i+3)] - 30.0*by[j * nx + (i+4)] + 8.0*by[j * nx + (i+5)] - by[j * nx + (i+6)])*(1.0/60.0);  
           }  
   
         i=nx-3;  
         for (j = 0; j <= ny-2; j++)  
           {  
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = (-2.0*by[j * nx + i] + 24.0*by[j * nx + (i-1)] + 35.0*by[j * nx + (i-2)] - 80.0*by[j * nx + (i-3)] + 30.0*by[j * nx + (i-4)] - 8.0*by[j * nx + (i-5)] + by[j * nx + (i-6)])*(1.0/60.0);  
           }  
   
   
         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;              if isnan(bx[j * nx + i]) continue;
              dery[j * nx + i] = (-147.0*bx[j * nx + i] + 360.0*bx[(j+1) * nx + i] - 450.0*bx[(j+2) * nx + i] + 400.0*bx[(j+3) * nx + i] - 225.0*bx[(j+4) * nx + i] + 72.0*bx[(j+5) * nx + i] - 10.0*bx[(j+6) * nx + i])*(1.0/60.0);               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;              if isnan(bx[j * nx + i]) continue;
              dery[j * nx + i] = ( 147.0*bx[j * nx + i] - 360.0*bx[(j-1) * nx + i] + 450.0*bx[(j-2) * nx + i] - 400.0*bx[(j-3) * nx + i] + 225.0*bx[(j-4) * nx + i] - 72.0*bx[(j-5) * nx + i] + 10.0*bx[(j-6) * nx + i])*(1.0/60.0);               dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
           }  
   
         j=1;  
         for (i = 0; i <= nx-2; i++)  
           {  
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = (-10.0*bx[j * nx + i] - 77.0*bx[(j+1) * nx + i] + 150.0*bx[(j+2) * nx + i] - 100.0*bx[(j+3) * nx + i] + 50.0*bx[(j+4) * nx + i] - 15.0*bx[(j+5) * nx + i] + 2.0*bx[(j+6) * nx + i])*(1.0/60.0);  
           }           }
  
         j=ny-2;  
         for (i = 0; i <= nx-2; i++)  
           {  
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( 10.0*bx[j * nx + i] + 77.0*bx[(j-1) * nx + i] - 150.0*bx[(j-2) * nx + i] + 100.0*bx[(j-3) * nx + i] - 50.0*bx[(j-4) * nx + i] + 15.0*bx[(j-5) * nx + i] - 2.0*bx[(j-6) * nx + i])*(1.0/60.0);  
           }  
   
         j=2;  
         for (i = 0; i <= nx-3; i++)  
           {  
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( 2.0*bx[j * nx + i] - 24.0*bx[(j+1) * nx + i] - 35.0*bx[(j+2) * nx + i] + 80.0*bx[(j+3) * nx + i] - 30.0*bx[(j+4) * nx + i] + 8.0*bx[(j+5) * nx + i] - bx[(j+6) * nx + i])*(1.0/60.0);  
           }  
   
         j=ny-3;  
         for (i = 0; i <= nx-3; i++)  
           {  
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = (-2.0*bx[j * nx + i] + 24.0*bx[(j-1) * nx + i] + 35.0*bx[(j-2) * nx + i] - 80.0*bx[(j-3) * nx + i] + 30.0*bx[(j-4) * nx + i] - 8.0*bx[(j-5) * nx + i] + bx[(j-6) * nx + i])*(1.0/60.0);  
           }  
  
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; 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]) +
                                               (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
                /* calculate the error in jz at all points */  
                jz_err[j * nx + i]=sqrt( (1.0/60.0)*bx_err[(j-3) * nx + i]*(1.0/60.0)*bx_err[(j-3) * nx + i] + (9.0/60.0)*bx_err[(j-2) * nx + i]*(9.0/60.0)*bx_err[(j-2) * nx + i] + (45.0/60.0)*bx_err[(j-1) * nx + i]*(45.0/60.0)*bx_err[(j-1) * nx + i] +  
                                         (1.0/60.0)*bx_err[(j+3) * nx + i]*(1.0/60.0)*bx_err[(j+3) * nx + i] + (9.0/60.0)*bx_err[(j+2) * nx + i]*(9.0/60.0)*bx_err[(j+2) * nx + i] + (45.0/60.0)*bx_err[(j+1) * nx + i]*(45.0/60.0)*bx_err[(j+1) * nx + i] +  
                                         (1.0/60.0)*by_err[j * nx + (i-3)]*(1.0/60.0)*by_err[j * nx + (i-3)] + (9.0/60.0)*by_err[j * nx + (i-2)]*(9.0/60.0)*by_err[j * nx + (i-2)] + (45.0/60.0)*by_err[j * nx + (i-1)]*(45.0/60.0)*by_err[j * nx + (i-1)] +  
                                         (1.0/60.0)*by_err[j * nx + (i+3)]*(1.0/60.0)*by_err[j * nx + (i+3)] + (9.0/60.0)*by_err[j * nx + (i+2)]*(9.0/60.0)*by_err[j * nx + (i+2)] + (45.0/60.0)*by_err[j * nx + (i+1)]*(45.0/60.0)*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++;
             }             }
             //printf("\n");  
           }           }
  
         return 0;         return 0;
Line 812  int computeJz(float *bx_err, float *by_e
Line 605  int computeJz(float *bx_err, float *by_e
 /*===========================================*/ /*===========================================*/
  
  
 /* Example function 9:  Compute quantities on smoothed Jz array */  /* Example function 9:  Compute quantities on Jz array */
   // Compute mean and total current on Jz array.
 // All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's  
 // fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels  
 // and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values  
 // of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course,  
 // give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels.  
  
 int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth, int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth,
                     float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,                     float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
Line 826  int computeJzsmooth(float *bx, float *by
Line 614  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;
           float curl=0.0;
           float us_i=0.0;
           float test_perimeter=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;
  
         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++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                //printf("%f ",jz_smooth[j * nx + i]);  
                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(derx[j * nx + i]) continue;                if isnan(derx[j * nx + i]) continue;
                if isnan(dery[j * nx + i]) continue;                if isnan(dery[j * nx + i]) continue;
                //if isnan(jz_smooth[j * nx + i]) continue;  
                //curl +=     (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */  
                //us_i += fabs(jz_smooth[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */  
                //jz_rms_err[j * nx + i] = sqrt(jz_err_squared_smooth[j * nx + i]);  
                //err += (jz_rms_err[j * nx + i]*jz_rms_err[j * nx + i]);  
                if isnan(jz[j * nx + i]) continue;                if isnan(jz[j * nx + i]) continue;
                curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
                us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */                us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
                err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);                err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);
                count_mask++;                count_mask++;
             }             }
             //printf("\n");  
           }           }
  
         /* 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_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
Line 904  int computeJzsmooth(float *bx, float *by
Line 690  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;
           float a=0.0;
           float b=0.0;
           float c=0.0;
           float d=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;
  
         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 945  int computeAlpha(float *jz_err, float *b
Line 742  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);
Line 973  int computeHelicity(float *jz_err, float
Line 767  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;
           float sum=0.0;
           float sum2=0.0;
           float 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 1000  int computeHelicity(float *jz_err, float
Line 798  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);
Line 1031  int computeSumAbsPerPolarity(float *jz_e
Line 829  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;
           float sum1=0.0;
           float sum2=0.0;
           float 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 1063  int computeSumAbsPerPolarity(float *jz_e
Line 865  int computeSumAbsPerPolarity(float *jz_e
 /*===========================================*/ /*===========================================*/
 /* Example function 13:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */ /* Example function 13:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
 // automatically yields erg per cubic centimeter for an input B in Gauss.  // automatically yields erg per cubic centimeter for an input B in Gauss. Note that the 8*PI can come out of the integral; thus,
   // the integral is over B^2 dV and the 8*PI is divided at the end.
 // //
 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm: // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
Line 1077  int computeFreeEnergy(float *bx_err, flo
Line 880  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;
           float sum=0.0;
           float sum1=0.0;
           float 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 1093  int computeFreeEnergy(float *bx_err, flo
Line 900  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])) ) / 8.*PI;                   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]);
                  //err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];  
                  count_mask++;                  count_mask++;
               }               }
           }           }
  
         *meanpotptr    = (sum) / (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)
         //*mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // 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*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/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);
Line 1123  int computeFreeEnergy(float *bx_err, flo
Line 929  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;
           float dotproduct = 0.0;
           float magnitude_potential = 0.0;
           float magnitude_vector=0.0;
           float shear_angle=0.0;
           float err=0.0;
           float sum = 0.0;
           float 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 1158  int computeShearAngle(float *bx_err, flo
Line 972  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);
Line 1282  void greenpot(float *bx, float *by, floa
Line 1096  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.9  
changed lines
  Added in v.1.14

Karen Tian
Powered by
ViewCVS 0.9.4