(file) Return to sw_functions.c CVS log (file) (dir) Up to [Development] / JSOC / proj / sharp / apps

File: [Development] / JSOC / proj / sharp / apps / sw_functions.c (download)
Revision: 1.6, Tue Feb 5 18:45:55 2013 UTC (9 years, 5 months ago) by mbobra
Branch: MAIN
Changes since 1.5: +15 -0 lines
Only changed documentation

/*===========================================

   The following 14 functions calculate the following spaceweather indices:

    USFLUX Total unsigned flux in Maxwells
    MEANGAM Mean inclination angle, gamma, in degrees
    MEANGBT Mean value of the total field gradient, in Gauss/Mm
    MEANGBZ Mean value of the vertical field gradient, in Gauss/Mm
    MEANGBH Mean value of the horizontal field gradient, in Gauss/Mm
    MEANJZD Mean vertical current density, in mA/m2
    TOTUSJZ Total unsigned vertical current, in Amperes
    MEANALP Mean twist parameter, alpha, in 1/Mm
    MEANJZH Mean current helicity in G2/m
    TOTUSJH Total unsigned current helicity in G2/m
    ABSNJZH Absolute value of the net current helicity in G2/m
    SAVNCPP Sum of the Absolute Value of the Net Currents Per Polarity in Amperes
    MEANPOT Mean photospheric excess 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  

   The indices are calculated on the pixels in which the conf_disambig segment is greater than 70 and
   pixels in which the bitmap segment is greater than 30. These ranges are selected because the CCD 
   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 
   Potential Field code (appended) written by Keiji Hayashi

===========================================*/
#include <math.h>

#define PI              (M_PI)
#define MUNAUGHT (0.0000012566370614) /* magnetic constant */

/*===========================================*/

/* Example function 1: Compute total unsigned flux in units of G/cm^2 */

//  To compute the unsigned flux, we simply calculate 
//  flux = surface integral [(vector Bz) dot (normal vector)],
//       = surface integral [(magnitude Bz)*(magnitude normal)*(cos theta)].
//  However, since the field is radial, we will assume cos theta = 1.
//  Therefore the pixels only need to be corrected for the projection.

//  To convert G to G*cm^2, simply multiply by the number of square centimeters per pixel.
//  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
//  (Gauss/pix^2)(CDELT1)^2(RSUN_REF/RSUN_OBS)^2(100.cm/m)^2
//  =(Gauss/pix^2)(0.5 arcsec/pix)^2(722500m/arcsec)^2(100cm/m)^2
//  =(1.30501e15)Gauss*cm^2

//  The disambig mask value selects only the pixels with values of 5 or 7 -- that is,
//  5: pixels for which the radial acute disambiguation solution was chosen
//  7: pixels for which the radial acute and NRWA disambiguation agree

int computeAbsFlux(float *bz, int *dims, float *absFlux, 
				   float *mean_vf_ptr, int *mask,  int *bitmask,
                                   float cdelt1, double rsun_ref, double rsun_obs)

{

    int nx = dims[0], ny = dims[1];
    int i, j, count_mask=0;
    double sum=0.0;
    
    if (nx <= 0 || ny <= 0) return 1;
	
    *absFlux = 0.0;
    *mean_vf_ptr =0.0;

	for (i = 0; i < nx; i++) 
	{
		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]));
                  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;
     return 0;
}

/*===========================================*/
/* Example function 2: Calculate Bh in units of Gauss */
// Native units of Bh are Gauss

int computeBh(float *bx, float *by, float *bz, float *bh, int *dims, 
			  float *mean_hf_ptr, int *mask, int *bitmask)

{

    int nx = dims[0], ny = dims[1];
    int i, j, count_mask=0;
    float sum=0.0;	
    *mean_hf_ptr =0.0;

    if (nx <= 0 || ny <= 0) return 1;

	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] );
                sum += bh[j * nx + i];
                count_mask++;
	      }	
	  }
     
    *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;
}

