(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.3 and 1.5

version 1.3, 2012/10/23 18:42:56 version 1.5, 2013/01/14 18:27:45
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;
                   sum += (fabs(bz[j * nx + i]));                   sum += (fabs(bz[j * nx + i]));
                   count_mask++;                   count_mask++;
                 }                 }
Line 113  int computeBh(float *bx, float *by, floa
Line 114  int computeBh(float *bx, float *by, floa
  
     if (nx <= 0 || ny <= 0) return 1;     if (nx <= 0 || ny <= 0) return 1;
  
         for (j = 0; j < ny; j++)  
           {  
             for (i = 0; i < nx; i++)             for (i = 0; i < nx; i++)
               {               {
               for (j = 0; j < ny; j++)
                 {
                   if isnan(bx[j * nx + i]) continue;
                   if isnan(by[j * nx + i]) continue;
                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );                 bh[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] );
                 sum += bh[j * nx + i];                 sum += bh[j * nx + i];
                 count_mask++;                 count_mask++;
Line 152  int computeGamma(float *bx, float *by, f
Line 155  int computeGamma(float *bx, float *by, f
                 if (bh[j * nx + i] > 100)                 if (bh[j * nx + i] > 100)
                   {                   {
                     if ( mask[j * nx + i] < 70 || 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;
                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));                     sum += (atan (fabs( bz[j * nx + i] / bh[j * nx + i] ))* (180./PI));
                     count_mask++;                     count_mask++;
                   }                   }
Line 179  int computeB_total(float *bx, float *by,
Line 183  int computeB_total(float *bx, float *by,
           {           {
             for (j = 0; j < ny; j++)             for (j = 0; j < ny; j++)
               {               {
                   if isnan(bx[j * nx + i]) continue;
                   if isnan(by[j * nx + i]) continue;
                   if isnan(bz[j * nx + i]) continue;
                 bt[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]);                 bt[j * nx + i] = sqrt( bx[j * nx + i]*bx[j * nx + i] + by[j * nx + i]*by[j * nx + i] + bz[j * nx + i]*bz[j * nx + i]);
               }               }
           }           }
Line 245  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 335  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 383  int computeBzderivative(float *bz, int *
Line 364  int computeBzderivative(float *bz, int *
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                   if isnan(bz[j * nx + i]) continue;
                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;                 derx_bz[j * nx + i] = (bz[j * nx + i+1] - bz[j * nx + i-1])*0.5;
               }               }
           }           }
Line 392  int computeBzderivative(float *bz, int *
Line 374  int computeBzderivative(float *bz, int *
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                   if isnan(bz[j * nx + i]) continue;
                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;                 dery_bz[j * nx + i] = (bz[(j+1) * nx + i] - bz[(j-1) * nx + i])*0.5;
               }               }
           }           }
Line 401  int computeBzderivative(float *bz, int *
Line 384  int computeBzderivative(float *bz, int *
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
                if isnan(bz[j * nx + i]) continue;
              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;              derx_bz[j * nx + i] = ( (-3*bz[j * nx + i]) + (4*bz[j * nx + (i+1)]) - (bz[j * nx + (i+2)]) )*0.5;
           }           }
  
         i=nx-1;         i=nx-1;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
                if isnan(bz[j * nx + i]) continue;
              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;              derx_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[j * nx + (i-1)]) - (-bz[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
                if isnan(bz[j * nx + i]) continue;
              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (-3*bz[j*nx + i]) + (4*bz[(j+1) * nx + i]) - (bz[(j+2) * nx + i]) )*0.5;
           }           }
  
         j=ny-1;         j=ny-1;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
                if isnan(bz[j * nx + i]) continue;
              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;              dery_bz[j * nx + i] = ( (3*bz[j * nx + i]) + (-4*bz[(j-1) * nx + i]) - (-bz[(j-2) * nx + i]) )*0.5;
           }           }
  
  
         /*Just some print statements  
         for (i = 0; i < nx; i++)  
           {  
              for (j = 0; j < ny; j++)  
               {  
               printf("j=%d\n",j);  
               printf("i=%d\n",i);  
               printf("dery_bz[j*nx+i]=%f\n",dery_bz[j*nx+i]);  
               printf("derx_bz[j*nx+i]=%f\n",derx_bz[j*nx+i]);  
               printf("bz[j*nx+i]=%f\n",bz[j*nx+i]);  
               }  
           }  
         */  
   
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;                // if ( (derx_bz[j * nx + i]-dery_bz[j * nx + i]) == 0) continue;
                if ( mask[j * nx + i] < 70 || 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(derx_bz[j * nx + i]) continue;
                  if isnan(dery_bz[j * nx + i]) continue;
                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */                sum += sqrt( derx_bz[j * nx + i]*derx_bz[j * nx + i]  + dery_bz[j * nx + i]*dery_bz[j * nx + i]  ); /* Units of Gauss */
                count_mask++;                count_mask++;
             }             }
