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

version 1.2, 2012/08/27 19:55:49 version 1.8, 2013/02/09 02:39:20
Line 1 
Line 1 
 /*=========================================== /*===========================================
  
    The following 13 functions calculate the following spaceweather indices:     The following 14 functions calculate the following spaceweather indices:
  
     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 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 66  int computeAbsFlux(float *bz, int *dims,
Line 81  int computeAbsFlux(float *bz, int *dims,
     *absFlux = 0.0;     *absFlux = 0.0;
     *mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;
  
         for (j = 0; j < ny; j++)  
         {  
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   if ( mask[j * nx + i] < 70 ) continue;                  for (j = 0; j < ny; j++)
                   {
                     if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                     if isnan(bz[j * nx + i]) continue;
                   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 103  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 99  int computeBh(float *bx, float *by, floa
Line 114  int computeBh(float *bx, float *by, floa
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
         for (j = 0; j < ny; j++)  
           {  
             for (i = 0; i < nx; i++)             for (i = 0; i < nx; i++)
               {               {
               for (j = 0; j < ny; j++)
                 {
                   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 137  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 154  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 171  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 183  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 195  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 231  int computeBtotalderivative(float *bt, i
Line 252  int computeBtotalderivative(float *bt, i
           }           }
  
  
         /* 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_bt[j*nx+i]=%f\n",dery_bt[j*nx+i]);  
               printf("derx_bt[j*nx+i]=%f\n",derx_bt[j*nx+i]);  
               printf("bt[j*nx+i]=%f\n",bt[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_bt[j * nx + i]-dery_bt[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_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 */
                count_mask++;                count_mask++;
             }             }
Line 265  int computeBtotalderivative(float *bt, i
Line 273  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 321  int computeBhderivative(float *bh, int *
Line 329  int computeBhderivative(float *bh, int *
           }           }
  
  
         /*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_bh[j*nx+i]=%f\n",dery_bh[j*nx+i]);  
               printf("derx_bh[j*nx+i]=%f\n",derx_bh[j*nx+i]);  
               printf("bh[j*nx+i]=%f\n",bh[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_bh[j * nx + i]-dery_bh[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_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 */
                count_mask++;                count_mask++;
             }             }
Line 353  int computeBhderivative(float *bh, int *
Line 348  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 364  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 374  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 384  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;
           }           }
  
  
         /*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_bz[j*nx+i]=%f\n",dery_bz[j*nx+i]);  
               printf("derx_bz[j*nx+i]=%f\n",derx_bz[j*nx+i]);  
               printf("bz[j*nx+i]=%f\n",bz[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_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 439  int computeBzderivative(float *bz, int *
Line 429  int computeBzderivative(float *bz, int *
 } }
  
 /*===========================================*/ /*===========================================*/
   
 /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */ /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */
  
 //  In discretized space like data pixels, //  In discretized space like data pixels,
