(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.2 and 1.4

version 1.2, 2012/08/27 19:55:49 version 1.4, 2012/12/18 23:05:10
Line 18 
Line 18 
     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
  
    The indices are calculated on the pixels in which the disambig bitmap equals 5 or 7:     The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
     5: pixels for which the radial acute disambiguation solution was chosen     pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD
     7: pixels for which the radial acute and NRWA disambiguation agree     coordinate bitmaps are interpolated.
   
      In the CCD coordinates, this means that we are selecting the pixels that equal 90 in conf_disambig
      and the pixels that equal 33 or 44 in bitmap. Here are the definitions of the pixel values:
   
      For conf_disambig:
      50 : not all solutions agree (weak field method applied)
      60 : not all solutions agree (weak field + annealed)
      90 : all solutions agree (strong field + annealed)
       0 : not disambiguated
   
      For bitmap:
      1  : weak field outside smooth bounding curve
      2  : strong field outside smooth bounding curve
      33 : weak field inside smooth bounding curve
      34 : strong field inside smooth bounding curve
  
    Written by Monica Bobra 15 August 2012    Written by Monica Bobra 15 August 2012
    Potential Field code (appended) written by Keiji Hayashi    Potential Field code (appended) written by Keiji Hayashi
Line 52 
Line 67 
 //  7: pixels for which the radial acute and NRWA disambiguation agree //  7: pixels for which the radial acute and NRWA disambiguation agree
  
 int computeAbsFlux(float *bz, int *dims, float *absFlux, int computeAbsFlux(float *bz, int *dims, float *absFlux,
                                    float *mean_vf_ptr, int *mask,                                     float *mean_vf_ptr, int *mask,  int *bitmask,
                                    float cdelt1, double rsun_ref, double rsun_obs)                                    float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 70  int computeAbsFlux(float *bz, int *dims,
Line 85  int computeAbsFlux(float *bz, int *dims,
         {         {
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   if ( mask[j * nx + i] < 70 ) continue;                    if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                     if isnan(bz[j * nx + i]) continue;
                     //printf("bz[j * nx + i]=%f\n",bz[j * nx + i]);
                   sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
                   //sum += (fabs(bz[j * nx + i]))*inverseMu[j * nx + i]; // use this with mu function  
                   count_mask++;                   count_mask++;
                 }                 }
         }         }
Line 88  int computeAbsFlux(float *bz, int *dims,
Line 104  int computeAbsFlux(float *bz, int *dims,
 // Native units of Bh are Gauss // Native units of Bh are Gauss
  
 int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, int computeBh(float *bx, float *by, float *bz, float *bh, int *dims,
                           float *mean_hf_ptr, int *mask)                            float *mean_hf_ptr, int *mask, int *bitmask)
  
 { {
  
Line 103  int computeBh(float *bx, float *by, floa
Line 119  int computeBh(float *bx, float *by, floa
           {           {
             for (i = 0; i < nx; i++)             for (i = 0; i < nx; i++)
               {               {
                   if isnan(bx[j * nx + i]) continue;
                   if isnan(by[j * nx + i]) continue;
                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
                 sum += bh[j * nx + i];                 sum += bh[j * nx + i];
                 count_mask++;                 count_mask++;
Line 120  int computeBh(float *bx, float *by, floa
Line 138  int computeBh(float *bx, float *by, floa
  
  
 int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims, int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims,
                                  float *mean_gamma_ptr, int *mask)                                   float *mean_gamma_ptr, int *mask, int *bitmask)
 { {
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
     int i, j, count_mask=0;     int i, j, count_mask=0;
Line 137  int computeGamma(float *bx, float *by, f
Line 155  int computeGamma(float *bx, float *by, f
               {               {
                 if (bh[j * nx + i] > 100)                 if (bh[j * nx + i] > 100)
                   {                   {
                     if (mask[j * nx + i] < 70 ) continue;                      if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                       if isnan(bz[j * nx + i]) continue;
                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
                     count_mask++;                     count_mask++;
                   }                   }
Line 153  int computeGamma(float *bx, float *by, f
Line 172  int computeGamma(float *bx, float *by, f
 /* Example function 4: Calculate B_Total*/ /* Example function 4: Calculate B_Total*/
 // Native units of B_Total are in gauss // Native units of B_Total are in gauss
  
 int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask)  int computeB_total(float *bx, float *by, float *bz, float *bt, int *dims, int *mask, int *bitmask)
 { {
  
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
Line 165  int computeB_total(float *bx, float *by,
Line 184  int computeB_total(float *bx, float *by,
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                   if isnan(bx[j * nx + i]) continue;
                   if isnan(by[j * nx + i]) continue;
                   if isnan(bz[j * nx + i]) continue;
                 bt[j * nx + i] = 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]);                 bt[j * nx + i] = 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]);
               }               }
           }           }
Line 174  int computeB_total(float *bx, float *by,
Line 196  int computeB_total(float *bx, float *by,
 /*===========================================*/ /*===========================================*/
 /* 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, float *derx_bt, float *dery_bt)  int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt)
 { {
  
     int nx = dims[0], ny = dims[1];     int nx = dims[0], ny = dims[1];
Line 250  int computeBtotalderivative(float *bt, i
Line 272  int computeBtotalderivative(float *bt, i
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;                // if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */
                count_mask++;                count_mask++;
             }             }
Line 265  int computeBtotalderivative(float *bt, i
Line 287  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, int *dims, float *mean_derivative_bh_ptr, int *mask, float *derx_bh, float *dery_bh)  int computeBhderivative(float *bh, int *dims, float *mean_derivative_bh_ptr, int *mask, int *bitmask, float *derx_bh, float *dery_bh)
 { {
  
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 340  int computeBhderivative(float *bh, int *
Line 362  int computeBhderivative(float *bh, int *
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;                // if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */
                count_mask++;                count_mask++;
             }             }
Line 353  int computeBhderivative(float *bh, int *
Line 375  int computeBhderivative(float *bh, int *
 /*===========================================*/ /*===========================================*/
 /* 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, int *dims, float *mean_derivative_bz_ptr, int *mask, float *derx_bz, float *dery_bz)  int computeBzderivative(float *bz, int *dims, float *mean_derivative_bz_ptr, int *mask, int *bitmask, float *derx_bz, float *dery_bz)
 { {
  
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 369  int computeBzderivative(float *bz, int *
Line 391  int computeBzderivative(float *bz, int *
           {           {
             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;
               }               }
           }           }
Line 378  int computeBzderivative(float *bz, int *
Line 401  int computeBzderivative(float *bz, int *
           {           {
             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;
               }               }
           }           }
Line 387  int computeBzderivative(float *bz, int *
Line 411  int computeBzderivative(float *bz, int *
         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 428  int computeBzderivative(float *bz, int *
Line 456  int computeBzderivative(float *bz, int *
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if isnan(bz[j * nx + i]) continue;
                  if isnan(derx_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 */
                count_mask++;                count_mask++;
             }             }
Line 459  int computeBzderivative(float *bz, int *
Line 490  int computeBzderivative(float *bz, int *
 //  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or //  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
 //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
 //  where a Tesla is represented as a Newton/Ampere*meter. //  where a Tesla is represented as a Newton/Ampere*meter.
   //
 //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
 //  In that case, we would have the following: //  In that case, we would have the following:
 //  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or //  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
Line 471  int computeBzderivative(float *bz, int *
Line 503  int computeBzderivative(float *bz, int *
 //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.) //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
  
 int computeJz(float *bx, float *by, int *dims, float *jz, int computeJz(float *bx, float *by, int *dims, float *jz,
                           float *mean_jz_ptr, float *us_i_ptr, int *mask,                            float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
  
  
Line 492  int computeJz(float *bx, float *by, int
Line 524  int computeJz(float *bx, float *by, int
           {           {
             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;
               }               }
           }           }
Line 501  int computeJz(float *bx, float *by, int
Line 534  int computeJz(float *bx, float *by, int
           {           {
             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;
               }               }
           }           }
Line 510  int computeJz(float *bx, float *by, int
Line 544  int computeJz(float *bx, float *by, int
         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;
           }           }
  
         /* Just some print statements  
         for (i = 0; i < nx; i++)  
           {  
              for (j = 0; j < ny; j++)  
               {  
               printf("j=%d\n",j);  
               printf("i=%d\n",i);  
               printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);  
               printf("derx[j*nx+i]=%f\n",derx[j*nx+i]);  
               printf("bx[j*nx+i]=%f\n",bx[j*nx+i]);  
               printf("by[j*nx+i]=%f\n",by[j*nx+i]);  
               }  
           }  
         */  
   
  
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                if (mask[j * nx + i] < 70 ) continue;                 if isnan(derx[j * nx + i]) continue;
                  if isnan(dery[j * nx + i]) continue;
                curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
                us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */                us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */
                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 */
   
                count_mask++;                count_mask++;
             }             }
           }           }
Line 566  int computeJz(float *bx, float *by, int
Line 589  int computeJz(float *bx, float *by, int
         printf("cdelt1, what is it?=%f\n",cdelt1);         printf("cdelt1, what is it?=%f\n",cdelt1);
         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */
         printf("count_mask=%d\n",count_mask);         printf("count_mask=%d\n",count_mask);
         *us_i_ptr        = (us_i);                   /* us_i gets populated as MEANJZD */          *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
         return 0;         return 0;
  
 } }
Line 582  int computeJz(float *bx, float *by, int
Line 605  int computeJz(float *bx, float *by, int
 //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
 //                               = 1/Mm //                               = 1/Mm
  
 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask,  int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask, int *bitmask,
                  float cdelt1, double rsun_ref, double rsun_obs)                  float cdelt1, double rsun_ref, double rsun_obs)
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 597  int computeAlpha(float *bz, int *dims, f
Line 620  int computeAlpha(float *bz, int *dims, f
           {           {
             for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
               {               {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                 if isnan(jz[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
Line 624  int computeAlpha(float *bz, int *dims, f
Line 647  int computeAlpha(float *bz, int *dims, f
  
  
 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr, int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,
                                         float *total_abs_ih_ptr, int *mask, float cdelt1, double rsun_ref, double rsun_obs)                                          float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
  
Line 640  int computeHelicity(float *bz, int *dims
Line 663  int computeHelicity(float *bz, int *dims
         {         {
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                 if isnan(jz[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
Line 670  int computeHelicity(float *bz, int *dims
Line 693  int computeHelicity(float *bz, int *dims
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
  
 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr, int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,
                                                          int *mask, float cdelt1, double rsun_ref, double rsun_obs)                                                           int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
Line 685  int computeSumAbsPerPolarity(float *bz,
Line 708  int computeSumAbsPerPolarity(float *bz,
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
Line 708  int computeSumAbsPerPolarity(float *bz,
Line 732  int computeSumAbsPerPolarity(float *bz,
 // = erg/cm(1/pix^2) // = erg/cm(1/pix^2)
  
 int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims,
                                           float *meanpotptr, float *totpotptr, int *mask,                                            float *meanpotptr, float *totpotptr, int *mask, int *bitmask,
                                           float cdelt1, double rsun_ref, double rsun_obs)                                           float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 725  int computeFreeEnergy(float *bx, float *
Line 749  int computeFreeEnergy(float *bx, float *
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  if (mask[j * nx + i] < 70 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                    if isnan(bx[j * nx + i]) continue;
                    if isnan(by[j * nx + i]) continue;
                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);
                  count_mask++;                  count_mask++;
               }               }
Line 742  int computeFreeEnergy(float *bx, float *
Line 768  int computeFreeEnergy(float *bx, float *
 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,
                                           float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,                                           float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,
                                           int *mask)                                            int *mask, int *bitmask)
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
         int i, j;         int i, j;
Line 758  int computeShearAngle(float *bx, float *
Line 784  int computeShearAngle(float *bx, float *
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  if (mask[j * nx + i] < 70 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if isnan(bpx[j * nx + i]) continue;                  if isnan(bpx[j * nx + i]) continue;
                  if isnan(bpy[j * nx + i]) continue;                  if isnan(bpy[j * nx + i]) continue;
                  if isnan(bpz[j * nx + i]) continue;                  if isnan(bpz[j * nx + i]) continue;
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                    if isnan(bx[j * nx + i]) continue;
                    if isnan(by[j * nx + i]) continue;
                  /* For mean 3D shear angle, area with shear greater than 45*/                  /* For mean 3D shear angle, area with shear greater than 45*/
                  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]));


Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4

Karen Tian
Powered by
ViewCVS 0.9.4