/*===========================================*/
/* 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)


int computeGamma(float *bx, float *by, float *bz, float *bh, int *dims,
				 float *mean_gamma_ptr, int *mask, int *bitmask)
{
    int nx = dims[0], ny = dims[1];
    int i, j, count_mask=0;

    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 (j = 0; j < ny; j++) 
	      {
		if (bh[j * nx + i] > 100)
		  {
                    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));
		    count_mask++;
		  }
	      }
	  }

     *mean_gamma_ptr = sum/count_mask;
     printf("*mean_gamma_ptr=%f\n",*mean_gamma_ptr);
     return 0;
}

/*===========================================*/
/* Example function 4: Calculate B_Total*/
// Native units of B_Total are in gauss

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 i, j, count_mask=0;
	
    if (nx <= 0 || ny <= 0) return 1;

	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;
                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]);
	      }	
	  }
     return 0;
}

/*===========================================*/
/* Example function 5:  Derivative of B_Total SQRT( (dBt/dx)^2 + (dBt/dy)^2 ) */

int computeBtotalderivative(float *bt, int *dims, float *mean_derivative_btotal_ptr, int *mask, int *bitmask, float *derx_bt, float *dery_bt)
{

    int nx = dims[0], ny = dims[1];
    int i, j, count_mask=0;

    if (nx <= 0 || ny <= 0) return 1;

    *mean_derivative_btotal_ptr = 0.0;
    float sum = 0.0;


        /* brute force method of calculating the derivative (no consideration for edges) */
      	for (i = 1; i <= nx-2; i++) 
	  {
	    for (j = 0; j <= ny-1; j++) 
	      {
		derx_bt[j * nx + i] = (bt[j * nx + i+1] - bt[j * nx + i-1])*0.5;
              }
          }

        /* brute force method of calculating the derivative (no consideration for edges) */
      	for (i = 0; i <= nx-1; i++) 
	  {
	    for (j = 1; j <= ny-2; j++) 
	      {
                dery_bt[j * nx + i] = (bt[(j+1) * nx + i] - bt[(j-1) * nx + i])*0.5;
              }
          }


        /* consider the edges */
        i=0; 
      	for (j = 0; j <= ny-1; j++) 
          {
             derx_bt[j * nx + i] = ( (-3*bt[j * nx + i]) + (4*bt[j * nx + (i+1)]) - (bt[j * nx + (i+2)]) )*0.5;
          }

        i=nx-1; 
      	for (j = 0; j <= ny-1; j++) 
          {
             derx_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[j * nx + (i-1)]) - (-bt[j * nx + (i-2)]) )*0.5; 
          }

        j=0;
        for (i = 0; i <= nx-1; i++) 
          {
             dery_bt[j * nx + i] = ( (-3*bt[j*nx + i]) + (4*bt[(j+1) * nx + i]) - (bt[(j+2) * nx + i]) )*0.5; 
          }

        j=ny-1;
        for (i = 0; i <= nx-1; i++) 
          {
             dery_bt[j * nx + i] = ( (3*bt[j * nx + i]) + (-4*bt[(j-1) * nx + i]) - (-bt[(j-2) * nx + i]) )*0.5;
          }


      	for (i = 0; i <= nx-1; i++) 
          {
            for (j = 0; j <= ny-1; j++) 
            {
	       if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */
               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
        printf("*mean_derivative_btotal_ptr=%f\n",*mean_derivative_btotal_ptr);
        return 0;
}


/*===========================================*/
/* 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, int *bitmask, float *derx_bh, float *dery_bh)
{

        int nx = dims[0], ny = dims[1];
        int i, j, count_mask=0;

        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) */
      	for (i = 1; i <= nx-2; i++) 
	  {
	    for (j = 0; j <= ny-1; j++) 
	      {
		derx_bh[j * nx + i] = (bh[j * nx + i+1] - bh[j * nx + i-1])*0.5;
              }
          }

        /* brute force method of calculating the derivative (no consideration for edges) */
      	for (i = 0; i <= nx-1; i++) 
	  {
	    for (j = 1; j <= ny-2; j++) 
	      {
                dery_bh[j * nx + i] = (bh[(j+1) * nx + i] - bh[(j-1) * nx + i])*0.5;
              }
          }


        /* consider the edges */
        i=0; 
      	for (j = 0; j <= ny-1; j++) 
          {
             derx_bh[j * nx + i] = ( (-3*bh[j * nx + i]) + (4*bh[j * nx + (i+1)]) - (bh[j * nx + (i+2)]) )*0.5;
          }

        i=nx-1; 
      	for (j = 0; j <= ny-1; j++) 
          {
             derx_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[j * nx + (i-1)]) - (-bh[j * nx + (i-2)]) )*0.5; 
          }

        j=0;
        for (i = 0; i <= nx-1; i++) 
          {
             dery_bh[j * nx + i] = ( (-3*bh[j*nx + i]) + (4*bh[(j+1) * nx + i]) - (bh[(j+2) * nx + i]) )*0.5; 
          }

        j=ny-1;
        for (i = 0; i <= nx-1; i++) 
          {
             dery_bh[j * nx + i] = ( (3*bh[j * nx + i]) + (-4*bh[(j-1) * nx + i]) - (-bh[(j-2) * nx + i]) )*0.5;
          }


      	for (i = 0; i <= nx-1; i++) 
          {
            for (j = 0; j <= ny-1; j++) 
            {
	       if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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 */
               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
        return 0;
}

/*===========================================*/
/* 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, int *bitmask, float *derx_bz, float *dery_bz)
{

	int nx = dims[0], ny = dims[1];
	int i, j, count_mask=0;

	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) */
      	for (i = 1; i <= nx-2; i++) 
	  {
	    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;
              }
          }

        /* brute force method of calculating the derivative (no consideration for edges) */
      	for (i = 0; i <= nx-1; i++) 
	  {
	    for (j = 1; j <= ny-2; j++) 
	      {
                if isnan(bz[j * nx + i]) continue;
                dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
              }
          }


        /* consider the edges */
        i=0; 
      	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;
          }

        i=nx-1; 
      	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; 
          }

        j=0;
        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; 
          }

        j=ny-1;
        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;
          }


      	for (i = 0; i <= nx-1; i++) 
          {
            for (j = 0; j <= ny-1; j++) 
            {
               // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) 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 */
               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
	return 0;
}

/*===========================================*/
/* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */

//  In discretized space like data pixels,
//  the current (or curl of B) is calculated as
//  the integration of the field Bx and By along
//  the circumference of the data pixel divided by the area of the pixel.
//  One form of differencing (a word for the differential operator
//  in the discretized space) of the curl is expressed as 
//  (dx * (Bx(i,j-1)+Bx(i,j)) / 2
//  +dy * (By(i+1,j)+By(i,j)) / 2
//  -dx * (Bx(i,j+1)+Bx(i,j)) / 2
//  -dy * (By(i-1,j)+By(i,j)) / 2) / (dx * dy) 
//
//  
//  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:
//  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
//  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
//  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).
//  In that case, we would have the following:
//  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
//  jz * (35.0)
//
//  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)(0.0010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
//  =(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 *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;
	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) */
      	for (i = 1; i <= nx-2; i++) 
	  {
	    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;
              }
          }

        /* brute force method of calculating the derivative (no consideration for edges) */
      	for (i = 0; i <= nx-1; i++) 
	  {
	    for (j = 1; j <= ny-2; j++) 
	      {
                 if isnan(bx[j * nx + i]) continue;
                 dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
              }
          }


        /* consider the edges */
        i=0; 
      	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;
          }

        i=nx-1; 
      	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;
          }

        j=0;
        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;
          }

        j=ny-1;
        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;
          }


      	for (i = 0; i <= nx-1; i++) 
          {
            for (j = 0; j <= ny-1; j++) 
            {
               /* calculate jz at all points */ 
               jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                              /* jz is in units of Gauss/pix */ 
               count_mask++;
            }	
	  }

	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 (j = 0; j <= ny-1; j++) 
            {
               if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
               if isnan(derx[j * nx + i]) continue;
               if isnan(dery[j * nx + i]) continue;
               if isnan(jz_smooth[j * nx + i]) continue;
               //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])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */
               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);
        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 */
        printf("count_mask=%d\n",count_mask);
        *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
	return 0;

}

