(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.22 and 1.36

version 1.22, 2013/11/11 23:21:21 version 1.36, 2015/10/30 18:53:02
Line 17 
Line 17 
     MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter     MEANPOT Mean photospheric excess magnetic energy density in ergs per cubic centimeter
     TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter     TOTPOT Total photospheric magnetic energy density in ergs per cubic centimeter
     MEANSHR Mean shear angle (measured using Btotal) in degrees     MEANSHR Mean shear angle (measured using Btotal) in degrees
    R_VALUE Karel Schrijver's R parameter
  
    The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and    The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
    pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD    pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
Line 195  int computeGamma(float *bz_err, float *b
Line 196  int computeGamma(float *bz_err, float *b
      *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);      *mean_gamma_err_ptr = (sqrt(err)/(count_mask))*(180./PI);
      //printf("MEANGAM=%f\n",*mean_gamma_ptr);      //printf("MEANGAM=%f\n",*mean_gamma_ptr);
      //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);      //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
       printf("sum,count_mask=%f,%d\n",sum,count_mask);
       printf("bh[200,200],bh_err[200,200]=%f,%f\n",bh[200,200],bh_err[200,200]);
   
      return 0;      return 0;
 } }
  
Line 245  int computeB_total(float *bx_err, float
Line 249  int computeB_total(float *bx_err, float
 /*===========================================*/ /*===========================================*/
 /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */ /* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */
  
 int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt, float *bt_err, float *mean_derivative_btotal_err_ptr, float *err_term1, float *err_term2)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 265  int computeBtotalderivative(float *bt, i
Line 269  int computeBtotalderivative(float *bt, i
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;                 derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
              err_term1[j * nx + i] = (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) ;
               }               }
           }           }
  
Line 274  int computeBtotalderivative(float *bt, i
Line 279  int computeBtotalderivative(float *bt, i
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;                 dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
              err_term2[j * nx + i] = (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) ;
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
Line 303  int computeBtotalderivative(float *bt, i
Line 310  int computeBtotalderivative(float *bt, i
              dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;              dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
           }           }
  
       // Calculate the sum only
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
Line 319  int computeBtotalderivative(float *bt, i
Line 326  int computeBtotalderivative(float *bt, i
                if isnan(derx_bt[j * nx + i]) continue;                if isnan(derx_bt[j * nx + i]) continue;
                if isnan(dery_bt[j * nx + i]) continue;                if isnan(dery_bt[j * nx + i]) continue;
                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */
                err += (((bt[(j+1) * nx + i]-bt[(j-1) * nx + i])*(bt[(j+1) * nx + i]-bt[(j-1) * nx + i])) * (bt_err[(j+1) * nx + i]*bt_err[(j+1) * nx + i] + bt_err[(j-1) * nx + i]*bt_err[(j-1) * nx + i])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+              err += err_term2[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ))+
                       (((bt[j * nx + (i+1)]-bt[j * nx + (i-1)])*(bt[j * nx + (i+1)]-bt[j * nx + (i-1)])) * (bt_err[j * nx + (i+1)]*bt_err[j * nx + (i+1)] + bt_err[j * nx + (i-1)]*bt_err[j * nx + (i-1)])) / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;                     err_term1[j * nx + i] / (16.0*( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
Line 337  int computeBtotalderivative(float *bt, i
Line 344  int computeBtotalderivative(float *bt, i
 /*===========================================*/ /*===========================================*/
 /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */ /* Example function 6:  Derivative of Bh SQRT( (dBh/dx)^2 + (dBh/dy)^2 ) */
  
 int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)  int computeBhderivative(float *bh, float *bh_err, int *dims, float *mean_derivative_bh_ptr, float *mean_derivative_bh_err_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh, float *err_term1, float *err_term2)
 { {
  
      int nx = dims[0];      int nx = dims[0];
Line 357  int computeBhderivative(float *bh, float
Line 364  int computeBhderivative(float *bh, float
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;                 derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
              err_term1[j * nx + i] = (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)]));
               }               }
           }           }
  
Line 366  int computeBhderivative(float *bh, float
Line 374  int computeBhderivative(float *bh, float
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;                 dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
             err_term2[j * nx + i] = (((bh[ (j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i]));
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
Line 411  int computeBhderivative(float *bh, float
Line 421  int computeBhderivative(float *bh, float
                if isnan(derx_bh[j * nx + i]) continue;                if isnan(derx_bh[j * nx + i]) continue;
                if isnan(dery_bh[j * nx + i]) continue;                if isnan(dery_bh[j * nx + i]) continue;
                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */
                err += (((bh[(j+1) * nx + i]-bh[(j-1) * nx + i])*(bh[(j+1) * nx + i]-bh[(j-1) * nx + i])) * (bh_err[(j+1) * nx + i]*bh_err[(j+1) * nx + i] + bh_err[(j-1) * nx + i]*bh_err[(j-1) * nx + i])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+              err += err_term2[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ))+
                       (((bh[j * nx + (i+1)]-bh[j * nx + (i-1)])*(bh[j * nx + (i+1)]-bh[j * nx + (i-1)])) * (bh_err[j * nx + (i+1)]*bh_err[j * nx + (i+1)] + bh_err[j * nx + (i-1)]*bh_err[j * nx + (i-1)])) / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;                     err_term1[j * nx + i] / (16.0*( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  )) ;
                count_mask++;                count_mask++;
             }             }
           }           }
Line 428  int computeBhderivative(float *bh, float
Line 438  int computeBhderivative(float *bh, float
 /*===========================================*/ /*===========================================*/
 /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */ /* Example function 7:  Derivative of B_vertical SQRT( (dBz/dx)^2 + (dBz/dy)^2 ) */
  
 int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)  int computeBzderivative(float *bz, float *bz_err, int *dims, float *mean_derivative_bz_ptr, float *mean_derivative_bz_err_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz, float *err_term1, float *err_term2)
 { {
  
         int nx = dims[0];         int nx = dims[0];
Line 447  int computeBzderivative(float *bz, float
Line 457  int computeBzderivative(float *bz, float
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                 if isnan(bz[j * nx + i]) continue;  
                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
              err_term1[j * nx + i] = (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)]));
               }               }
           }           }
  
Line 457  int computeBzderivative(float *bz, float
Line 467  int computeBzderivative(float *bz, float
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                 if isnan(bz[j * nx + i]) continue;  
                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
              err_term2[j * nx + i] = (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i]));
               }               }
           }           }
  
       /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
         /* consider the edges */      ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
           }           }
  
         i=nx-1;         i=nx-1;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
           }           }
  
         j=ny-1;         j=ny-1;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bz[j * nx + i]) continue;  
              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
           }           }
  