Line 453  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 473  int computeBzderivative(float *bz, int *
Line 448  int computeBzderivative(float *bz, int *
 //  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or //  (Gauss/pix)(pix/arcsec)(arcsec/meter)(Newton/Gauss*Ampere*meter)(Ampere^2/Newton)(milliAmpere/Ampere), or
 //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere), //  (Gauss/pix)(1/CDELT1)(RSUN_OBS/RSUN_REF)(1 T / 10^4 Gauss)(1 / 4*PI*10^-7)( 10^3 milliAmpere/Ampere),
 //  where a Tesla is represented as a Newton/Ampere*meter. //  where a Tesla is represented as a Newton/Ampere*meter.
   //
 //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS). //  As an order of magnitude estimate, we can assign 0.5 to CDELT1 and 722500m/arcsec to (RSUN_REF/RSUN_OBS).
 //  In that case, we would have the following: //  In that case, we would have the following:
 //  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or //  (Gauss/pix)(1/0.5)(1/722500)(10^-4)(4*PI*10^7)(10^3), or
Line 485  int computeBzderivative(float *bz, int *
Line 461  int computeBzderivative(float *bz, int *
 //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.) //  =(Gauss/pix)(1/CDELT1)(0.00010)(1/MUNAUGHT)(RSUN_REF/RSUN_OBS)(1000.)
  
 int computeJz(float *bx, float *by, int *dims, float *jz, int computeJz(float *bx, float *by, int *dims, float *jz,
                           float *mean_jz_ptr, float *us_i_ptr, int *mask, 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 496  int computeJz(float *bx, float *by, int
Line 472  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 506  int computeJz(float *bx, float *by, int
Line 480  int computeJz(float *bx, float *by, int
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                    if isnan(by[j * nx + i]) continue;
                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;                  derx[j * nx + i] = (by[j * nx + i+1] - by[j * nx + i-1])*0.5;
               }               }
           }           }
Line 515  int computeJz(float *bx, float *by, int
Line 490  int computeJz(float *bx, float *by, int
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                    if isnan(bx[j * nx + i]) continue;
                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;                  dery[j * nx + i] = (bx[(j+1) * nx + i] - bx[(j-1) * nx + i])*0.5;
               }               }
           }           }
Line 524  int computeJz(float *bx, float *by, int
Line 500  int computeJz(float *bx, float *by, int
         i=0;         i=0;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
                if isnan(by[j * nx + i]) continue;
              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;              derx[j * nx + i] = ( (-3*by[j * nx + i]) + (4*by[j * nx + (i+1)]) - (by[j * nx + (i+2)]) )*0.5;
           }           }
  
         i=nx-1;         i=nx-1;
         for (j = 0; j <= ny-1; j++)         for (j = 0; j <= ny-1; j++)
           {           {
                if isnan(by[j * nx + i]) continue;
              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;              derx[j * nx + i] = ( (3*by[j * nx + i]) + (-4*by[j * nx + (i-1)]) - (-by[j * nx + (i-2)]) )*0.5;
           }           }
  
         j=0;         j=0;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
                if isnan(bx[j * nx + i]) continue;
              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (-3*bx[j*nx + i]) + (4*bx[(j+1) * nx + i]) - (bx[(j+2) * nx + i]) )*0.5;
           }           }
  
         j=ny-1;         j=ny-1;
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
                if isnan(bx[j * nx + i]) continue;
              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;              dery[j * nx + i] = ( (3*bx[j * nx + i]) + (-4*bx[(j-1) * nx + i]) - (-bx[(j-2) * nx + i]) )*0.5;
           }           }
  
         /* Just some print statements  
         for (i = 0; i < nx; i++)          for (i = 0; i <= nx-1; i++)
           {           {
              for (j = 0; j < ny; j++)              for (j = 0; j <= ny-1; j++)
               {               {
               printf("j=%d\n",j);                 /* calculate jz at all points */
               printf("i=%d\n",i);                 jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                              /* jz is in units of Gauss/pix */
               printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);                 count_mask++;
               printf("derx[j*nx+i]=%f\n",derx[j*nx+i]);  
               printf("bx[j*nx+i]=%f\n",bx[j*nx+i]);  
               printf("by[j*nx+i]=%f\n",by[j*nx+i]);  
               }               }
           }           }
         */  
           return 0;
   }
   
   /*===========================================*/
   
   /* Example function 9:  Compute quantities on smoothed Jz array */
   
   int computeJzsmooth(float *bx, float *by, int *dims, float *jz_smooth,
                             float *mean_jz_ptr, float *us_i_ptr, int *mask, int *bitmask,
                             float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)
   
   {
   
           int nx = dims[0], ny = dims[1];
           int i, j, count_mask=0;
   
           if (nx <= 0 || ny <= 0) return 1;
   
           *mean_jz_ptr = 0.0;
           float curl=0.0, us_i=0.0,test_perimeter=0.0,mean_curl=0.0;
  
  
           /* At this point, use the smoothed Jz array with a Gaussian (FWHM of 4 pix and truncation width of 12 pixels) but keep the original array dimensions*/
         for (i = 0; i <= nx-1; i++)         for (i = 0; i <= nx-1; i++)
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
             {             {
                // if ( (derx[j * nx + i]-dery[j * nx + i]) == 0) continue;  
                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) continue;                if ( mask[j * nx + i] < 70 || bitmask[j * nx + i] < 30 ) 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(derx[j * nx + i]) continue;
                us_i += fabs(derx[j * nx + i]-dery[j * nx + i])*(1/cdelt1)*(rsun_ref/rsun_obs)*(0.00010)*(1/MUNAUGHT);         /* us_i is in units of A  / m^2 */                 if isnan(dery[j * nx + i]) continue;
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                 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++;                count_mask++;
             }             }
           }           }
  
           /* Calculate mean vertical current density (mean_curl) and total unsigned vertical current (us_i) using smoothed Jz array and continue conditions above */
         mean_curl        = (curl/count_mask);         mean_curl        = (curl/count_mask);
         printf("mean_curl=%f\n",mean_curl);         printf("mean_curl=%f\n",mean_curl);
         printf("cdelt1, what is it?=%f\n",cdelt1);         printf("cdelt1, what is it?=%f\n",cdelt1);
         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */         *mean_jz_ptr     = curl/(count_mask);        /* mean_jz gets populated as MEANJZD */
         printf("count_mask=%d\n",count_mask);         printf("count_mask=%d\n",count_mask);
         *us_i_ptr        = (us_i);                   /* us_i gets populated as MEANJZD */          *us_i_ptr        = (us_i);                   /* us_i gets populated as TOTUSJZ */
         return 0;         return 0;
  
 } }
  
   
 /*===========================================*/ /*===========================================*/
 /* Example function 9:  Twist Parameter, alpha */  
  
 // The twist parameter, alpha, is defined as alpha = Jz/Bz and the units are in 1/Mm  /* Example function 10:  Twist Parameter, alpha */
   
   // The twist parameter, alpha, is defined as alpha = Jz/Bz. In this case, the calculation
   // for alpha is calculated in the following way (different from Leka and Barnes' approach):
   
          // (sum of all positive Bz + abs(sum of all negative Bz)) = avg Bz
          // (abs(sum of all Jz at positive Bz) + abs(sum of all Jz at negative Bz)) = avg Jz
          // avg alpha = avg Jz / avg Bz
   
   // The 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 637  int computeAlpha(float *bz, int *dims, f