/*===========================================*/

/* 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.
//
// 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)
//                               = 1/Mm

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)

{
	int nx = dims[0], ny = dims[1];
	int i, j=0;

	if (nx <= 0 || ny <= 0) return 1;

	*mean_alpha_ptr = 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 (j = 1; j < ny-1; j++) 
	      {
                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 (bz[j * nx + i] == 0.0) continue;
                if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i]);
                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];
	      }	
	  }
        
        sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */

        /* Determine the sign of alpha */
        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;
}

/*===========================================*/
/* 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 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) 
//                                                      = (Gauss^2/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF) 
//                                                      = G^2 / m.


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, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{

	int nx = dims[0], ny = dims[1];
	int i, j, count_mask=0;
	
	if (nx <= 0 || ny <= 0) return 1;

	*mean_ih_ptr = 0.0;
	float sum=0.0, sum2=0.0;

	for (i = 0; i < nx; i++) 
	{
		for (j = 0; j < ny; j++) 
		{
                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 (bz[j * nx + i] == 0.0) continue;
                if (jz_smooth[j * nx + i] == 0.0) continue;
                sum  +=     (jz_smooth[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++;
                }	
	 }

            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 */ 
	    *total_us_ih_ptr = sum2;           /* Units are G^2 / m */
	    *total_abs_ih_ptr= fabs(sum);      /* Units are G^2 / m */

	return 0;
}

