(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.4 and 1.7

version 1.4, 2012/12/18 23:05:10 version 1.7, 2013/02/08 23:41:13
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 81  int computeAbsFlux(float *bz, int *dims,
Line 81  int computeAbsFlux(float *bz, int *dims,
     *absFlux = 0.0;     *absFlux = 0.0;
     *mean_vf_ptr =0.0;     *mean_vf_ptr =0.0;
  
         for (j = 0; j < ny; j++)  
         {  
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   for (j = 0; j < ny; j++)
                   {
                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                   if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                   if isnan(bz[j * nx + i]) continue;                   if isnan(bz[j * nx + i]) continue;
                   //printf("bz[j * nx + i]=%f\n",bz[j * nx + i]);  
                   sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
                   count_mask++;                   count_mask++;
                 }                 }
Line 115  int computeBh(float *bx, float *by, floa
Line 114  int computeBh(float *bx, float *by, floa
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
         for (j = 0; j < ny; j++)  
           {  
             for (i = 0; i < nx; i++)             for (i = 0; i < nx; i++)
               {               {
               for (j = 0; j < ny; j++)
                 {
                 if isnan(bx[j * nx + i]) continue;                 if isnan(bx[j * nx + i]) continue;
                 if isnan(by[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] );
Line 253  int computeBtotalderivative(float *bt, i
Line 252  int computeBtotalderivative(float *bt, i
           }           }
  
  
         /* Just some print statements  
         for (i = 0; i < nx; i++)  
           {  
              for (j = 0; j < ny; j++)  
               {  
               printf("j=%d\n",j);  
               printf("i=%d\n",i);  
               printf("dery_bt[j*nx+i]=%f\n",dery_bt[j*nx+i]);  
               printf("derx_bt[j*nx+i]=%f\n",derx_bt[j*nx+i]);  
               printf("bt[j*nx+i]=%f\n",bt[j*nx+i]);  
               }  
           }  
         */  
   
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bt[j * nx + i]-dery_bt[j * nx + i]) == 0) continue;  
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || 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 */                sum += sqrt( derx_bt[j * nx + i]*derx_bt[j * nx + i]  + dery_bt[j * nx + i]*dery_bt[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 343  int computeBhderivative(float *bh, int *
Line 329  int computeBhderivative(float *bh, int *
           }           }
  
  
         /*Just some print statements  
         for (i = 0; i < nx; i++)  
           {  
              for (j = 0; j < ny; j++)  
               {  
               printf("j=%d\n",j);  
               printf("i=%d\n",i);  
               printf("dery_bh[j*nx+i]=%f\n",dery_bh[j*nx+i]);  
               printf("derx_bh[j*nx+i]=%f\n",derx_bh[j*nx+i]);  
               printf("bh[j*nx+i]=%f\n",bh[j*nx+i]);  
               }  
           }  
         */  
   
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bh[j * nx + i]-dery_bh[j * nx + i]) == 0) continue;  
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || 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 */                sum += sqrt( derx_bh[j * nx + i]*derx_bh[j * nx + i]  + dery_bh[j * nx + i]*dery_bh[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 437  int computeBzderivative(float *bz, int *
Line 410  int computeBzderivative(float *bz, 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_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++)
Line 470  int computeBzderivative(float *bz, int *
Line 429  int computeBzderivative(float *bz, int *
 } }
  
 /*===========================================*/ /*===========================================*/
   
 /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */ /* Example function 8:  Current Jz = (dBy/dx) - (dBx/dy) */
  
 //  In discretized space like data pixels, //  In discretized space like data pixels,
Line 498  int computeBzderivative(float *bz, int *
Line 456  int computeBzderivative(float *bz, int *
 // //
 //  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)(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/0.5)(10^-4)(4*PI*10^7)(722500)(1000.)
 //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.) //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
  
 int computeJz(float *bx, float *by, int *dims, float *jz, int computeJz(float *bx, float *by, int *dims, float *jz,
                           float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,                           int *mask, int *bitmask,
                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                           float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
  
  
Line 514  int computeJz(float *bx, float *by, int
Line 471  int computeJz(float *bx, float *by, int
         int i, j, count_mask=0;         int i, j, count_mask=0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
   
         *mean_jz_ptr = 0.0;  
         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;         float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
  
  
Line 574  int computeJz(float *bx, float *by, int
Line 529  int computeJz(float *bx, float *by, int
           {           {
             for (j = 0; j <= ny-1; j++)             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 ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                if isnan(derx[j * nx + i]) continue;                if isnan(derx[j * nx + i]) continue;
                if isnan(dery[j * nx + i]) continue;                if isnan(dery[j * nx + i]) continue;
                curl +=     (derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref)*(0.00010)*(1/MUNAUGHT)*(1000.); /* curl is in units of mA / m^2 */                 if isnan(jz_smooth[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 */                 //printf("%d,%d,%f\n",i,j,jz_smooth[j * nx + i]);
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                 curl +=     (jz_smooth[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(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++;                count_mask++;
             }             }
           }           }
  
           /* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
         mean_curl        = (curl/count_mask);         mean_curl        = (curl/count_mask);
         printf("mean_curl=%f\n",mean_curl);         printf("mean_curl=%f\n",mean_curl);
         printf("cdelt1, what is it?=%f\n",cdelt1);         printf("cdelt1, what is it?=%f\n",cdelt1);
Line 594  int computeJz(float *bx, float *by, int
Line 590  int computeJz(float *bx, float *by, int
  
 } }
  
   
 /*===========================================*/ /*===========================================*/
 /* Example function 9:  Twist Parameter, alpha */  
  
 // The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm  /* Example function 10:  Twist Parameter, alpha */
   
   // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
   // for alpha is calculated in the following way (different from Leka and Barnes' approach):
   
          // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
          // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
          // avg alpha = avg Jz / avg Bz
   
   // The sign is assigned as follows:
   // If the sum of all Bz is greater than 0, then evaluate the sum of Jz at the positive Bz pixels.
   // If this value is > 0, then alpha is > 0.
   // If this value is < 0, then alpha is <0.
   //
   // If the sum of all Bz is less than 0, then evaluate the sum of Jz at the negative Bz pixels.
   // If this value is > 0, then alpha is < 0.
   // If this value is < 0, then alpha is > 0.
   
   // The units of alpha are in 1/Mm
 // The units of Jz are in Gauss/pix; the units of Bz are in Gauss. // The units of Jz are in Gauss/pix; the units of Bz are in Gauss.
 // //
 // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or // Therefore, the units of Jz/Bz = (Gauss/pix)(1/Gauss)(pix/arcsec)(arsec/meter)(meter/Mm), or
 //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6) //                               = (Gauss/pix)(1/Gauss)(1/CDELT1)(RSUN_OBS/RSUN_REF)(10^6)
 //                               = 1/Mm //                               = 1/Mm
  
 int computeAlpha(float *bz, int *dims, float *jz, float *mean_alpha_ptr, int *mask, int *bitmask,  int computeAlpha(float *bz, int *dims, float *jz_smooth, float *mean_alpha_ptr, int *mask, int *bitmask,
                  float cdelt1, double rsun_ref, double rsun_obs)                  float cdelt1, double rsun_ref, double rsun_obs)
   
 { {
         int nx = dims[0], ny = dims[1];         int nx = dims[0], ny = dims[1];
         int i, j, count_mask=0;          int i, j=0;
  
         if (nx <= 0 || ny <= 0) return 1;         if (nx <= 0 || ny <= 0) return 1;
  
         *mean_alpha_ptr = 0.0;         *mean_alpha_ptr = 0.0;
         float aa, bb, cc, bznew, alpha2, sum=0.0;          float aa, bb, cc, bznew, alpha2, sum1, sum2, sum3, sum4, sum, sum5, sum6=0.0;
  
         for (i = 1; i < nx-1; i++)         for (i = 1; i < nx-1; i++)
           {           {
             for (j = 1; j < ny-1; j++)             for (j = 1; j < ny-1; j++)
               {               {
                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                  if isnan(jz_smooth[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
                 sum += (jz[j * nx + i] / bz[j * nx + i])*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */                  if (bz[j * nx + i] >  0) sum1 += ( bz[j * nx + i]);
                 count_mask++;                  if (bz[j * nx + i] <= 0) sum2 += ( bz[j * nx + i]);
                   if (bz[j * nx + i] >  0) sum3 += ( jz_smooth[j * nx + i]);
                   if (bz[j * nx + i] <= 0) sum4 += ( jz_smooth[j * nx + i]);
                   sum5 += bz[j * nx + i];
               }               }
           }           }
  
         printf("cdelt1=%f,rsun_ref=%f,rsun_obs=%f\n",cdelt1,rsun_ref,rsun_obs);          sum = (((fabs(sum3))+(fabs(sum4)))/((fabs(sum2))+sum1))*((1/cdelt1)*(rsun_obs/rsun_ref)*(1000000.)) ; /* the units for (jz/bz) are 1/Mm */
         printf("count_mask=%d\n",count_mask);  
         printf("sum=%f\n",sum);          /* Determine the sign of alpha */
         *mean_alpha_ptr = sum/count_mask; /* Units are 1/Mm */          if ((sum5 > 0) && (sum3 >  0)) sum=sum;
           if ((sum5 > 0) && (sum3 <= 0)) sum=-sum;
           if ((sum5 < 0) && (sum4 <= 0)) sum=sum;
           if ((sum5 < 0) && (sum4 >  0)) sum=-sum;
   
           *mean_alpha_ptr = sum; /* Units are 1/Mm */
         return 0;         return 0;
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 10:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */  /* Example function 11:  Helicity (mean current helicty, mean unsigned current helicity, and mean absolute current helicity) */
  
 //  The current helicity is defined as Bz*Jz and the units are G^2 / m //  The current helicity is defined as Bz*Jz and the units are G^2 / m
 //  The units of Jz are in G/pix; the units of Bz are in G. //  The units of Jz are in G/pix; the units of Bz are in G.
Line 646  int computeAlpha(float *bz, int *dims, f
Line 667  int computeAlpha(float *bz, int *dims, f
 //                                                      = G^2 / m. //                                                      = G^2 / m.
  
  
 int computeHelicity(float *bz, int *dims, float *jz, float *mean_ih_ptr, float *total_us_ih_ptr,  int computeHelicity(float *bz, int *dims, float *jz_smooth, float *mean_ih_ptr, float *total_us_ih_ptr,
                                         float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                                         float *total_abs_ih_ptr, int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 659  int computeHelicity(float *bz, int *dims
Line 680  int computeHelicity(float *bz, int *dims
         *mean_ih_ptr = 0.0;         *mean_ih_ptr = 0.0;
         float sum=0.0, sum2=0.0;         float sum=0.0, sum2=0.0;
  
         for (j = 0; j < ny; j++)  
         {  
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   for (j = 0; j < ny; j++)
                   {
                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(jz[j * nx + i]) continue;                  if isnan(jz_smooth[j * nx + i]) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] == 0.0) continue;                 if (bz[j * nx + i] == 0.0) continue;
                 if (jz[j * nx + i] == 0.0) continue;                  if (jz_smooth[j * nx + i] == 0.0) continue;
                 sum  +=     (jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum  +=     (jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
                 sum2 += fabs(jz[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);                  sum2 += fabs(jz_smooth[j * nx + i]*bz[j * nx + i])*(1/cdelt1)*(rsun_obs/rsun_ref);
                 count_mask++;                 count_mask++;
                 }                 }
          }          }
Line 684  int computeHelicity(float *bz, int *dims
Line 705  int computeHelicity(float *bz, int *dims
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 11:  Sum of Absolute Value per polarity  */  /* Example function 12:  Sum of Absolute Value per polarity  */
  
 //  The Sum of the Absolute Value per polarity is defined as the following: //  The Sum of the Absolute Value per polarity is defined as the following:
 //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes. //  fabs(sum(jz gt 0)) + fabs(sum(jz lt 0)) and the units are in Amperes.
Line 692  int computeHelicity(float *bz, int *dims
Line 713  int computeHelicity(float *bz, int *dims
 //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF), //  Jz = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(RSUN_REF/RSUN_OBS)(RSUN_OBS/RSUN_REF),
 //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS) //     = (Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)
  
 int computeSumAbsPerPolarity(float *bz, float *jz, int *dims, float *totaljzptr,  int computeSumAbsPerPolarity(float *bz, float *jz_smooth, int *dims, float *totaljzptr,
                                                          int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)                                                          int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs)
  
 { {
Line 710  int computeSumAbsPerPolarity(float *bz,
Line 731  int computeSumAbsPerPolarity(float *bz,
               {               {
                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                 if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;
                 if isnan(bz[j * nx + i]) continue;                 if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                  if (bz[j * nx + i] >  0) sum1 += ( jz_smooth[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_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
           }           }
  
Line 720  int computeSumAbsPerPolarity(float *bz,
Line 741  int computeSumAbsPerPolarity(float *bz,
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 12:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */  /* Example function 13:  Mean photospheric excess magnetic energy and total photospheric excess magnetic energy density */
 // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV // The units for magnetic energy density in cgs are ergs per cubic centimeter. The formula B^2/8*PI integrated over all space, dV
 // automatically yields erg per cubic centimeter for an input B in Gauss. // automatically yields erg per cubic centimeter for an input B in Gauss.
 // //
Line 763  int computeFreeEnergy(float *bx, float *
Line 784  int computeFreeEnergy(float *bx, float *
 } }
  
 /*===========================================*/ /*===========================================*/
 /* Example function 13:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */  /* Example function 14:  Mean 3D shear angle, area with shear greater than 45, mean horizontal shear angle, area with horizontal shear angle greater than 45 */
  
 int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims, int computeShearAngle(float *bx, float *by, float *bz, float *bpx, float *bpy, float *bpz, int *dims,
                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,                                           float *meanshear_angleptr, float *area_w_shear_gt_45ptr,


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

Karen Tian
Powered by
ViewCVS 0.9.4