(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.30 and 1.33

version 1.30, 2014/06/02 19:46:44 version 1.33, 2015/01/23 01:18:13
Line 563  int computeBzderivative(float *bz, float
Line 563  int computeBzderivative(float *bz, float
 //              float *noiseby, float *noisebz) //              float *noiseby, float *noisebz)
  
 int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared, int computeJz(float *bx_err, float *by_err, float *bx, float *by, int *dims, float *jz, float *jz_err, float *jz_err_squared,
               int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery)                int *mask, int *bitmask, float cdelt1, double rsun_ref, double rsun_obs,float *derx, float *dery, float *err_term1, float *err_term2)
  
  
 { {
Line 584  int computeJz(float *bx_err, float *by_e
Line 584  int computeJz(float *bx_err, float *by_e
         {         {
             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;
               err_term1[j * 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]);
         }         }
     }     }
  
Line 593  int computeJz(float *bx_err, float *by_e
Line 594  int computeJz(float *bx_err, float *by_e
         {         {
             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;
               err_term2[j * nx + i] = (bx_err[(j+1) * nx + i])*(bx_err[(j+1) * nx + i]) + (bx_err[(j-1) * nx + i])*(bx_err[(j-1) * nx + i]);
         }         }
     }     }
  
Line 602  int computeJz(float *bx_err, float *by_e
Line 604  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(by[j * nx + i]) continue;         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;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     i=nx-1;     i=nx-1;
Line 609  int computeJz(float *bx_err, float *by_e
Line 612  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(by[j * nx + i]) continue;         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;
           err_term1[j * nx + i] = ( (3*by_err[j * nx + i])*(3*by_err[j * nx + i]) + (4*by_err[j * nx + (i+1)])*(4*by_err[j * nx + (i+1)]) + (by_err[j * nx + (i+2)])*(by_err[j * nx + (i+2)]) );
     }     }
  
     j=0;     j=0;
Line 616  int computeJz(float *bx_err, float *by_e
Line 620  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(bx[j * nx + i]) continue;         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;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
     }     }
  
     j=ny-1;     j=ny-1;
Line 623  int computeJz(float *bx_err, float *by_e
Line 628  int computeJz(float *bx_err, float *by_e
     {     {
         if isnan(bx[j * nx + i]) continue;         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;
           err_term2[j * nx + i] = ( (3*bx_err[j*nx + i])*(3*bx_err[j*nx + i]) + (4*bx_err[(j+1) * nx + i])*(4*bx_err[(j+1) * nx + i]) + (bx_err[(j+2) * nx + i])*(bx_err[(j+2) * nx + i]) );
   
     }     }
  
     for (i = 1; i <= nx-2; i++)  
       for (i = 0; i <= nx-1; i++)
     {     {
         for (j = 1; j <= ny-2; j++)          for (j = 0; j <= ny-1; j++)
         {         {
             // 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
             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( err_term1[j * nx + i] + err_term2[j * 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)]) ) ;  
             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]);
             count_mask++;             count_mask++;
   
         }         }
     }     }
         return 0;         return 0;
Line 832  int computeHelicity(float *jz_err, float
Line 837  int computeHelicity(float *jz_err, float
 /* Example function 12:  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 per arcsecond.
 //  The units of jz are in G/pix. In this case, we would have the following: //  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), //  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)
Line 869  int computeSumAbsPerPolarity(float *jz_e
Line 874  int computeSumAbsPerPolarity(float *jz_e
         }         }
     }     }
  
         *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are A */      *totaljzptr    = fabs(sum1) + fabs(sum2);  /* Units are Amperes per arcsecond */
     *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));     *totaljz_err_ptr = sqrt(err)*(1/cdelt1)*fabs((0.00010)*(1/MUNAUGHT)*(rsun_ref/rsun_obs));
     //printf("SAVNCPP=%g\n",*totaljzptr);     //printf("SAVNCPP=%g\n",*totaljzptr);
     //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);     //printf("SAVNCPP_err=%g\n",*totaljz_err_ptr);