/*===========================================*/
/* Example function 12:  Sum of Absolute Value per polarity  */

//  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.
//  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),
//     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)

int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr, 
							 int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)

{	
	int nx = dims[0], ny = dims[1];
	int i, j, count_mask=0;

	if (nx <= 0 || ny <= 0) return 1;
	
	*totaljzptr=0.0;
	float sum1=0.0, sum2=0.0;
     
	for (i = 0; i < nx; i++) 
	  {
	    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;
		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);
       	      }
	  }
	
	*totaljzptr = fabs(sum1) + fabs(sum2);  /* Units are A */
	return 0;
}

/*===========================================*/
/* 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
// automatically yields erg per cubic centimeter for an input B in Gauss.
//
// 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:
// 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(1.30501e15)
// = erg/cm(1/pix^2)

int computeFreeEnergy(float *bx, float *by, float *bpx, float *bpy, int *dims, 
					  float *meanpotptr, float *totpotptr, int *mask, int *bitmask, 
					  float cdelt1, double rsun_ref, double rsun_obs)

{
	int nx = dims[0], ny = dims[1];
	int i, j, count_mask=0;
	
	if (nx <= 0 || ny <= 0) return 1;
	
        *totpotptr=0.0;
	*meanpotptr=0.0;
	float sum=0.0;

	for (i = 0; i < nx; i++) 
	  {
	    for (j = 0; j < ny; j++) 
	      {
                 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);
                 count_mask++;
	      }
	  }

	*meanpotptr = (sum)/(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 */
	return 0;
}

/*===========================================*/
/* 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,
					  float *meanshear_angleptr, float *area_w_shear_gt_45ptr, 
					  float *meanshear_anglehptr, float *area_w_shear_gt_45hptr, 
					  int *mask, int *bitmask)
{	
	int nx = dims[0], ny = dims[1];
	int i, j;
	
	if (nx <= 0 || ny <= 0) return 1;
	
        *area_w_shear_gt_45ptr=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;

	for (i = 0; i < nx; i++) 
	  {
	    for (j = 0; j < ny; j++) 
	      {
                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(bpx[j * nx + i]) continue;                
                 if isnan(bpy[j * nx + i]) continue;                
                 if isnan(bpz[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*/
                 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_vector      = 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]) );
                 shear_angle           = acos(dotproduct/(magnitude_potential*magnitude_vector))*(180./PI);
                 count ++;
                 sum += shear_angle ;
                 if (shear_angle > 45) count_mask ++;
	      }
	  }
	
        /* For mean 3D shear angle, area with shear greater than 45*/
	*meanshear_angleptr = (sum)/(count);              /* Units are degrees */
        printf("count_mask=%f\n",count_mask);
        printf("count=%f\n",count);
        *area_w_shear_gt_45ptr = (count_mask/(count))*(100.);  /* The area here is a fractional area -- the % of the total area */

	return 0;
}