Line 653  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 650  int computeHelicity(float *bz, int *dims
Line 666  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 675  int computeHelicity(float *bz, int *dims
Line 691  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 683  int computeHelicity(float *bz, int *dims
Line 699  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 700  int computeSumAbsPerPolarity(float *bz,
Line 716  int computeSumAbsPerPolarity(float *bz,
             for (j = 0; j < ny; j++)             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 (bz[j * nx + i] >  0) sum1 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                  if isnan(bz[j * nx + i]) continue;
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                  if (bz[j * nx + i] >  0) sum1 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                   if (bz[j * nx + i] <= 0) sum2 += ( jz_smooth[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
           }           }
  
Line 710  int computeSumAbsPerPolarity(float *bz,
Line 727  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 740  int computeFreeEnergy(float *bx, float *
Line 757  int computeFreeEnergy(float *bx, float *
             for (j = 0; j < ny; j++)             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(bx[j * nx + i]) continue;
                    if isnan(by[j * nx + i]) continue;
                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);                  sum += ((    ((bx[j * nx + i])*(bx[j * nx + i]) + (by[j * nx + i])*(by[j * nx + i]) ) -  ((bpx[j * nx + i])*(bpx[j * nx + i]) + (bpy[j * nx + i])*(bpy[j * nx + i]))  )/8.*PI);
                  count_mask++;                  count_mask++;
               }               }
Line 751  int computeFreeEnergy(float *bx, float *
Line 770  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,
Line 777  int computeShearAngle(float *bx, float *
Line 796  int computeShearAngle(float *bx, float *
                  if isnan(bpy[j * nx + i]) continue;                  if isnan(bpy[j * nx + i]) continue;
                  if isnan(bpz[j * nx + i]) continue;                  if isnan(bpz[j * nx + i]) continue;
                  if isnan(bz[j * nx + i]) continue;                  if isnan(bz[j * nx + i]) continue;
                    if isnan(bx[j * nx + i]) continue;
                    if isnan(by[j * nx + i]) continue;
                  /* For mean 3D shear angle, area with shear greater than 45*/                  /* For mean 3D shear angle, area with shear greater than 45*/
                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);                  dotproduct            = (bpx[j * nx + i])*(bx[j * nx + i]) + (bpy[j * nx + i])*(by[j * nx + i]) + (bpz[j * nx + i])*(bz[j * nx + i]);
                  magnitude_potential   = sqrt((bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));                  magnitude_potential   = sqrt((bpx[j * nx + i]*bpx[j * nx + i]) + (bpy[j * nx + i]*bpy[j * nx + i]) + (bpz[j * nx + i]*bpz[j * nx + i]));

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