(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.31 and 1.40

version 1.31, 2014/06/05 21:27:04 version 1.40, 2021/05/24 22:17:06
Line 1 
Line 1 
   
 /*=========================================== /*===========================================
  
  The following 14 functions calculate the following spaceweather indices:   The following functions calculate these spaceweather indices from the vector magnetic field data:
  
  USFLUX Total unsigned flux in Maxwells  USFLUX Total unsigned flux in Maxwells
  MEANGAM Mean inclination angle, gamma, in degrees  MEANGAM Mean inclination angle, gamma, in degrees
Line 17 
Line 18 
  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
    CMASK The total number of pixels that contributed to the calculation of all the indices listed above
   
    And these spaceweather indices from the line-of-sight magnetic field data:
    USFLUXL Total unsigned flux in Maxwells
    MEANGBL Mean value of the line-of-sight field gradient, in Gauss/Mm
    CMASKL The total number of pixels that contributed to the calculation of USFLUXL and MEANGBL
  R_VALUE Karel Schrijver's R parameter  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
Line 246  int computeB_total(float *bx_err, float
Line 253  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 273  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 283  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])) ;
         }         }
     }     }
  
       /* 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 304  int computeBtotalderivative(float *bt, i
Line 314  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 330  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 348  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 368  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 378  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]));
         }         }
     }     }
  
       /* 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 412  int computeBhderivative(float *bh, float
Line 425  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 442  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 461  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 471  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 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 509  int computeBzderivative(float *bz, float
Line 519  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 571  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 590  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 599  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 832  int computeHelicity(float *jz_err, float
Line 837  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 869  int computeSumAbsPerPolarity(float *jz_e
Line 874  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 1023  int computeShearAngle(float *bx_err, flo
Line 1028  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 = (sumsum)/(count);                 /* Units are degrees */     *meanshear_angleptr = (sumsum)/(count);                 /* Units are degrees */
   
       // For the error in the mean 3D shear angle:
       // If count_mask is 0, then we run into a divide by zero error. In this case, set *meanshear_angle_err_ptr to NAN
       // If count_mask is greater than zero, then compute the error.
       if (count_mask == 0)
           *meanshear_angle_err_ptr = NAN;
       else
     *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);     *meanshear_angle_err_ptr = (sqrt(err)/count_mask)*(180./PI);
  
     /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */     /* The area here is a fractional area -- the % of the total area. This has no error associated with it. */
     *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);     *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);
  
     //printf("MEANSHR=%f\n",*meanshear_angleptr);     //printf("MEANSHR=%f\n",*meanshear_angleptr);
     //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);      //printf("ERRMSHA=%f\n",*meanshear_angle_err_ptr);
     //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);     //printf("SHRGT45=%f\n",*area_w_shear_gt_45ptr);
   
         return 0;         return 0;
 } }
  
Line 1160  int computeR(float *bz_err, float *los,
Line 1171  int computeR(float *bz_err, float *los,
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
         {         {
             index = j * nx1 + i;             index = j * nx1 + i;
               if isnan(pmapn[index]) continue;
               if isnan(rim[index]) continue;
             sum += pmapn[index]*abs(rim[index]);             sum += pmapn[index]*abs(rim[index]);
         }         }
     }     }
Line 1169  int computeR(float *bz_err, float *los,
Line 1182  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);
  