/*==================KEIJI'S CODE =========================*/

// #include <omp.h>
#include <math.h>

void greenpot(float *bx, float *by, float *bz, int nnx, int nny)
{
/* local workings */
  int inx, iny, i, j, n;
/* local array */
  float *pfpot, *rdist;
  pfpot=(float *)malloc(sizeof(float) *nnx*nny);
  rdist=(float *)malloc(sizeof(float) *nnx*nny);
  float *bztmp;
  bztmp=(float *)malloc(sizeof(float) *nnx*nny);
/* make nan */
//  unsigned long long llnan = 0x7ff0000000000000;
//  float NAN = (float)(llnan);

// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){pfpot[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){rdist[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){bx[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++){by[nnx*iny+inx] = 0.0;}}
// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
  {
    float val0 = bz[nnx*iny + inx];
    if (isnan(val0)){bztmp[nnx*iny + inx] = 0.0;}else{bztmp[nnx*iny + inx] = val0;}
  }}

  // dz is the monopole depth
  float dz = 0.001;

// #pragma omp parallel for private (inx)
  for (iny=0; iny < nny; iny++){for (inx=0; inx < nnx; inx++)
  {
    float rdd, rdd1, rdd2;
    float r;
    rdd1 = (float)(inx);
    rdd2 = (float)(iny);
    rdd = rdd1 * rdd1 + rdd2 * rdd2 + dz * dz;
    rdist[nnx*iny+inx] = 1.0/sqrt(rdd);
  }}

  int iwindow;
  if (nnx > nny) {iwindow = nnx;} else {iwindow = nny;}
  float rwindow;
  rwindow = (float)(iwindow);
  rwindow = rwindow * rwindow + 0.01; // must be of square

  rwindow = 1.0e2; // limit the window size to be 10.

  rwindow = sqrt(rwindow);
  iwindow = (int)(rwindow);

// #pragma omp parallel for private(inx)
  for (iny=0;iny<nny;iny++){for (inx=0;inx<nnx;inx++)
  {
    float val0 = bz[nnx*iny + inx];
    if (isnan(val0))
    {
      pfpot[nnx*iny + inx] = 0.0; // hmmm.. NAN;
    }
    else
    {
      float sum;
      sum = 0.0;
      int j2, i2;
      int j2s, j2e, i2s, i2e;
      j2s = iny - iwindow;
      j2e = iny + iwindow;
      if (j2s <   0){j2s =   0;}
      if (j2e > nny){j2e = nny;}
      i2s = inx - iwindow;
      i2e = inx + iwindow;
      if (i2s <   0){i2s =   0;}
      if (i2e > nnx){i2e = nnx;}

      for (j2=j2s;j2<j2e;j2++){for (i2=i2s;i2<i2e;i2++)
      {
        float val1 = bztmp[nnx*j2 + i2];
        float rr, r1, r2;
//        r1 = (float)(i2 - inx);
//        r2 = (float)(j2 - iny);
//        rr = r1*r1 + r2*r2;
//        if (rr < rwindow)
//        {
          int   di, dj;
          di = abs(i2 - inx);
          dj = abs(j2 - iny);
          sum = sum + val1 * rdist[nnx * dj + di] * dz;
//        }
      } }
      pfpot[nnx*iny + inx] = sum; // Note that this is a simplified definition.
    }
  } } // end of OpenMP parallelism

// #pragma omp parallel for private(inx)
  for (iny=1; iny < nny - 1; iny++){for (inx=1; inx < nnx - 1; inx++)
  {
    bx[nnx*iny + inx] = -(pfpot[nnx*iny + (inx+1)]-pfpot[nnx*iny + (inx-1)]) * 0.5;
    by[nnx*iny + inx] = -(pfpot[nnx*(iny+1) + inx]-pfpot[nnx*(iny-1) + inx]) * 0.5;
  } } // end of OpenMP parallelism

  free(rdist);
  free(pfpot);
  free(bztmp);
} // end of void func. greenpot


/*===========END OF KEIJI'S CODE =========================*/
/* ---------------- end of this file ----------------*/

Karen Tian
Powered by
ViewCVS 0.9.4