Line 508  int computeBzderivative(float *bz, float
Line 515  int computeBzderivative(float *bz, float
                if isnan(derx_bz[j * nx + i]) continue;                if isnan(derx_bz[j * nx + i]) continue;
                if isnan(dery_bz[j * nx + i]) continue;                if isnan(dery_bz[j * nx + i]) continue;
                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */
                err += (((bz[(j+1) * nx + i]-bz[(j-1) * nx + i])*(bz[(j+1) * nx + i]-bz[(j-1) * nx + i])) * (bz_err[(j+1) * nx + i]*bz_err[(j+1) * nx + i] + bz_err[(j-1) * nx + i]*bz_err[(j-1) * nx + i])) /              err += err_term2[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +
                         (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) +                     err_term1[j * nx + i] / (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;
                       (((bz[j * nx + (i+1)]-bz[j * nx + (i-1)])*(bz[j * nx + (i+1)]-bz[j * nx + (i-1)])) * (bz_err[j * nx + (i+1)]*bz_err[j * nx + (i+1)] + bz_err[j * nx + (i-1)]*bz_err[j * nx + (i-1)])) /  
                         (16.0*( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  )) ;  
                count_mask++;                count_mask++;
             }             }
           }           }
Line 562  int computeBzderivative(float *bz, float
Line 567  int computeBzderivative(float *bz, float
 //              float *noiseby, float *noisebz) //              float *noiseby, float *noisebz)
  
 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared, int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
               int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)
  
  
 { {
Line 581  int computeJz(float *bx_err, float *by_e
Line 586  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  if isnan(by[j * nx + i]) continue;  
                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
              err_term1[j * nx + i] = (by_err[j * nx + i+1])*(by_err[j * nx + i+1]) + (by_err[j * nx + i-1])*(by_err[j * nx + i-1]);
               }               }
           }           }
  
