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

version 1.9, 2013/05/20 23:55:32 version 1.19, 2013/10/18 23:36:02
Line 63 
Line 63 
 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
 //  =(1.30501e15)Gauss*cm^2 //  =(1.30501e15)Gauss*cm^2
  
 //  The disambig mask value selects only the pixels with values of 5 or 7 -- that is,  
 //  5: pixels for which the radial acute disambiguation solution was chosen  
 //  7: pixels for which the radial acute and NRWA disambiguation agree  
   
 int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux, int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
                    float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,                    float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
                    int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                    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];
     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 89  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 98  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 117  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 153  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 166  int computeGamma(float *bz_err, float *b
Line 175  int computeGamma(float *bz_err, float *b
                     if isnan(bz_err[j * nx + i]) continue;                     if isnan(bz_err[j * nx + i]) continue;
                     if isnan(bh_err[j * nx + i]) continue;                     if isnan(bh_err[j * nx + i]) continue;
                     if (bz[j * nx + i] == 0) continue;                     if (bz[j * nx + i] == 0) continue;
                     sum += (atan(fabs(bz[j * nx + i]/bh[j * nx + i] )))*(180./PI);                      sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
                     err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i])))  * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);                     err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i])))  * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);
                     count_mask++;                     count_mask++;
                   }                   }
Line 174  int computeGamma(float *bz_err, float *b
Line 183  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 196  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 224  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 = 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 332  int computeBtotalderivative(float *bt, i
Line 295  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 344  int computeBtotalderivative(float *bt, i
Line 307  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 = 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 463  int computeBhderivative(float *bh, float
Line 378  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 475  int computeBhderivative(float *bh, float
Line 390  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;
           double sum = 0.0;
           double 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 612  int computeBzderivative(float *bz, float
Line 470  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;
 } }
  
 /*===========================================*/ /*===========================================*/
   
 /* 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 484  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 513  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;
   
         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 = 1; i <= nx-2; i++)
         for (i = 0; i <= nx-2; i++)  
           {           {
              if isnan(bx[j * nx + i]) continue;              for (j = 1; j <= ny-2; j++)
              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;                 // calculate jz at all points
              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 (j = 0; j <= ny-1; j++)  
             {  
                /* calculate jz at all points */  
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);       /* jz is in units of Gauss/pix */  
   
   
                /* 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[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)]) ) ;
                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;
 } }
  
 /*===========================================*/ /*===========================================*/
  
  
 /* 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 605  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++)
           {           {
             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_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 879  int computeJzsmooth(float *bx, float *by
Line 654  int computeJzsmooth(float *bx, float *by
 /* Example function 10:  Twist Parameter, alpha */ /* Example function 10:  Twist Parameter, alpha */
  
 // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
 // for alpha is calculated in the following way (different from Leka and Barnes' approach):  // for alpha is weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
  
        // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz         // numerator   = sum of all Jz*Bz
        // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz         // denominator = sum of Bz*Bz
        // avg alpha = avg Jz / avg Bz         // alpha       = numerator/denominator
   
 // The sign is assigned as follows:  
 // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.  
 // If this value is > 0, then alpha is > 0.  
 // If this value is < 0, then alpha is <0.  
 //  
 // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.  
 // If this value is > 0, then alpha is < 0.  
 // If this value is < 0, then alpha is > 0.  
  
 // The units of alpha are in 1/Mm // The units of alpha are in 1/Mm
 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
Line 904  int computeJzsmooth(float *bx, float *by
Line 670  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;
           double alpha_total         = 0.0;
           double C                   = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
           double total               = 0.0;
           double A                   = 0.0;
           double B                   = 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[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 (jz_err[j * nx + i] > abs(jz[j * nx + i]) ) continue;
                 if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i] ); b++;                  //if (bz_err[j * nx + i] > abs(bz[j * nx + i]) ) continue;
                 //if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);                  A += jz[j*nx+i]*bz[j*nx+i];
                 //if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);                  B += bz[j*nx+i]*bz[j*nx+i];
                 if (bz[j * nx + i] >  0) sum3 += ( jz[j * nx + i] ); c++;  
                 if (bz[j * nx + i] <= 0) sum4 += ( jz[j * nx + i] ); d++;  
                 sum5    += bz[j * nx + i];  
                 /* sum_err is a fractional uncertainty */  
                 sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs( ( (jz[j * nx + i]) / (bz[j * nx + i]) ) *(1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));  
                 count_mask++;  
               }               }
           }           }
  
         sum     = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)); /* the units for (jz/bz) are 1/Mm */          for (i = 1; i < nx-1; i++)
             {
         /* Determine the sign of alpha */              for (j = 1; j < ny-1; j++)
         if ((sum5 > 0) && (sum3 >  0)) sum=sum;                {
         if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
         if ((sum5 < 0) && (sum4 <= 0)) sum=sum;                  if isnan(jz[j * nx + i])   continue;
         if ((sum5 < 0) && (sum4 >  0)) sum=-sum;                  if isnan(bz[j * nx + i])   continue;
                   if (jz[j * nx + i] == 0.0) continue;
         *mean_alpha_ptr = sum; /* Units are 1/Mm */                  if (bz[j * nx + i] == 0.0) continue;
         *mean_alpha_err_ptr    = (sqrt(sum_err*sum_err)) / ((a+b+c+d)*100.); // error in the quantity (sum)/(count_mask); factor of 100 comes from converting percent                  total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
                 }
         printf("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);          /* Determine the absolute value of alpha. The units for alpha are 1/Mm */
         printf("MEANALP_err=%f\n",*mean_alpha_err_ptr);          alpha_total              = ((A/B)*C);
           *mean_alpha_ptr          = alpha_total;
           *mean_alpha_err_ptr      = (C/B)*(sqrt(total));
  
         return 0;         return 0;
 } }
Line 973  int computeHelicity(float *jz_err, float
Line 733  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 1000  int computeHelicity(float *jz_err, float
Line 764  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 1031  int computeSumAbsPerPolarity(float *jz_e
Line 795  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++)
               {               {
                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                   if isnan(jz[j * nx + i]) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);                 err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
Line 1054  int computeSumAbsPerPolarity(float *jz_e
Line 823  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 1063  int computeSumAbsPerPolarity(float *jz_e
Line 832  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 847  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 1093  int computeFreeEnergy(float *bx_err, flo
Line 867  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);
  
         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 1123  int computeFreeEnergy(float *bx_err, flo
Line 896  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 1158  int computeShearAngle(float *bx_err, flo
Line 939  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 1282  void greenpot(float *bx, float *by, floa
Line 1063  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.19
Karen Tian
 Powered by
ViewCVS 0.9.4