Line 1160  int computeR(float *bz_err, float *los,
Line 1165  int computeR(float *bz_err, float *los,
         for (j = 0; j < ny1; j++)         for (j = 0; j < ny1; j++)
         {         {
             index = j * nx1 + i;             index = j * nx1 + i;
               if isnan(pmapn[index]) continue;
               if isnan(rim[index]) continue;
             sum += pmapn[index]*abs(rim[index]);             sum += pmapn[index]*abs(rim[index]);
         }         }
     }     }
Line 1169  int computeR(float *bz_err, float *los,
Line 1176  int computeR(float *bz_err, float *los,
     else     else
         *Rparam = log10(sum);         *Rparam = log10(sum);
  
       //printf("R_VALUE=%f\n",*Rparam);
   
     free_fresize(&fresboxcar);     free_fresize(&fresboxcar);
     free_fresize(&fresgauss);     free_fresize(&fresgauss);
  
Line 1176  int computeR(float *bz_err, float *los,
Line 1185  int computeR(float *bz_err, float *los,
  
 } }
  
   /*===========================================*/
   /* Example function 16: Lorentz force as defined in Fisher, 2012 */
   //
   // This calculation is adapted from Xudong's code
   // at /proj/cgem/lorentz/apps/lorentz.c
   
   int computeLorentz(float *bx,  float *by, float *bz, float *fx, float *fy, float *fz, int *dims,
                      float *totfx_ptr, float *totfy_ptr, float *totfz_ptr, float *totbsq_ptr,
                      float *epsx_ptr, float *epsy_ptr, float *epsz_ptr, int *mask, int *bitmask,
                      float cdelt1, double rsun_ref, double rsun_obs)
   
   {
   
       int nx = dims[0];
       int ny = dims[1];
       int nxny = nx*ny;
       int j = 0;
       int index;
       double totfx = 0, totfy = 0, totfz = 0;
       double bsq = 0, totbsq = 0;
       double epsx = 0, epsy = 0, epsz = 0;
       double area = cdelt1*cdelt1*(rsun_ref/rsun_obs)*(rsun_ref/rsun_obs)*100.0*100.0;
       double k_h = -1.0 * area / (4. * PI) / 1.0e20;
       double k_z = area / (8. * PI) / 1.0e20;
   
       if (nx <= 0 || ny <= 0) return 1;
   
       for (int i = 0; i < nxny; i++)
       {
          if ( mask[i] < 70 || bitmask[i] < 30 ) continue;
          if isnan(bx[i]) continue;
          if isnan(by[i]) continue;
          if isnan(bz[i]) continue;
          fx[i]  = bx[i] * bz[i] * k_h;
          fy[i]  = by[i] * bz[i] * k_h;
          fz[i]  = (bx[i] * bx[i] + by[i] * by[i] - bz[i] * bz[i]) * k_z;
          bsq    = bx[i] * bx[i] + by[i] * by[i] + bz[i] * bz[i];
          totfx  += fx[i]; totfy += fy[i]; totfz += fz[i];
          totbsq += bsq;
       }
   
       *totfx_ptr  = totfx;
       *totfy_ptr  = totfy;
       *totfz_ptr  = totfz;
       *totbsq_ptr = totbsq;
       *epsx_ptr   = (totfx / k_h) / totbsq;
       *epsy_ptr   = (totfy / k_h) / totbsq;
       *epsz_ptr   = (totfz / k_z) / totbsq;
   
       //printf("TOTBSQ=%f\n",*totbsq_ptr);
   
       return 0;
   
   }
   
 /*==================KEIJI'S CODE =========================*/ /*==================KEIJI'S CODE =========================*/
  
 // #include <omp.h> // #include <omp.h>


Legend:
Removed from v.1.30  
changed lines
  Added in v.1.33

Karen Tian
Powered by
ViewCVS 0.9.4