Line 590  int computeJz(float *bx_err, float *by_e
Line 595  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                  if isnan(bx[j * nx + i]) continue;  
                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
              err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
               }               }
           }           }
  
         // consider the edges      /* consider the edges for the arrays that contribute to the variable "sum" in the computation below.
       ignore the edges for the error terms as those arrays have been initialized to zero.
       this is okay because the error term will ultimately not include the edge pixels as they are selected out by the mask and bitmask arrays.*/
   
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
           }           }
  
         i=nx-1;         i=nx-1;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
              if isnan(by[j * nx + i]) continue;  
              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
           }           }
  
         j=ny-1;         j=ny-1;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
              if isnan(bx[j * nx + i]) continue;  
              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
           }           }
  
         for (i = 1; i <= nx-2; i++)  
       for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 1; j <= ny-2; j++)          for (j = 0; j <= ny-1; j++)
             {             {
                // calculate jz at all points                // calculate jz at all points
   
                jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix                jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix
                jz_err[j * nx + i]        = 0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +              jz_err[j * nx + i]        = 0.5*sqrt( err_term1[j * nx + i] + err_term2[j * nx + i] ) ;
                                                      (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;  
                jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);                jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
                count_mask++;                count_mask++;
   
             }             }
           }           }
         return 0;         return 0;
Line 831  int computeHelicity(float *jz_err, float
Line 833  int computeHelicity(float *jz_err, float
 /* Example function 12:  Sum of Absolute Value per polarity  */ /* Example function 12:  Sum of Absolute Value per polarity  */
  
 //  The Sum of the Absolute Value per polarity is defined as the following: //  The Sum of the Absolute Value per polarity is defined as the following:
 //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.  //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes per square arcsecond.
 //  The units of jz are in G/pix. In this case, we would have the following: //  The units of jz are in G/pix. In this case, we would have the following:
 //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
Line 868  int computeSumAbsPerPolarity(float *jz_e
Line 870  int computeSumAbsPerPolarity(float *jz_e
               }               }
           }           }
  
         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */      *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are Amperes per arcsecond */
         *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));         *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
         //printf("SAVNCPP=%g\n",*totaljzptr);         //printf("SAVNCPP=%g\n",*totaljzptr);
         //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);         //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
Line 990  int computeShearAngle(float *bx_err, flo
Line 992  int computeShearAngle(float *bx_err, flo
                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));                  magnitude_potential   = sqrt( (bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));
                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );                  magnitude_vector      = sqrt( (bx[j * nx + i]*bx[j * nx + i])   + (by[j * nx + i]*by[j * nx + i])   + (bz[j * nx + i]*bz[j * nx + i]) );
   
               //printf("i,j=%d,%d\n",i,j);
                  //printf("dotproduct=%f\n",dotproduct);                  //printf("dotproduct=%f\n",dotproduct);
                  //printf("magnitude_potential=%f\n",magnitude_potential);                  //printf("magnitude_potential=%f\n",magnitude_potential);
                  //printf("magnitude_vector=%f\n",magnitude_vector);                  //printf("magnitude_vector=%f\n",magnitude_vector);
               //printf("dotproduct/(magnitude_potential*magnitude_vector)=%f\n",dotproduct/(magnitude_potential*magnitude_vector));
               //printf("acos(dotproduct/(magnitude_potential*magnitude_vector))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector)));
               //printf("acos(dotproduct/(magnitude_potential*magnitude_vector)*(180./PI))=%f\n",acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI));
  
                  shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);                  shear_angle            = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
                  sumsum                  += shear_angle;                  sumsum                  += shear_angle;
