(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.11 and 1.13

version 1.11, 2013/05/30 23:26:02 version 1.13, 2013/05/31 22:47:15
Line 513  int computeJz(float *bx_err, float *by_e
Line 513  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 0; j <= ny-1; j++)             for (j = 0; j <= ny-1; j++)
               {               {
                  //if isnan(by[j * nx + i]) continue;                   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 522  int computeJz(float *bx_err, float *by_e
Line 522  int computeJz(float *bx_err, float *by_e
           {           {
             for (j = 1; j <= ny-2; j++)             for (j = 1; j <= ny-2; j++)
               {               {
                  //if isnan(bx[j * nx + i]) continue;                   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 563  int computeJz(float *bx_err, float *by_e
Line 563  int computeJz(float *bx_err, float *by_e
             {             {
                // calculate jz at all points                // calculate jz at all points
                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
   
                //printf("jz[j * nx + i]=%f,i=%d,j=%d\n,",jz[j * nx + i],i,j); // tmp.txt  
                //printf("i=%d,j=%d,jz[j * nx + i]=%f\n,",i,j,jz[j * nx + i]); // tmp1.txt  
   
   
                jz_err[j * nx + i]=0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +                jz_err[j * nx + i]=0.5*sqrt( (bx_err[(j+1) * nx + i]*bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i]*bx_err[(j-1) * nx + i]) +
                                             (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;                                             (by_err[j * nx + (i+1)]*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i-1)]*by_err[j * nx + (i-1)]) ) ;
                jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);                jz_err_squared[j * nx + i]=(jz_err[j * nx + i]*jz_err[j * nx + i]);
Line 842  int computeFreeEnergy(float *bx_err, flo
Line 837  int computeFreeEnergy(float *bx_err, flo
  
         *totpotptr=0.0;         *totpotptr=0.0;
         *meanpotptr=0.0;         *meanpotptr=0.0;
         float sum,err=0.0;          float sum,sum1,err=0.0;
  
         for (i = 0; i < nx; i++)         for (i = 0; i < nx; i++)
           {           {
Line 852  int computeFreeEnergy(float *bx_err, flo
Line 847  int computeFreeEnergy(float *bx_err, flo
                  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;
                  sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);                  sum += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) )*(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0);
                    sum1 += ( ((bpx[j * nx + i] - bx[j * nx + i])*(bpx[j * nx + i] - bx[j * nx + i])) + ((bpy[j * nx + i] - by[j * nx + i])*(bpy[j * nx + i] - by[j * nx + i])) );
                  err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);                  err += (4.0*bx[j * nx + i]*bx[j * nx + i]*bx_err[j * nx + i]*bx_err[j * nx + i]) + (4.0*by[j * nx + i]*by[j * nx + i]*by_err[j * nx + i]*by_err[j * nx + i]);
                  count_mask++;                  count_mask++;
               }               }
           }           }
  
         *meanpotptr      = (sum/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */          *meanpotptr      = (sum1/(8.*PI)) / (count_mask);     /* Units are ergs per cubic centimeter */
         *meanpot_err_ptr = (sqrt(err)) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)          *meanpot_err_ptr = (sqrt(err))*fabs(cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0) / (count_mask*8.*PI); // error in the quantity (sum)/(count_mask)
  
         /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */         /* Units of sum are ergs/cm^3, units of factor are cm^2/pix^2; therefore, units of totpotptr are ergs per centimeter */
         *totpotptr       = (sum)/(8.*PI);         *totpotptr       = (sum)/(8.*PI);


Legend:
Removed from v.1.11  
changed lines
  Added in v.1.13

Karen Tian
Powered by
ViewCVS 0.9.4