grid.h

pic 模拟程序！面向对象

      return	Vector3(v.e1()*dt, v.e2()*dt, v.e3()*dt);
      
    case ZRGEOM:
    default:
      {
        Vector3	xf(x0.e1() + v.e1()*dt, x0.e2() + v.e2()*dt, v.e3()*dt);
        //	get new radius:
        Scalar r2 = sqrt(xf.e2()*xf.e2() + xf.e3()*xf.e3());
        
        //	alpha is the angle of rotation, phi2 - phi1, for this timestep
        Scalar	sinAlpha;
        Scalar	cosAlpha;
        if (r2 > 0)
          {
            sinAlpha = xf.e3()/r2;
            cosAlpha = xf.e2()/r2;
          }
        else
          {
            sinAlpha = 0;
            cosAlpha = 1;
          }
        
        //	rotate u into new coordinate system:
        Scalar	u2r = cosAlpha*u.e2() + sinAlpha*u.e3();
        Scalar	u2phi = -sinAlpha*u.e2() + cosAlpha*u.e3();
        u.set_e2(u2r);
        u.set_e3(u2phi);
        
        //	compute and return dx:
        return Vector3(xf.e1() - x0.e1(), r2 - x0.e2(), xf.e3());
      }
    }
}

//--------------------------------------------------------------------
//	Return MKS coordinates from grid coordinates.

inline Vector2 Grid::getMKS(const Vector2& x)
{
  int j = (int)x.e1();
  int k = (int)x.e2();
  Scalar w1 = x.e1() - j;
  Scalar w2 = x.e2() - k;
  if (j == J)
    {
      j--;
      w1 += 1;
    }
  if (k == K)
    {
      k--;
      w2 += 1;
    }	
  return (1-w1)*(1-w2)*X[j][k] + (1-w1)*w2*X[j][k+1] + w1*(1-w2)*X[j+1][k]
    + w1*w2*X[j+1][k+1];
}


//--------------------------------------------------------------------
//	Translate from the location x in grid coordinates to the nearest cell
//	side intersection along the trajectory given by dxMKS in MKS units.
//	Final position is stored in grid coordinates in x.  Assumes orthogonal
//	coordinate system.

#define EDGE_DELTA 1e-6