Line 1034  int computeShearAngle(float *bx_err, flo
Line 1041  int computeShearAngle(float *bx_err, flo
         return 0;         return 0;
 } }
  
   /*===========================================*/
   /* Example function 15: R parameter as defined in Schrijver, 2007 */
   //
   // Note that there is a restriction on the function fsample()
   // If the following occurs:
   //      nx_out > floor((ny_in-1)/scale + 1)
   //      ny_out > floor((ny_in-1)/scale + 1),
   // where n*_out are the dimensions of the output array and n*_in
   // are the dimensions of the input array, fsample() will usually result
   // in a segfault (though not always, depending on how the segfault was accessed.)
   
   int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1,
                float *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
                float *pmap, int nx1, int ny1,
                int scale, float *p1pad, int nxp, int nyp, float *pmapn)
   
   {
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int index, index1;
       double sum = 0.0;
       double err = 0.0;
       *Rparam = 0.0;
       struct fresize_struct fresboxcar, fresgauss;
       struct fint_struct fints;
       float sigma = 10.0/2.3548;
   
       // set up convolution kernels
       init_fresize_boxcar(&fresboxcar,1,1);
       init_fresize_gaussian(&fresgauss,sigma,20,1);
   
       // =============== [STEP 1] ===============
       // bin the line-of-sight magnetogram down by a factor of scale
       fsample(los, rim, nx, ny, nx, nx1, ny1, nx1, scale, 0, 0, 0.0);
   
       // =============== [STEP 2] ===============
       // identify positive and negative pixels greater than +/- 150 gauss
       // and label those pixels with a 1.0 in arrays p1p0 and p1n0
       for (i = 0; i < nx1; i++)
       {
           for (j = 0; j < ny1; j++)
           {
               index = j * nx1 + i;
               if (rim[index] > 150)
                   p1p0[index]=1.0;
               else
                   p1p0[index]=0.0;
               if (rim[index] < -150)
                   p1n0[index]=1.0;
               else
                   p1n0[index]=0.0;
           }
       }
   
       // =============== [STEP 3] ===============
       // smooth each of the negative and positive pixel bitmaps
       fresize(&fresboxcar, p1p0, p1p, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
       fresize(&fresboxcar, p1n0, p1n, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 0.0);
   
       // =============== [STEP 4] ===============
       // find the pixels for which p1p and p1n are both equal to 1.
       // this defines the polarity inversion line
       for (i = 0; i < nx1; i++)
       {
           for (j = 0; j < ny1; j++)
           {
               index = j * nx1 + i;
               if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
                   p1[index]=1.0;
               else
                   p1[index]=0.0;
           }
       }
   
       // pad p1 with zeroes so that the gaussian colvolution in step 5
       // does not cut off data within hwidth of the edge
   
       // step i: zero p1pad
       for (i = 0; i < nxp; i++)
       {
           for (j = 0; j < nyp; j++)
           {
               index = j * nxp + i;
               p1pad[index]=0.0;
           }
       }
   
       // step ii: place p1 at the center of p1pad
       for (i = 0; i < nx1; i++)
       {
          for (j = 0; j < ny1; j++)
          {
               index  = j * nx1 + i;
               index1 = (j+20) * nxp + (i+20);
               p1pad[index1]=p1[index];
           }
       }
   
       // =============== [STEP 5] ===============
       // convolve the polarity inversion line map with a gaussian
       // to identify the region near the plarity inversion line
       // the resultant array is called pmap
       fresize(&fresgauss, p1pad, pmap, nxp, nyp, nxp, nxp, nyp, nxp, 0, 0, 0.0);
   
   
      // select out the nx1 x ny1 non-padded array  within pmap
       for (i = 0; i < nx1; i++)
       {
          for (j = 0; j < ny1; j++)
          {
               index  = j * nx1 + i;
               index1 = (j+20) * nxp + (i+20);
               pmapn[index]=pmap[index1];
           }
       }
   
       // =============== [STEP 6] ===============
       // the R parameter is calculated
       for (i = 0; i < nx1; i++)
       {
           for (j = 0; j < ny1; j++)
           {
               index = j * nx1 + i;
               if isnan(pmapn[index]) continue;
               if isnan(rim[index]) continue;
               sum += pmapn[index]*abs(rim[index]);
           }
       }
   
       if (sum < 1.0)
           *Rparam = 0.0;
       else
           *Rparam = log10(sum);
   
       //printf("R_VALUE=%f\n",*Rparam);
   
       free_fresize(&fresboxcar);
       free_fresize(&fresgauss);
   
       return 0;
   
   }
   
   /*===========================================*/
   /* Example function 16: Lorentz force as defined in Fisher, 2012 */
   //
   // This calculation is adapted from Xudong's code
   // at /proj/cgem/lorentz/apps/lorentz.c
   
   int computeLorentz(float *bx,  float *by, float *bz, float *fx, float *fy, float *fz, int *dims,
                      float *totfx_ptr, float *totfy_ptr, float *totfz_ptr, float *totbsq_ptr,
                      float *epsx_ptr, float *epsy_ptr, float *epsz_ptr, int *mask, int *bitmask,
                      float cdelt1, double rsun_ref, double rsun_obs)
   
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int nxny = nx*ny;
       int j = 0;
       int index;
       double totfx = 0, totfy = 0, totfz = 0;
       double bsq = 0, totbsq = 0;
       double epsx = 0, epsy = 0, epsz = 0;
       double area = cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       double k_h = -1.0 * area / (4. * PI) / 1.0e20;
       double k_z = area / (8. * PI) / 1.0e20;
   
       if (nx <= 0 || ny <= 0) return 1;
   
       for (int i = 0; i < nxny; i++)
       {
          if ( mask[i] < 70 || bitmask[i] < 30 ) continue;
          if isnan(bx[i]) continue;
          if isnan(by[i]) continue;
          if isnan(bz[i]) continue;
          fx[i]  = bx[i] * bz[i] * k_h;
          fy[i]  = by[i] * bz[i] * k_h;
          fz[i]  = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
          bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
          totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];
          totbsq += bsq;
       }
   
       *totfx_ptr  = totfx;
       *totfy_ptr  = totfy;
       *totfz_ptr  = totfz;
       *totbsq_ptr = totbsq;
       *epsx_ptr   = (totfx / k_h) / totbsq;
       *epsy_ptr   = (totfy / k_h) / totbsq;
       *epsz_ptr   = (totfz / k_z) / totbsq;
   
       //printf("TOTBSQ=%f\n",*totbsq_ptr);
   
       return 0;
   
   }
  
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/
  