Line 1201  int computeLorentz(float *bx, float *by
Line 1216  int computeLorentz(float *bx, float *by
     double k_h = -1.0 * area / (4. * PI) / 1.0e20;     double k_h = -1.0 * area / (4. * PI) / 1.0e20;
     double k_z = area / (8. * PI) / 1.0e20;     double k_z = area / (8. * PI) / 1.0e20;
  
     /* Multiplier */  
     float vectorMulti[] = {1.,-1.,1.};  
   
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
     for (int i = 0; i < nxny; i++)     for (int i = 0; i < nxny; i++)
Line 1212  int computeLorentz(float *bx, float *by
Line 1224  int computeLorentz(float *bx, float *by
        if isnan(bx[i]) continue;        if isnan(bx[i]) continue;
        if isnan(by[i]) continue;        if isnan(by[i]) continue;
        if isnan(bz[i]) continue;        if isnan(bz[i]) continue;
        fx[i]  = (bx[i] * vectorMulti[0]) * (bz[i] * vectorMulti[2]) * k_h;         fx[i]  = bx[i] * bz[i] * k_h;
        fy[i]  = (by[i] * vectorMulti[1]) * (bz[i] * vectorMulti[2]) * 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;        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];        bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
        totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];        totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];
Line 1228  int computeLorentz(float *bx, float *by
Line 1240  int computeLorentz(float *bx, float *by
     *epsy_ptr   = (totfy / k_h) / totbsq;     *epsy_ptr   = (totfy / k_h) / totbsq;
     *epsz_ptr   = (totfz / k_z) / totbsq;     *epsz_ptr   = (totfz / k_z) / totbsq;
  
       //printf("TOTBSQ=%f\n",*totbsq_ptr);
   
       return 0;
   
   }
   
   /*===========================================*/
   
   /* Example function 17: Compute total unsigned flux in units of G/cm^2 on the LOS field */
   
   //  To compute the unsigned flux, we simply calculate
   //  flux = surface integral [(vector LOS) dot (normal vector)],
   //       = surface integral [(magnitude LOS)*(magnitude normal)*(cos theta)].
   //  However, since the field is radial, we will assume cos theta = 1.
   //  Therefore the pixels only need to be corrected for the projection.
   
   //  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
   //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
   //  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
   //  =Gauss*cm^2
   
   int computeAbsFlux_los(float *los, int *dims, float *absFlux_los,
                          float *mean_vf_los_ptr, float *count_mask_los_ptr,
                          int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
   
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask_los = 0;
       double sum = 0.0;
       *absFlux_los = 0.0;
       *mean_vf_los_ptr = 0.0;
   
   
       if (nx <= 0 || ny <= 0) return 1;
   
           for (i = 0; i < nx; i++)
           {
              for (j = 0; j < ny; j++)
              {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if isnan(los[j * nx + i]) continue;
               sum += (fabs(los[j * nx + i]));
               count_mask_los++;
              }
           }
   
       *mean_vf_los_ptr     = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       *count_mask_los_ptr  = count_mask_los;
   
       printf("USFLUXL=%f\n",*mean_vf_los_ptr);
       printf("CMASKL=%f\n",*count_mask_los_ptr);
   
     return 0;     return 0;
   }
   
   /*===========================================*/
   /* Example function 18:  Derivative of B_LOS (approximately B_vertical) = SQRT( ( dLOS/dx )^2 + ( dLOS/dy )^2 ) */
   
   int computeLOSderivative(float *los, int *dims, float *mean_derivative_los_ptr, int *bitmask, float *derx_los, float *dery_los)
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       *mean_derivative_los_ptr = 0.0;
   
       if (nx <= 0 || ny <= 0) return 1;
   
       /* brute force method of calculating the derivative (no consideration for edges) */
       for (i = 1; i <= nx-2; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
              derx_los[j * nx + i] = (los[j * nx + i+1] - los[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 (j = 1; j <= ny-2; j++)
           {
              dery_los[j * nx + i] = (los[(j+1) * nx + i] - los[(j-1) * nx + i])*0.5;
           }
       }
  
       /* 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;
       for (j = 0; j <= ny-1; j++)
       {
           derx_los[j * nx + i] = ( (-3*los[j * nx + i]) + (4*los[j * nx + (i+1)]) - (los[j * nx + (i+2)]) )*0.5;
       }
   
       i=nx-1;
       for (j = 0; j <= ny-1; j++)
       {
           derx_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[j * nx + (i-1)]) - (-los[j * nx + (i-2)]) )*0.5;
       }
   
       j=0;
       for (i = 0; i <= nx-1; i++)
       {
           dery_los[j * nx + i] = ( (-3*los[j*nx + i]) + (4*los[(j+1) * nx + i]) - (los[(j+2) * nx + i]) )*0.5;
       }
   
       j=ny-1;
       for (i = 0; i <= nx-1; i++)
       {
           dery_los[j * nx + i] = ( (3*los[j * nx + i]) + (-4*los[(j-1) * nx + i]) - (-los[(j-2) * nx + i]) )*0.5;
       }
   
   
       for (i = 0; i <= nx-1; i++)
       {
           for (j = 0; j <= ny-1; j++)
           {
               if ( bitmask[j * nx + i] < 30 ) continue;
               if ( (derx_los[j * nx + i] + dery_los[j * nx + i]) == 0) continue;
               if isnan(los[j * nx + i])      continue;
               if isnan(los[(j+1) * nx + i])  continue;
               if isnan(los[(j-1) * nx + i])  continue;
               if isnan(los[j * nx + i-1])    continue;
               if isnan(los[j * nx + i+1])    continue;
               if isnan(derx_los[j * nx + i]) continue;
               if isnan(dery_los[j * nx + i]) continue;
               sum += sqrt( derx_los[j * nx + i]*derx_los[j * nx + i]  + dery_los[j * nx + i]*dery_los[j * nx + i] ); /* Units of Gauss */
               count_mask++;
           }
       }
   
       *mean_derivative_los_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
   
       printf("MEANGBL=%f\n",*mean_derivative_los_ptr);
   
           return 0;
 } }
  
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/

Legend:
Removed from v.1.31  
changed lines
  Added in v.1.40
Karen Tian
 Powered by
ViewCVS 0.9.4