(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.29 and 1.34

version 1.29, 2014/05/16 21:55:54 version 1.34, 2015/02/27 19:49:43
Line 246  int computeB_total(float *bx_err, float
Line 246  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_termAt, float *err_termBt)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 266  int computeBtotalderivative(float *bt, i
Line 266  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_termAt[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 275  int computeBtotalderivative(float *bt, i
Line 276  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_termBt[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])) ;
         }         }
     }     }
  
Line 304  int computeBtotalderivative(float *bt, i
Line 306  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 320  int computeBtotalderivative(float *bt, i
Line 322  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_termBt[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_termAt[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 338  int computeBtotalderivative(float *bt, i
Line 340  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_termAh, float *err_termBh)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 358  int computeBhderivative(float *bh, float
Line 360  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_termAh[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 367  int computeBhderivative(float *bh, float
Line 370  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_termBh[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]));
         }         }
     }     }
  
Line 412  int computeBhderivative(float *bh, float
Line 416  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_termBh[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_termAh[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 429  int computeBhderivative(float *bh, float
Line 433  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_termA, float *err_termB)
 { {
  
     int nx = dims[0];     int nx = dims[0];
Line 448  int computeBzderivative(float *bz, float
Line 452  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_termA[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 458  int computeBzderivative(float *bz, float
Line 462  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_termB[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 err_termA and err_term B 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 509  int computeBzderivative(float *bz, float
Line 510  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_termB[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_termA[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 563  int computeBzderivative(float *bz, float
Line 562  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 582  int computeJz(float *bx_err, float *by_e
Line 581  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 591  int computeJz(float *bx_err, float *by_e
Line 590  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]);
         }         }
     }     }
  
Line 600  int computeJz(float *bx_err, float *by_e
Line 599  int computeJz(float *bx_err, float *by_e
     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;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     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;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     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;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
     }     }
  
     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;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
   
     }     }
  
     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 832  int computeHelicity(float *jz_err, float
Line 830  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 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 869  int computeSumAbsPerPolarity(float *jz_e
Line 867  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 1048  int computeShearAngle(float *bx_err, flo
Line 1046  int computeShearAngle(float *bx_err, flo
  
 int computeR(float *bz_err, float *los, int *dims, float *Rparam, float cdelt1, 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 *rim, float *p1p0, float *p1n0, float *p1p, float *p1n, float *p1,
              float *pmap, int nx1, int ny1)               float *pmap, int nx1, int ny1,
                int scale, float *p1pad, int nxp, int nyp, float *pmapn)
   
 { {
     int nx = dims[0];     int nx = dims[0];
     int ny = dims[1];     int ny = dims[1];
     int i = 0;     int i = 0;
     int j = 0;     int j = 0;
     int index;      int index, index1;
     double sum = 0.0;     double sum = 0.0;
     double err = 0.0;     double err = 0.0;
     *Rparam = 0.0;     *Rparam = 0.0;
     struct fresize_struct fresboxcar, fresgauss;     struct fresize_struct fresboxcar, fresgauss;
     int scale = round(2.0/cdelt1);      struct fint_struct fints;
     float sigma = 10.0/2.3548;     float sigma = 10.0/2.3548;
  
       // set up convolution kernels
     init_fresize_boxcar(&fresboxcar,1,1);     init_fresize_boxcar(&fresboxcar,1,1);
   
     // set up convolution kernel  
     init_fresize_gaussian(&fresgauss,sigma,20,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);     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 (i = 0; i < nx1; i++)
     {     {
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
Line 1085  int computeR(float *bz_err, float *los,
Line 1089  int computeR(float *bz_err, float *los,
         }         }
     }     }
  
       // =============== [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, 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);     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 (i = 0; i < nx1; i++)
     {     {
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
         {         {
             index = j * nx1 + i;             index = j * nx1 + i;
             if (p1p[index] > 0 && p1n[index] > 0)              if ((p1p[index] > 0.0) && (p1n[index] > 0.0))
                 p1[index]=1.0;                 p1[index]=1.0;
             else             else
                 p1[index]=0.0;                 p1[index]=0.0;
         }         }
     }     }
  
     fresize(&fresgauss, p1, pmap, nx1, ny1, nx1, nx1, ny1, nx1, 0, 0, 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 (i = 0; i < nx1; i++)
     {     {
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
         {         {
             index = j * nx1 + i;             index = j * nx1 + i;
             sum += pmap[index]*abs(rim[index]);              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]);
         }         }
     }     }
  
Line 1116  int computeR(float *bz_err, float *los,
Line 1169  int computeR(float *bz_err, float *los,
     else     else
         *Rparam = log10(sum);         *Rparam = log10(sum);
  
       //printf("R_VALUE=%f\n",*Rparam);
   
     free_fresize(&fresboxcar);     free_fresize(&fresboxcar);
     free_fresize(&fresgauss);     free_fresize(&fresgauss);
  
     return 0;     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 =========================*/
  


Legend:
Removed from v.1.29  
changed lines
  Added in v.1.34

Karen Tian
Powered by
ViewCVS 0.9.4