Boundary*	Grid::translate(Vector2& x, Vector3& dxMKS)
{
#ifdef DEBUG
  int BRANCH;
#endif
  int	j = (int)x.e1();
  int	k = (int)x.e2();
  //	if (x.e1()==(Scalar)j && dxMKS.e1() < 0) j--; // which cell when on cell edge
  //	if (x.e2()==(Scalar)k && dxMKS.e2() < 0) k--; // which cell when on cell edge
  if (x.e1()==(Scalar)j && dxMKS.e1()<0. && j!=0) j--; // which cell when on cell edge
  if (x.e2()==(Scalar)k && dxMKS.e2()<0. && k!=0) k--; // which cell when on cell edge
  int j1 = j+1, k1 = k+1;
  Scalar cell1 = dl1(j,k);
  Scalar dx1 = dxMKS.e1()/cell1;// /cellSize.e1();
  Scalar cell2 = dl2(j,k);
  Scalar dx2 = dxMKS.e2()/cell2;// /cellSize.e2();
  volatile Scalar tmp1 = x.e1() + dx1; // volatile ensures 32-bit truncation
  volatile Scalar tmp2 = x.e2() + dx2;
  Scalar x1New = tmp1;
  Scalar x2New = tmp2;
  double m=0.; // need double precision to ensure proper slope for small moves
  BOOL	mInfinite = FALSE;
  if (x1New != x.e1()){ // (fabs(dx1) > 0) {
    m = ((double)dx2)/dx1;	//	optimization: move m inside tests
    if (fabs(m) > 1e5) {
      mInfinite = TRUE; // slope too large -> loss of precision
      dx1 = 0;
      dxMKS.set_e1(0);
    }
    if (fabs(m) < 1e-5) { // slope too small -> loss of precision
      m = dx2 = 0;
      dxMKS.set_e2(0);
    }
  }
  else mInfinite = TRUE;
  Scalar x2Intersect;
  Scalar x1Intersect;
  
  if (x1New >= j1)
    {
      if (mInfinite) x2Intersect = x2New; // handles vertical move
      else x2Intersect = m*(j1 - x.e1()) + x.e2();
      if (x2Intersect >= k1)					//	top edge
        {
          if (mInfinite) x1Intersect = j1; // handles vertical edge move
          else if (m > 0) x1Intersect = MIN(j1,(k1 - x.e2())/m + x.e1());
          else x1Intersect = MIN(j1,x1New); // handles horizontal move
          if (fabs(x1New - j1) <= EDGE_DELTA && fabs(x2New - k1) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k1 - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k1);
#ifdef DEBUG
          BRANCH = 1;
#endif
          if (bc2[j1][k] && x1Intersect >= j1) return bc2[j1][k]; 
          else return bc1[j][k1];
        }
      else if (x2Intersect <= k)			//	bottom edge
        {
          if (mInfinite) x1Intersect = j1; // handles vertical edge move
          else if (m < 0) x1Intersect = MIN(j1, (k - x.e2())/m + x.e1());
          else x1Intersect = MIN(j1,x1New); // handles horizontal move
          if (fabs(x1New - j1) <= EDGE_DELTA && fabs(x2New - k) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k);
#ifdef DEBUG
          BRANCH = 2;
#endif
          if (bc2[j1][k] && x1Intersect >= j1) return bc2[j1][k];
          else return bc1[j][k];
        }
      else										//	right edge
        {
          if (fabs(x1New - j1) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((j1 - x.e1())*cell1, (x2Intersect - x.e2())*cell2, 0);
          x.set_e1(j1);
          x.set_e2(x2Intersect);
#ifdef DEBUG
          BRANCH = 3;
#endif
          return bc2[j1][k];
        }
    }
  else if (x1New <= j)
    {
      if (mInfinite) x2Intersect = x2New; // handles vertical move
      else x2Intersect = m*(j - x.e1()) + x.e2();
      if (x2Intersect >= k1)					//	top edge
        {
          if (mInfinite) x1Intersect = j; // handles vertical edge move
          else if (m < 0) x1Intersect = MAX(j, (k1 - x.e2())/m + x.e1());
          else x1Intersect = MAX(j, x1New); // handles horizontal move
          if (fabs(x1New - j) <= EDGE_DELTA && fabs(x2New - k1) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k1 - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k1);
#ifdef DEBUG
          BRANCH = 4;
#endif
          if (bc2[j][k] && x1Intersect <= j) return bc2[j][k];
          else return bc1[j][k1];
        }
      else if (x2Intersect <= k)			//	bottom edge
        {
          if (mInfinite) x1Intersect = j; // handles vertical edge move
          if (m > 0) x1Intersect = MAX(j, (k - x.e2())/m + x.e1());
          else x1Intersect = MAX(j,x1New); // handles horizontal move
          if (fabs(x1New - j) <= EDGE_DELTA && fabs(x2New - k) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k);
#ifdef DEBUG
          BRANCH = 5;
#endif
          if (bc2[j][k] && x1Intersect <= j) return bc2[j][k];
          else return bc1[j][k];
        }
      else										//	left edge
        {
          if (fabs(x1New - j) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((j - x.e1())*cell1, (x2Intersect - x.e2())*cell2, 0);
          x.set_e1(j);
          x.set_e2(x2Intersect);
#ifdef DEBUG
          BRANCH = 6;
#endif
          return bc2[j][k];
        }
    }
  else
    {
      if (x2New >= k1)						//	up
        {
          if (mInfinite) x1Intersect = x.e1();	//	vertical
          else if (m > 0) x1Intersect = MIN(j1, MAX(j, (k1 - x.e2())/m + x.e1()));
          else if (m < 0) x1Intersect = MAX(x1New, (k1 - x.e2())/m + x.e1());
          else x1Intersect = x1New; // handles horizontal move (j < x1New < j1)
          if (fabs(x2New - k1) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k1 - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k1);
#ifdef DEBUG
          BRANCH = 7;
#endif
          return bc1[j][k1];
        }
      else if (x2New <= k)					//	down
        {
          if (mInfinite) x1Intersect = x.e1();	//	vertical
          else if (m < 0) x1Intersect = MIN(x1New, (k - x.e2())/m + x.e1()); // right
          else if (m < 0) x1Intersect = MIN(j1, MAX(j, (k - x.e2())/m + x.e1())); // left
          else x1Intersect = x1New; // handles horizontal move (j < x1New < j1)
          if (fabs(x2New - k) <= EDGE_DELTA) {
            dxMKS.set_e1(0);
            dxMKS.set_e2(0);
          }
          else dxMKS -= Vector3((x1Intersect - x.e1())*cell1, (k - x.e2())*cell2, 0);
          x.set_e1(x1Intersect);
          x.set_e2(k);
#ifdef DEBUG
          BRANCH = 8;
#endif
          return bc1[j][k];
        }
      else										//	same cell
        {
          dxMKS.set_e1(0);
          dxMKS.set_e2(0);
          x.set_e1(x1New);
          x.set_e2(x2New);
#ifdef DEBUG
          BRANCH = 9;
#endif
          return NULL;
        }
    }
}


#endif											//	ifndef __GRID_H