Line 456  int computeBzderivative(float *bz, int *
Line 445  int computeBzderivative(float *bz, int *
 // //
 //  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003), //  To change units from Gauss/pixel to mA/m^2 (the units for Jz in Leka and Barnes, 2003),
 //  one must perform the following unit conversions: //  one must perform the following unit conversions:
 //  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or  //  (Gauss)(1/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
 //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),  //  (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), or
   //  (Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(1000.),
 //  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
 //  jz * (35.0) //  jz * (35.0)
 // //
 //  The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following: //  The units of total unsigned vertical current (us_i) are simply in A. In this case, we would have the following:
 //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(1000.)  //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(0.00010)(1/MUNAUGHT)(CDELT1)(CDELT1)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)
 //  =(Gauss/pix)(1/CDELT1)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)  //  = (Gauss/pix)(0.00010)(1/MUNAUGHT)(CDELT1)(RSUN_REF/RSUN_OBS)
 //  =(Gauss/pix)(1/0.5)(10^-4)(4*PI*10^7)(722500)(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,                           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 482  int computeJz(float *bx, float *by, int
Line 471  int computeJz(float *bx, float *by, int
         int i, j, count_mask=0;         int i, j, count_mask=0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
   
         *mean_jz_ptr = 0.0;  
         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
  
  
Line 492  int computeJz(float *bx, float *by, int
Line 479  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 489  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 499  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 (i = 0; i <= nx-1; i++)
           {           {
              for (j = 0; j < ny; j++)              for (j = 0; j <= ny-1; j++)
               {               {
               printf("j=%d\n",j);                 /* calculate jz at all points */
               printf("i=%d\n",i);                 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                              /* jz is in units of Gauss/pix */
               printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);                 count_mask++;
               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]);  
               }               }
           }           }
         */  
           return 0;
   }
   
   /*===========================================*/
   
   /* Example function 9:  Compute quantities on smoothed Jz array */
   
   // All of the subsequent functions, including this one, use a smoothed Jz array. The smoothing is performed by Jesper's
   // fresize routines. These routines are located at /cvs/JSOC/proj/libs/interpolate. A Gaussian with a FWHM of 4 pixels
   // and truncation width of 12 pixels is used to smooth the array; however, a quick analysis shows that the mean values
   // of qualities like Jz and helicity do not change much as a result of smoothing. The smoothed array will, of course,
   // give a lower total Jz as the stron field pixels have been smoothed out to neighboring weaker field pixels.
   
   int computeJzsmooth(float *bx, float *by, int *dims, float *jz_smooth,
                             float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
                             float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
   
   {
   
           int nx = dims[0], ny = dims[1];
           int i, j, count_mask=0;
   
           if (nx <= 0 || ny <= 0) return 1;
   
           *mean_jz_ptr = 0.0;
           float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
  
  
           /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions*/
         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;
                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 */                 if isnan(dery[j * nx + i]) continue;
                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 */                 if isnan(jz_smooth[j * nx + i]) continue;
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                 //printf("%d,%d,%f\n",i,j,jz_smooth[j * nx + i]);
                  curl +=     (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
                  us_i += fabs(jz_smooth[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
                count_mask++;                count_mask++;
             }             }
           }           }
  
           /* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
         mean_curl        = (curl/count_mask);         mean_curl        = (curl/count_mask);
         printf("mean_curl=%f\n",mean_curl);         printf("mean_curl=%f\n",mean_curl);
         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;
  
 } }
  
   
 /*===========================================*/ /*===========================================*/
 /* Example function 9:  Twist Parameter, alpha */  
  
 // The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm  /* Example function 10:  Twist Parameter, alpha */
   
   // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
   // for alpha is calculated in the following way (different from Leka and Barnes' approach):
   
          // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
          // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
          // avg alpha = avg Jz / avg Bz
   
   // The sign is assigned as follows:
   // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
   // If this value is > 0, then alpha is > 0.
   // If this value is < 0, then alpha is <0.
   //
   // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
   // If this value is > 0, then alpha is < 0.
   // If this value is < 0, then alpha is > 0.
   
   // The units of alpha are in 1/Mm
 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
 // //
 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
 //                               = (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_smooth, 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];
         int i, j, count_mask=0;          int i, j=0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_alpha_ptr = 0.0;         *mean_alpha_ptr = 0.0;
         float aa, bb, cc, bznew, alpha2, sum=0.0;          float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=0.0;
  
         for (i = 1; i < nx-1; i++)         for (i = 1; i < nx-1; i++)
           {           {
             for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
               {               {
                 if (mask[j * nx + i] < 70 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                  if isnan(jz_smooth[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;
                 sum += (jz[j * nx + i] / bz[j * nx + i])*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */                  if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i]);
                 count_mask++;                  if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i]);
                   if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);
                   if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
                   sum5 += bz[j * nx + i];
               }               }
           }           }
  
         printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);          sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
         printf("count_mask=%d\n",count_mask);  
         printf("sum=%f\n",sum);          /* Determine the sign of alpha */
         *mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */          if ((sum5 > 0) && (sum3 >  0)) sum=sum;
           if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
           if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
           if ((sum5 < 0) && (sum4 >  0)) sum=-sum;
   
           *mean_alpha_ptr = sum; /* Units are 1/Mm */
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 10:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */  /* Example function 11:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
  
 //  The current helicity is defined as Bz*Jz and the units are G^2 / m //  The current helicity is defined as Bz*Jz and the units are G^2 / m
 //  The units of Jz are in G/pix; the units of Bz are in G. //  The units of Jz are in G/pix; the units of Bz are in G.
Line 623  int computeAlpha(float *bz, int *dims, f
Line 667  int computeAlpha(float *bz, int *dims, f
 //                                                      = G^2 / m. //                                                      = G^2 / m.
  
  
 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_smooth, 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 636  int computeHelicity(float *bz, int *dims
Line 680  int computeHelicity(float *bz, int *dims
         *mean_ih_ptr = 0.0;         *mean_ih_ptr = 0.0;
         float sum=0.0, sum2=0.0;         float sum=0.0, sum2=0.0;
  
         for (j = 0; j < ny; j++)  
         {  
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                 if (mask[j * nx + i] < 70 ) continue;                  for (j = 0; j < ny; j++)
                 if isnan(jz[j * nx + i]) continue;                  {
                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(jz_smooth[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;
                 if (jz[j * nx + i] == 0.0) continue;                  if (jz_smooth[j * nx + i] == 0.0) continue;
                 sum  +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum  +=     (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
                 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
                 count_mask++;                 count_mask++;
                 }                 }
          }          }
Line 661  int computeHelicity(float *bz, int *dims
Line 705  int computeHelicity(float *bz, int *dims
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 11:  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.
Line 669  int computeHelicity(float *bz, int *dims
Line 713  int computeHelicity(float *bz, int *dims
 //  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)
  
 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,  int computeSumAbsPerPolarity(float *bz, float *jz_smooth, 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 729  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 (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                  if isnan(bz[j * nx + i]) continue;
                 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) sum1 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                   if (bz[j * nx + i] <= 0) sum2 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
           }           }
  
Line 696  int computeSumAbsPerPolarity(float *bz,
Line 741  int computeSumAbsPerPolarity(float *bz,
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 12:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */  /* Example function 13:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
 // automatically yields erg per cubic centimeter for an input B in Gauss. // automatically yields erg per cubic centimeter for an input B in Gauss.
 // //
Line 708  int computeSumAbsPerPolarity(float *bz,
Line 753  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 770  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 737  int computeFreeEnergy(float *bx, float *
Line 784  int computeFreeEnergy(float *bx, float *
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 13:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */  /* Example function 14:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
  
 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 805  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.8
Karen Tian
 Powered by
ViewCVS 0.9.4