(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.1 and 1.19

version 1.1, 2012/08/23 08:38:33 version 1.19, 2013/10/18 23:36:02
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
  
 ===========================================*/ ===========================================*/
 #include <math.h> #include <math.h>
   #include <mkl.h>
  
 #define PI              (M_PI) #define PI              (M_PI)
 #define MUNAUGHT (0.0000012566370614) /* magnetic constant */ #define MUNAUGHT (0.0000012566370614) /* magnetic constant */
Line 47 
Line 63 
 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2 //  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
 //  =(1.30501e15)Gauss*cm^2 //  =(1.30501e15)Gauss*cm^2
  
 //  The disambig mask value selects only the pixels with values of 5 or 7 -- that is,  int computeAbsFlux(float *bz_err, float *bz, int *dims, float *absFlux,
 //  5: pixels for which the radial acute disambiguation solution was chosen                     float *mean_vf_ptr, float *mean_vf_err_ptr, float *count_mask_ptr, int *mask,
 //  7: pixels for which the radial acute and NRWA disambiguation agree                     int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
   
 int computeAbsFlux(float *bz, int *dims, float *absFlux,  
                                    float *mean_vf_ptr, int *mask,  
                                    float cdelt1, double rsun_ref, double rsun_obs)  
  
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
     double sum=0.0;     double sum=0.0;
       double err = 0.0;
     if (nx <= 0 || ny <= 0) return 1;  
   
     *absFlux = 0.0;     *absFlux = 0.0;
     *mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;
  
         for (j = 0; j < ny; j++)  
         {      if (nx <= 0 || ny <= 0) return 1;
   
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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                    //printf("i,j,bz[j * nx + i]=%d,%d,%f\n",i,j,bz[j * nx + i]);
                     err += bz_err[j * nx + i]*bz_err[j * nx + i];
                   count_mask++;                   count_mask++;
                 }                 }
         }         }
  
      printf("nx=%d,ny=%d,count_mask=%d,sum=%f\n",nx,ny,count_mask,sum);  
      printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);  
      *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;      *mean_vf_ptr = sum*cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
        *mean_vf_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0); // error in the unsigned flux
        *count_mask_ptr  = count_mask;
        //printf("cdelt1=%f\n",cdelt1);
        //printf("rsun_obs=%f\n",rsun_obs);
        //printf("rsun_ref=%f\n",rsun_ref);
        //printf("CMASK=%g\n",*count_mask_ptr);
        //printf("USFLUX=%g\n",*mean_vf_ptr);
        //printf("sum=%f\n",sum);
        //printf("USFLUX_err=%g\n",*mean_vf_err_ptr);
      return 0;      return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 2: Calculate Bh in units of Gauss */  /* Example function 2: Calculate Bh, the horizontal field, in units of Gauss */
 // 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_err, float *by_err, float *bh_err, float *bx, float *by, float *bz, float *bh, int *dims,
                           float *mean_hf_ptr, int *mask)                            float *mean_hf_ptr, int *mask, int *bitmask)
  
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
     float sum=0.0;      int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
     *mean_hf_ptr =0.0;     *mean_hf_ptr =0.0;
  
     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];
                   bh_err[j * nx + i]=sqrt( bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i])/ bh[j * nx + i];
                 count_mask++;                 count_mask++;
               }               }
           }           }
  
     *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram     *mean_hf_ptr = sum/(count_mask); // would be divided by nx*ny if shape of count_mask = shape of magnetogram
     printf("*mean_hf_ptr=%f\n",*mean_hf_ptr);  
     return 0;     return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 3: Calculate Gamma in units of degrees */ /* Example function 3: Calculate Gamma in units of degrees */
 // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI) // Native units of atan(x) are in radians; to convert from radians to degrees, multiply by (180./PI)
   // Redo calculation in radians for error analysis (since derivatives are only true in units of radians).
  
   int computeGamma(float *bz_err, float *bh_err, 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, float *mean_gamma_err_ptr, int *mask, int *bitmask)
                                  float *mean_gamma_ptr, int *mask)  
 { {
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       double err = 0.0;
       *mean_gamma_ptr = 0.0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
     *mean_gamma_ptr=0.0;  
     float sum=0.0;  
     int count=0;  
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 if (bh[j * nx + i] > 100)                 if (bh[j * nx + i] > 100)
                   {                   {
                     //if (mask[j * nx + i] != 90 ) continue;                      if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                     if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;                      if isnan(bz[j * nx + i]) continue;
                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));                      if isnan(bz_err[j * nx + i]) continue;
                       if isnan(bh_err[j * nx + i]) continue;
                       if (bz[j * nx + i] == 0) continue;
                       sum += fabs(atan(bh[j * nx + i]/fabs(bz[j * nx + i])))*(180./PI);
                       err += (( sqrt ( ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i])) + ((bh_err[j * nx + i]*bh_err[j * nx + i])/(bh[j * nx + i]*bh[j * nx + i])))  * fabs(bz[j * nx + i]/bh[j * nx + i]) ) / (1 + (bz[j * nx + i]/bh[j * nx + i])*(bz[j * nx + i]/bh[j * nx + i]))) *(180./PI);
                     count_mask++;                     count_mask++;
                   }                   }
               }               }
           }           }
  
      *mean_gamma_ptr = sum/count_mask;      *mean_gamma_ptr = sum/count_mask;
      printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);       *mean_gamma_err_ptr = (sqrt(err*err))/(count_mask*100.0); // error in the quantity (sum)/(count_mask)
        //printf("MEANGAM=%f\n",*mean_gamma_ptr);
        //printf("MEANGAM_err=%f\n",*mean_gamma_err_ptr);
      return 0;      return 0;
 } }
  
