TankGeometry.h

Go to the documentation of this file.

#ifndef _CachedDirectInjectorOG_TankGeometry_h_
#define _CachedDirectInjectorOG_TankGeometry_h_

/* RB :code is a modified version of code taken from ParticleInjectorNEU, to be merged (utl?) */

#include <utl/GeometryUtilities.h>
#include <utl/Particle.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/Line.h>

using namespace utl;
using namespace std;


namespace CachedDirectInjectorOG {

  class TankGeometry {
  public:
    TankGeometry(const double bottom, const double top, const double radius,
                 const CoordinateSystemPtr cs) :
      fBottomPlane(Point(0, 0, bottom, cs), Vector(0,0,1, cs)),
      fTopPlane(Point(0, 0, top, cs), Vector(0,0,1, cs)),
      fRadius2(Sqr(radius)),
      fCS(cs)
    { }

    bool
    IsHit(const Particle& particle, Point& where)
      const
    {
      const Point& particlePoint = particle.GetPosition();
      const Vector& particleDirection = particle.GetDirection();
      const Line particleLine(particlePoint, particleDirection);

      // try top
      where = Intersection(fTopPlane, particleLine);
      if ((where - fTopPlane.GetAnchor()).GetMag2() < fRadius2)
        return true;

      const double ax = particleDirection.GetX(fCS);
      const double ay = particleDirection.GetY(fCS);
      const double px = particlePoint.GetX(fCS);
      const double py = particlePoint.GetY(fCS);
      const double a = Sqr(ax) + Sqr(ay);
      const double b = 2*(ax*px + ay*py);
      const double c = Sqr(px) + Sqr(py) - fRadius2;
      const double disc = Sqr(b) - 4*a*c;
      if (disc >= 0) {
        const double ta = 2*a;
        const double s0 = -b/ta;

        const double delta = sqrt(disc) / ta;

        const double s = s0 - delta;

        where = particleLine.Propagate(s);
        if (Distance(fBottomPlane, where) >= 0 && Distance(fTopPlane, where) <= 0)
          return true;
      }

      // try bottom, should not be neccesary, as long as particles come from above ...
      where = Intersection(fBottomPlane, particleLine);
      if ((where - fBottomPlane.GetAnchor()).GetMag2() < fRadius2)
        return true;

      return false;
    }

  private:
    const Plane fBottomPlane;
    const Plane fTopPlane;
    const double fRadius2;
    const CoordinateSystemPtr fCS;
  };

}


#endif