Line 1116  void greenpot(float *bx, float *by, floa
Line 1322  void greenpot(float *bx, float *by, floa
       i2e = inx + iwindow;       i2e = inx + iwindow;
       if (i2s <   0){i2s =   0;}       if (i2s <   0){i2s =   0;}
       if (i2e > nnx){i2e = nnx;}       if (i2e > nnx){i2e = nnx;}
               //            for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
       for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)              for (j2=j2s;j2<2;j2++){for (i2=i2s;i2<2;i2++)
       {       {
         float val1 = bztmp[nnx*j2 + i2];         float val1 = bztmp[nnx*j2 + i2];
         float rr, r1, r2;  
 //        r1 = (float)(i2 - inx);  
 //        r2 = (float)(j2 - iny);  
 //        rr = r1*r1 + r2*r2;  
 //        if (rr < rwindow)  
 //        {  
           int   di, dj;           int   di, dj;
           di = abs(i2 - inx);           di = abs(i2 - inx);
           dj = abs(j2 - iny);           dj = abs(j2 - iny);
           sum = sum + val1 * rdist[nnx * dj + di] * dz;           sum = sum + val1 * rdist[nnx * dj + di] * dz;
 //        }                  printf("sum,val1,rdist[nnx * dj + di]=%f,%f,%f\n",sum,val1,rdist[nnx * dj + di]);
       } }       } }
       pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.       pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
     }     }
   } } // end of OpenMP parallelism   } } // end of OpenMP parallelism
  
       printf("bztmp[0]=%f\n",bztmp[0]);
       printf("bztmp[91746]=%f\n",bztmp[91746]);
       printf("pfpot[0]=%f\n",pfpot[0]);
       printf("pfpot[91746]=%f\n",pfpot[91746]);
   
 // #pragma omp parallel for private(inx) // #pragma omp parallel for private(inx)
   for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)   for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
   {   {
Line 1143  void greenpot(float *bx, float *by, floa
Line 1348  void greenpot(float *bx, float *by, floa
     by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;     by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
   } } // end of OpenMP parallelism   } } // end of OpenMP parallelism
  
       printf("bx[0]=%f\n",bx[0]);
       printf("by[0]=%f\n",by[0]);
       printf("bx[91746]=%f\n",bx[91746]);
       printf("by[91746]=%f\n",by[91746]);
   free(rdist);   free(rdist);
   free(pfpot);   free(pfpot);
   free(bztmp);   free(bztmp);
Line 1153  void greenpot(float *bx, float *by, floa
Line 1362  void greenpot(float *bx, float *by, floa
  
 char *sw_functions_version() // Returns CVS version of sw_functions.c char *sw_functions_version() // Returns CVS version of sw_functions.c
 { {
   return strdup("$Id$");      return strdup("$Id");
 } }
  
   
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/


Legend:
Removed from v.1.22  
changed lines
  Added in v.1.36

Karen Tian
Powered by
ViewCVS 0.9.4