Line 154  int computeGamma(float *bx, float *by, f
Line 193  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_err, float *by_err, float *bz_err, float *bt_err, 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];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
       int j = 0;
       int count_mask = 0;
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
Line 166  int computeB_total(float *bx, float *by,
Line 208  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]);
                   bt_err[j * nx + i]=sqrt(bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i] + by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]*bz_err[j * nx + i]*bz_err[j * nx + i] )/bt[j * nx + i];
               }               }
           }           }
      return 0;      return 0;
Line 175  int computeB_total(float *bx, float *by,
Line 221  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, float *bt_err, float *mean_derivative_btotal_err_ptr)
 { {
  
     int nx = dims[0], ny = dims[1];      int nx = dims[0];
     int i, j, count_mask=0;      int ny = dims[1];
       int i = 0;
     if (nx <= 0 || ny <= 0) return 1;      int j = 0;
       int count_mask = 0;
       double sum = 0.0;
       double err = 0.0;
     *mean_derivative_btotal_ptr = 0.0;     *mean_derivative_btotal_ptr = 0.0;
     float sum = 0.0;  
  
       if (nx <= 0 || ny <= 0) return 1;
  
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
Line 232  int computeBtotalderivative(float *bt, i
Line 280  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] != 90 ) continue;                 if isnan(derx_bt[j * nx + i]) continue;
                if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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 += (2.0)*bt_err[j * nx + i]*bt_err[j * nx + i];
                count_mask++;                count_mask++;
             }             }
           }           }
  
         *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram         *mean_derivative_btotal_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
         printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);          *mean_derivative_btotal_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
           //printf("MEANGBT=%f\n",*mean_derivative_btotal_ptr);
           //printf("MEANGBT_err=%f\n",*mean_derivative_btotal_err_ptr);
         return 0;         return 0;
 } }
  
