Cylinder.cxx Source File

00001 // $Id: Cylinder.cxx,v 1.2 2000/12/14 23:18:11 burnett Exp $
00002 //
00003 #include "geometry/Cylinder.h"
00004 
00005 inline static double sqr(double x){return x*x;}
00006 #include <float.h>
00007 
00008 
00009 Cylinder::Cylinder( const Point& origin, const Vector& a, double r )
00010 : Surface( origin, a.unit() )
00011 ,_radius(r)
00012 {}
00013 Cylinder::Cylinder(const Point& origin, double r)
00014 :Surface(origin, Vector(0,0,1))
00015 , _radius(r)
00016 {}
00017 
00018 double
00019 Cylinder::how_near( const Point& x ) const
00020 {  //  Distance from the point x to the infinite Cylinder.
00021    //  The distance will be positive if the point is inside the Cylinder,
00022    //  negative if the point is outside.
00023         Vector r = x - origin();
00024         double rho = sqrt( r.mag2() - sqr(r*axis()) );
00025         return radius()>0?  radius()-rho  : rho+radius() ;
00026 }
00027 
00028 
00029 Vector
00030 Cylinder::normal( const Point& x )const
00031 {
00032    Vector r = x-origin();
00033    Vector rhohat = (r - (r*axis())*axis()).unit(); // unit vector in rho direction
00034    return radius()>0 ? rhohat : -rhohat;
00035 }
00036 double Cylinder::distance(const Point& x, const Vector& vhat, int inout) const
00037 {
00038 //   solve: | rperp(s) | = R,
00039 //   where rperp(s) = r(s)-ahat*(r(s)*ahat)
00040 //         r(s) = x-origin() + vhat*s
00041 
00042         Vector ahat = axis(),
00043                r = x - origin(),
00044                rperp = r - (ahat*r)*ahat,
00045                vperp = vhat-(ahat*vhat)*ahat;
00046 
00047         double a = vperp.mag2(),
00048                b = vperp*rperp,
00049                c = rperp.mag2() - sqr(radius()),
00050                disc = b*b-a*c,
00051                s;
00052         // now want appropriate root of a*s^2 + 2*b*s + c = 0
00053 
00054         if( disc<=0 ) return FLT_MAX;  // misses
00055         if( radius()<0 ) inout *=-1;    // reverse enter/leave sense
00056         int leaving = inout==1;
00057         if( b<0 ) // velocity toward axis
00058         {  s = (leaving) ? (-b+sqrt(disc))/a :
00059                           c/(-b+sqrt(disc));
00060         }
00061         else // heading away from axis: can only leave
00062         {  s = (leaving) ? c/(-b-sqrt(disc)) : FLT_MAX;
00063         }
00064         return s ;
00065 }
00066 
00067 void Cylinder::printOn( std::ostream& os ) const
00068 {       os << "Cylinder surface with origin: " << origin() << "\t"
00069            << "radius: " << radius() << "\tand axis " << axis() << "\n";
00070 }
00071