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

version 1.3, 2012/10/23 18:42:56 version 1.4, 2012/12/18 23:05:10
Line 86  int computeAbsFlux(float *bz, int *dims,
Line 86  int computeAbsFlux(float *bz, int *dims,
                 for (i = 0; i < nx; i++)                 for (i = 0; i < nx; i++)
                 {                 {
                   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;
                     //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 117  int computeBh(float *bx, float *by, floa
Line 119  int computeBh(float *bx, float *by, floa
           {           {
             for (i = 0; i < nx; i++)             for (i = 0; i < nx; i++)
               {               {
                   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 156  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 184  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 383  int computeBzderivative(float *bz, int *
Line 391  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 401  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 411  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;
           }           }
  
Line 443  int computeBzderivative(float *bz, int *
Line 457  int computeBzderivative(float *bz, int *
             {             {
                // 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 473  int computeBzderivative(float *bz, int *
Line 490  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 506  int computeJz(float *bx, float *by, int
Line 524  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 534  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 544  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 (j = 0; j < ny; j++)  
               {  
               printf("j=%d\n",j);  
               printf("i=%d\n",i);  
               printf("dery[j*nx+i]=%f\n",dery[j*nx+i]);  
               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]);  
               }  
           }  
         */  
   
  
         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;
                  if isnan(derx[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 */                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 */
                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 */                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 */
                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */                jz[j * nx + i] = (derx[j * nx + i]-dery[j * nx + i]);                                                          /* jz is in units of Gauss/pix */
   
                count_mask++;                count_mask++;
             }             }
           }           }
Line 580  int computeJz(float *bx, float *by, int
Line 589  int computeJz(float *bx, float *by, int
         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;
  
 } }
Line 700  int computeSumAbsPerPolarity(float *bz,
Line 709  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 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[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);                 if (bz[j * nx + i] <= 0) sum2 += ( jz[j * nx + i])*(1/cdelt1)*(0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs);
               }               }
Line 740  int computeFreeEnergy(float *bx, float *
Line 750  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 777  int computeShearAngle(float *bx, float *
Line 789  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.4
Karen Tian
 Powered by
ViewCVS 0.9.4