Line 267  int computeBtotalderivative(float *bt, i
Line 304  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, 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 nx = dims[0], ny = dims[1];       int nx = dims[0];
         int i, j, count_mask=0;       int ny = dims[1];
        int i = 0;
        int j = 0;
        int count_mask = 0;
        double sum= 0.0;
        double err =0.0;
        *mean_derivative_bh_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_derivative_bh_ptr = 0.0;  
         float sum = 0.0;  
   
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
Line 323  int computeBhderivative(float *bh, int *
Line 363  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] != 90 ) continue;                 if isnan(derx_bh[j * nx + i]) continue;
                if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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 += (2.0)*bh_err[j * nx + i]*bh_err[j * nx + i];
                count_mask++;                count_mask++;
             }             }
           }           }
  
         *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram         *mean_derivative_bh_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
           *mean_derivative_bh_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
           //printf("MEANGBH=%f\n",*mean_derivative_bh_ptr);
           //printf("MEANGBH_err=%f\n",*mean_derivative_bh_err_ptr);
   
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* 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, 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 nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double err = 0.0;
           *mean_derivative_bz_ptr = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_derivative_bz_ptr = 0.0;  
         float sum = 0.0;  
   
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
             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 381  int computeBzderivative(float *bz, int *
Line 416  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 390  int computeBzderivative(float *bz, int *
Line 426  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] != 90 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;                 if isnan(bz[j * nx + i]) continue;
                  //if isnan(bz_err[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 */
                  err += 2.0*bz_err[j * nx + i]*bz_err[j * nx + i];
                count_mask++;                count_mask++;
             }             }
           }           }
  
         *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram         *mean_derivative_bz_ptr = (sum)/(count_mask); // would be divided by ((nx-2)*(ny-2)) if shape of count_mask = shape of magnetogram
           *mean_derivative_bz_err_ptr = (sqrt(err))/(count_mask); // error in the quantity (sum)/(count_mask)
           //printf("MEANGBZ=%f\n",*mean_derivative_bz_ptr);
           //printf("MEANGBZ_err=%f\n",*mean_derivative_bz_err_ptr);
   
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
   
 /* 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 460  int computeBzderivative(float *bz, int *
Line 493  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,  
                           float *mean_jz_ptr, float *us_i_ptr, int *mask,  
                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)  
  
   // Comment out random number generator, which can only run on solar3
   //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, float *noisebx,
   //              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 *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;          int nx = dims[0];
           int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_jz_ptr = 0.0;          /* Calculate the derivative*/
         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;  
   
   
         /* brute force method of calculating the derivative (no consideration for edges) */         /* brute force method of calculating the derivative (no consideration for edges) */
   
         for (i = 1; i <= nx-2; i++)         for (i = 1; i <= nx-2; i++)
           {           {
             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;
               }               }
           }           }
  
         /* brute force method of calculating the derivative (no consideration for edges) */  
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             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;
               }               }
           }           }
  
           // consider the edges
         /* consider the edges */  
         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 = 1; i <= nx-2; i++)
         for (i = 0; i < nx; i++)  
           {           {
              for (j = 0; j < ny; j++)              for (j = 1; j <= ny-2; j++)
               {               {
               printf("j=%d\n",j);                 // calculate jz at all points
               printf("i=%d\n",i);  
               printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);                 jz[j * nx + i]            = (derx[j * nx + i]-dery[j * nx + i]);       // jz is in units of Gauss/pix
               printf("derx[j*nx+i]=%f\n",derx[j*nx+i]);                 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]) +
               printf("bx[j*nx+i]=%f\n",bx[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)]) ) ;
               printf("by[j*nx+i]=%f\n",by[j*nx+i]);                 jz_err_squared[j * nx + i]= (jz_err[j * nx + i]*jz_err[j * nx + i]);
                  count_mask++;
   
               }
               }               }
           return 0;
           }           }
         */  
  
   /*===========================================*/
   
   
   /* Example function 9:  Compute quantities on Jz array */
   // Compute mean and total current on Jz array.
   
   int computeJzsmooth(float *bx, float *by, int *dims, float *jz, float *jz_smooth, float *jz_err, float *jz_rms_err, float *jz_err_squared_smooth,
                       float *mean_jz_ptr, float *mean_jz_err_ptr, float *us_i_ptr, float *us_i_err_ptr, int *mask, int *bitmask,
                       float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
  
   {
   
           int nx = dims[0];
           int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double curl = 0.0;
           double us_i = 0.0;
           double err = 0.0;
   
           if (nx <= 0 || ny <= 0) return 1;
   
           /* 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] != 7) && (mask[j * nx + i] != 5) ) continue;                 if isnan(derx[j * nx + i]) continue;
                //if (mask[j * nx + i] != 90 ) 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 */                 if isnan(jz[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 */                 curl +=     (jz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                 us_i += fabs(jz[j * nx + i])*(cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A */
                  err  += (jz_err[j * nx + i]*jz_err[j * nx + i]);
                count_mask++;                count_mask++;
             }             }
           }           }
  
         mean_curl        = (curl/count_mask);          /* Calculate mean vertical current density (mean_jz) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
         printf("mean_curl=%f\n",mean_curl);  
         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);          *mean_jz_err_ptr = (sqrt(err))*fabs(((rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.))/(count_mask)); // error in the quantity MEANJZD
         *us_i_ptr        = (us_i);                   /* us_i gets populated as MEANJZD */  
           *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
           *us_i_err_ptr    = (sqrt(err))*fabs((cdelt1/1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT)); // error in the quantity TOTUSJZ
   
           //printf("MEANJZD=%f\n",*mean_jz_ptr);
           //printf("MEANJZD_err=%f\n",*mean_jz_err_ptr);
   
           //printf("TOTUSJZ=%g\n",*us_i_ptr);
           //printf("TOTUSJZ_err=%g\n",*us_i_err_ptr);
   
         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 weighted by Bz (following Hagino et al., http://adsabs.harvard.edu/abs/2004PASJ...56..831H):
   
          // numerator   = sum of all Jz*Bz
          // denominator = sum of Bz*Bz
          // alpha       = numerator/denominator
   
   // 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 *jz_err, float *bz_err, float *bz, int *dims, float *jz, float *jz_smooth, float *mean_alpha_ptr, float *mean_alpha_err_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];
         int i, j, count_mask=0;          int ny                     = dims[1];
           int i                      = 0;
           int j                      = 0;
           double alpha_total         = 0.0;
           double C                   = ((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.));
           double total               = 0.0;
           double A                   = 0.0;
           double B                   = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
           for (i = 1; i < nx-1; i++)
         *mean_alpha_ptr = 0.0;            {
         float aa, bb, cc, bznew, alpha2, sum=0.0;              for (j = 1; j < ny-1; j++)
                 {
                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(jz[j * nx + i])   continue;
                   if isnan(bz[j * nx + i])   continue;
                   if (jz[j * nx + i] == 0.0) continue;
                   if (bz[j * nx + i] == 0.0) continue;
                   //if (jz_err[j * nx + i] > abs(jz[j * nx + i]) ) continue;
                   //if (bz_err[j * nx + i] > abs(bz[j * nx + i]) ) continue;
                   A += jz[j*nx+i]*bz[j*nx+i];
                   B += bz[j*nx+i]*bz[j*nx+i];
                 }
               }
  
         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] != 90 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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 (jz[j * nx + i] == 0.0) 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 */                  total += bz[j*nx+i]*bz[j*nx+i]*jz_err[j*nx+i]*jz_err[j*nx+i] + (jz[j*nx+i]-2*bz[j*nx+i]*A/B)*(jz[j*nx+i]-2*bz[j*nx+i]*A/B)*bz_err[j*nx+i]*bz_err[j*nx+i];
                 count_mask++;  
               }               }
           }           }
  
         printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);          /* Determine the absolute value of alpha. The units for alpha are 1/Mm */
         printf("count_mask=%d\n",count_mask);          alpha_total              = ((A/B)*C);
         printf("sum=%f\n",sum);          *mean_alpha_ptr          = alpha_total;
         *mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */          *mean_alpha_err_ptr      = (C/B)*(sqrt(total));
   
         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, total unsigned current helicity, absolute value of net 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.
 //  Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/m)  //  Therefore, the units of Bz*Jz = (Gauss)*(Gauss/pix) = (Gauss^2/pix)(pix/arcsec)(arcsec/meter)
 //                                                      = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF) //                                                      = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)
 //                                                      = G^2 / m. //                                                      = G^2 / m.
  
   int computeHelicity(float *jz_err, float *jz_rms_err, float *bz_err, float *bz, int *dims, float *jz, float *mean_ih_ptr,
 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,                      float *mean_ih_err_ptr, float *total_us_ih_ptr, float *total_abs_ih_ptr,
                                         float *total_abs_ih_ptr, int *mask, float cdelt1, double rsun_ref, double rsun_obs)                      float *total_us_ih_err_ptr, float *total_abs_ih_err_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
  
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
           int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double sum2 = 0.0;
           double sum_err = 0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_ih_ptr = 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] != 90 ) continue;                  for (j = 0; j < ny; j++)
                 if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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;
                 if (jz[j * nx + i] == 0.0) continue;                 if (jz[j * nx + i] == 0.0) continue;
                 sum  +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                    sum     +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref); // contributes to MEANJZH and ABSNJZH
                 sum2 += fabs(jz[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); // contributes to TOTUSJH
                     sum_err += sqrt(((jz_err[j * nx + i]*jz_err[j * nx + i])/(jz[j * nx + i]*jz[j * nx + i])) + ((bz_err[j * nx + i]*bz_err[j * nx + i])/(bz[j * nx + i]*bz[j * nx + i]))) * fabs(jz[j * nx + i]*bz[j * nx + i]*(1/cdelt1)*(rsun_obs/rsun_ref));
                 count_mask++;                 count_mask++;
                 }                 }
          }          }
  
             printf("sum/count_mask=%f\n",sum/count_mask);  
             printf("(1/cdelt1)*(rsun_obs/rsun_ref)=%f\n",(1/cdelt1)*(rsun_obs/rsun_ref));  
             *mean_ih_ptr     = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */             *mean_ih_ptr     = sum/count_mask; /* Units are G^2 / m ; keyword is MEANJZH */
             *total_us_ih_ptr = sum2;           /* Units are G^2 / m */          *total_us_ih_ptr      = sum2           ; /* Units are G^2 / m ; keyword is TOTUSJH */
             *total_abs_ih_ptr= fabs(sum);      /* Units are G^2 / m */          *total_abs_ih_ptr     = fabs(sum)      ; /* Units are G^2 / m ; keyword is ABSNJZH */
   
           *mean_ih_err_ptr      = (sqrt(sum_err*sum_err)) / (count_mask*100.0)    ;  // error in the quantity MEANJZH
           *total_us_ih_err_ptr  = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity TOTUSJH
           *total_abs_ih_err_ptr = (sqrt(sum_err*sum_err)) / (100.0)               ;  // error in the quantity ABSNJZH
   
           //printf("MEANJZH=%f\n",*mean_ih_ptr);
           //printf("MEANJZH_err=%f\n",*mean_ih_err_ptr);
   
           //printf("TOTUSJH=%f\n",*total_us_ih_ptr);
           //printf("TOTUSJH_err=%f\n",*total_us_ih_err_ptr);
   
           //printf("ABSNJZH=%f\n",*total_abs_ih_ptr);
           //printf("ABSNJZH_err=%f\n",*total_abs_ih_err_ptr);
  
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* 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.
 //  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)
   //
   //  The error in this quantity is the same as the error in the mean vertical current (mean_jz_err).
  
 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,  int computeSumAbsPerPolarity(float *jz_err, float *bz_err, float *bz, float *jz, int *dims, float *totaljzptr, float *totaljz_err_ptr,
                                                          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];
         int i, j, count_mask=0;          int ny = dims[1];
           int i=0;
           int j=0;
           int count_mask=0;
           double sum1=0.0;
           double sum2=0.0;
           double err=0.0;
           *totaljzptr=0.0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *totaljzptr=0.0;  
         float sum1=0.0, sum2=0.0;  
   
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                 //if (mask[j * nx + i] != 90 ) continue;                  if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;                  if isnan(bz[j * nx + i]) continue;
                   if isnan(jz[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);
                   err += (jz_err[j * nx + i]*jz_err[j * nx + i]);
                   count_mask++;
               }               }
           }           }
  
         *totaljzptr = fabs(sum1) + fabs(sum2);  /* Units are A */         *totaljzptr = fabs(sum1) + fabs(sum2);  /* Units are A */
           *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
           //printf("SAVNCPP=%g\n",*totaljzptr);
           //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
   
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* 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. Note that the 8*PI can come out of the integral; thus,
   // the integral is over B^2 dV and the 8*PI is divided at the end.
 // //
 // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert // Total magnetic energy is the magnetic energy density times dA, or the area, and the units are thus ergs/cm. To convert
 // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm: // ergs per centimeter cubed to ergs per centimeter, simply multiply by the area per pixel in cm:
 // erg/cm^3(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.)^2  //   erg/cm^3*(CDELT1^2)*(RSUN_REF/RSUN_OBS ^2)*(100.^2)
 // = erg/cm^3(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2  // = erg/cm^3*(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
 // = erg/cm^3(1.30501e15)  // = erg/cm^3*(1.30501e15)
 // = 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_err, float *by_err, float *bx, float *by, float *bpx, float *bpy, int *dims,
                                           float *meanpotptr, float *totpotptr, int *mask,                        float *meanpotptr, float *meanpot_err_ptr, float *totpotptr, float *totpot_err_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];
         int i, j, count_mask=0;          int ny = dims[1];
           int i = 0;
         if (nx <= 0 || ny <= 0) return 1;          int j = 0;
           int count_mask = 0;
           double sum = 0.0;
           double sum1 = 0.0;
           double err = 0.0;
         *totpotptr=0.0;         *totpotptr=0.0;
         *meanpotptr=0.0;         *meanpotptr=0.0;
         float sum=0.0;  
           if (nx <= 0 || ny <= 0) return 1;
  
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  //if (mask[j * nx + i] != 90 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) continue;                   if isnan(bx[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);                   if isnan(by[j * nx + i]) continue;
                    sum  += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
                    sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );
                    err  += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
                  count_mask++;                  count_mask++;
               }               }
           }           }
  
         *meanpotptr = (sum)/(count_mask);              /* Units are ergs per cubic centimeter */          *meanpotptr      = (sum1/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */
         *totpotptr  = sum*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0)*(count_mask);   /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2, units of count_mask are pix^2; therefore, units of totpotptr are ergs per centimeter */          *meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
   
           /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
           *totpotptr       = (sum)/(8.*PI);
           *totpot_err_ptr  = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0*(1/(8.*PI)));
   
           //printf("MEANPOT=%g\n",*meanpotptr);
           //printf("MEANPOT_err=%g\n",*meanpot_err_ptr);
   
           //printf("TOTPOT=%g\n",*totpotptr);
           //printf("TOTPOT_err=%g\n",*totpot_err_ptr);
   
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* 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_err, float *by_err, float *bh_err, 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 *meanshear_angle_err_ptr, float *area_w_shear_gt_45ptr, int *mask, int *bitmask)
                                           float *meanshear_anglehptr, float *area_w_shear_gt_45hptr,  
                                           int *mask)  
 { {
         int nx = dims[0], ny = dims[1];          int nx = dims[0];
         int i, j;          int ny = dims[1];
           int i = 0;
         if (nx <= 0 || ny <= 0) return 1;          int j = 0;
           int count_mask = 0;
           double dotproduct = 0.0;
           double magnitude_potential = 0.0;
           double magnitude_vector = 0.0;
           double shear_angle = 0.0;
           double err = 0.0;
           double sum = 0.0;
           double count = 0.0;
         *area_w_shear_gt_45ptr=0.0;         *area_w_shear_gt_45ptr=0.0;
         *meanshear_angleptr=0.0;         *meanshear_angleptr=0.0;
         float dotproduct, magnitude_potential, magnitude_vector, shear_angle=0.0, sum = 0.0, count=0.0, count_mask=0.0;  
         float dotproducth, magnitude_potentialh, magnitude_vectorh, shear_angleh=0.0, sum1 = 0.0, counth = 0.0;          if (nx <= 0 || ny <= 0) return 1;
  
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                  //if (mask[j * nx + i] != 90 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                  if ( (mask[j * nx + i] != 7) && (mask[j * nx + i] != 5) ) 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]));
Line 780  int computeShearAngle(float *bx, float *
Line 931  int computeShearAngle(float *bx, float *
                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);                  shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
                  count ++;                  count ++;
                  sum += shear_angle ;                  sum += shear_angle ;
                    err += -(1./(1.- sqrt(bx_err[j * nx + i]*bx_err[j * nx + i]+by_err[j * nx + i]*by_err[j * nx + i]+bh_err[j * nx + i]*bh_err[j * nx + i])));
                  if (shear_angle > 45) count_mask ++;                  if (shear_angle > 45) count_mask ++;
               }               }
           }           }
  
         /* For mean 3D shear angle, area with shear greater than 45*/         /* For mean 3D shear angle, area with shear greater than 45*/
         *meanshear_angleptr = (sum)/(count);              /* Units are degrees */         *meanshear_angleptr = (sum)/(count);              /* Units are degrees */
         printf("count_mask=%f\n",count_mask);          *meanshear_angle_err_ptr = (sqrt(err*err))/(count);  // error in the quantity (sum)/(count_mask)
         printf("count=%f\n",count);          *area_w_shear_gt_45ptr   = (count_mask/(count))*(100.0);/* The area here is a fractional area -- the % of the total area */
         *area_w_shear_gt_45ptr = (count_mask/(count))*(100.);  /* The area here is a fractional area -- the % of the total area */  
           //printf("MEANSHR=%f\n",*meanshear_angleptr);
           //printf("MEANSHR_err=%f\n",*meanshear_angle_err_ptr);
  
         return 0;         return 0;
 } }
Line 909  void greenpot(float *bx, float *by, floa
Line 1063  void greenpot(float *bx, float *by, floa
  
  
 /*===========END OF KEIJI'S CODE =========================*/ /*===========END OF KEIJI'S CODE =========================*/
   
   char *sw_functions_version() // Returns CVS version of sw_functions.c
   {
     return strdup("$Id$");
   }
   
 /* ---------------- end of this file ----------------*/ /* ---------------- end of this file ----------------*/

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