#include "munga.h"
#pragma hdrstop

#include "boxsolid.h"
#include "line.h"
#include "plane.h"

//#############################################################################
//##########################    BoxedXAxisCylinder    #########################
//#############################################################################

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedXAxisCylinder::BoxedXAxisCylinder(
	const ExtentBox &extents,
	BoxedSolid::Material material,
	Simulation *owner,
	BoxedSolid *next_solid
):
	BoxedSolid(extents, XAxisCylinderType, material, owner, next_solid)
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedXAxisCylinder::~BoxedXAxisCylinder()
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedXAxisCylinder::IntersectsBounded(const ExtentBox &extents)
{
	Check(this);
	Check(&extents);

	Verify(minX <= extents.minX);
	Verify(maxX >= extents.maxX);
	Verify(minY <= extents.minY);
	Verify(maxY >= extents.maxY);
	Verify(minZ <= extents.minZ);
	Verify(maxZ >= extents.maxZ);

	//
	//---------------------------------------------------------------------
	// Find the center point in the YZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the Y axis.  The X value will
	// automatically be within the cylinder if Y & Z are
	//---------------------------------------------------------------------
	//
	Scalar y = (minY + maxY) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxY - y;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	Scalar y2 = y;
	Scalar z2 = z;
	Clamp(y2, extents.minY, extents.maxY);
	Clamp(z2, extents.minZ, extents.maxZ);
	y2 -= y;
	z2 -= z;
	return y2*y2 + z2*z2 <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedXAxisCylinder::ContainsBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);
	Verify(maxX >= point.x);
	Verify(minY <= point.y);
	Verify(maxY >= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the YZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the Y axis.  The X value will
	// automatically be within the cylinder if Y & Z are
	//---------------------------------------------------------------------
	//
	Scalar y = (minY + maxY) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxY - y;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	y -= point.y;
	z -= point.z;
	return y*y + z*z <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Scalar
	BoxedXAxisCylinder::FindDistanceBelowBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);
	Verify(maxX >= point.x);
	Verify(minY <= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxZ - z;

	//
	//--------------------------------------------------------------------------
	// figure out the thickness of cylinder slice where a plane perpendicular to
	// X drops through the cylinder and the test point
	//--------------------------------------------------------------------------
	//
	z = point.z - z;
	Scalar y = radius - Sqrt(radius*radius - z*z);
	if (point.y > maxY - y)
	{
		return point.y - maxY + y;
	}
	else if (point.y < minY + y)
	{
		return -1.0f;
	}
	else
	{
		return 0.0f;
	}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedXAxisCylinder::HitByBounded(
		Line *line,
		Scalar enters,
		Scalar leaves
	)
{
	Check(this);
	Check(line);

	Verify(enters <= leaves);
	Verify(leaves >= 0.0f);

	Scalar y = (minY + maxY) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxY - y;

	Scalar a =
		line->direction.y*line->direction.y
		 + line->direction.z*line->direction.z;

	//
	//--------------------------------------------------------------------------
	// If the line is parallel to the cylinder, see if it hits the bottom/top of
	// the cylinder in the X-Z plane.  If it does, it hits the cylinder when it
	// hits the box, otherwise it misses altogether
	//--------------------------------------------------------------------------
	//
	if (Small_Enough(a))
	{
		y -= line->origin.y;
		z -= line->origin.z;
		if (y*y + z*z > radius*radius)
		{
			return False;
		}
		line->length = Max(enters, 0.0f);
		return True;
	}

	//
	//-----------------------------------------------------------------------
	// The line is not parallel to the cylinder, so solve the equation giving
	// the intersection of the line with the cylinder using the quadratic
	// equation
	//-----------------------------------------------------------------------
	//
	y = line->origin.y - y;
	z = line->origin.z - z;

	Scalar b = 2.0f * (line->direction.y*y + line->direction.z*z);
	Scalar c = y*y + z*z - radius*radius;
	Scalar i = b*b - 4.0f*a*c;

	if (i < SMALL)
	{
		return False;
	}

	//
	//-------------------------------------------------------------------------
	// If the line does hit the cylinder, update the enter/leaving distances to
	// take the cylinder surface into account
	//-------------------------------------------------------------------------
	//
	Verify(a > SMALL);
	Scalar ratio = Sqrt(a)/(a*a);
	i = Sqrt(i);
	b *= -ratio;
	i *= ratio;
	Scalar enter = (b - i) * 0.5f;
	if (enter > enters)
	{
		enters = enter;
	}

	Scalar leave = enter + i;
	if (leave < leaves)
	{
		leaves = leave;
	}

	//
	//------------------------------------------------------------------------
	// If we have pushed the entering distance after the leaving distance, the
	// cylinder was missed, otherwise it is hit at the new entering distance
	//------------------------------------------------------------------------
	//
	if (enters > leaves || enters > line->length || leaves < 0.0f)
	{
		return False;
	}

	line->length = Max(enters, 0.0f);
	return True;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedXAxisCylinder::TestInstance() const
{
	return solidType == XAxisCylinderType;
}

//#############################################################################
//##########################    BoxedYAxisCylinder    #########################
//#############################################################################

enum {
	Min_X_Bit = 0x01,
	Max_X_Bit = 0x02,
	X_Bits = Min_X_Bit|Max_X_Bit,
	Min_Y_Bit = 0x04,
	Max_Y_Bit = 0x08,
	Y_Bits = Min_Y_Bit|Max_Y_Bit,
	Min_Z_Bit = 0x10,
	Max_Z_Bit = 0x20,
	Z_Bits = Min_Z_Bit|Max_Z_Bit,
	X_Axis_Bit = 0x04,
	Y_Axis_Bit = 0x02,
	Z_Axis_Bit = 0x01
};

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedYAxisCylinder::BoxedYAxisCylinder(
	const ExtentBox &extents,
	BoxedSolid::Material material,
	Simulation *owner,
	BoxedSolid *next_solid
):
	BoxedSolid(extents, YAxisCylinderType, material, owner, next_solid)
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedYAxisCylinder::~BoxedYAxisCylinder()
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::VerifyCollision(BoxedSolidCollision &collision)
{
	Check(this);
	Check(&collision);

	Verify(minX <= collision.collisionSlice.minX);
	Verify(maxX >= collision.collisionSlice.maxX);
	Verify(minY <= collision.collisionSlice.minY);
	Verify(maxY >= collision.collisionSlice.maxY);
	Verify(minZ <= collision.collisionSlice.minZ);
	Verify(maxZ >= collision.collisionSlice.maxZ);

	Scalar
		x,
		z,
		their_x,
		their_z,
		radius,
		their_radius;

	BoxedSolid* solid = collision.GetTreeVolume();

	//
	//--------------------------------------------------------------------------
	// See which type of collision we are going to have to verify, and branch to
	// it
	//--------------------------------------------------------------------------
	//
	switch (solid->solidType)
	{

	//
	//-------------------------------------------------------------------------
	// When colliding with a ramp or horizontal cylinder, go ahead and have the
	// collision slice assume it is from a block.  The error thus generated is
	// actually very small and not really noticed by anyone
	//-------------------------------------------------------------------------
	//
	case BlockType:
	case SphereType:
	case ReducibleBlockType:
	case RampFacingNegativeZType:
	case RampFacingNegativeXType:
	case RampFacingPositiveZType:
	case RampFacingPositiveXType:
	case InvertedRampFacingNegativeZType:
	case InvertedRampFacingNegativeXType:
	case InvertedRampFacingPositiveZType:
	case InvertedRampFacingPositiveXType:
	case WedgeFacingNegativeZAndPositiveXType:
	case WedgeFacingNegativeZAndNegativeXType:
	case WedgeFacingPositiveZAndNegativeXType:
	case WedgeFacingPositiveZAndPositiveXType:
	case XAxisCylinderType:
	case ZAxisCylinderType:
	case RightHandedTileType:
	case LeftHandedTileType:
		return IntersectsBounded(collision.collisionSlice);

	//
	//--------------------------------------------------------------------------
	// When colliding with another upright cylinder, simply calculate the
	// distance between the two center points and see if it is less than the sum
	// of the two radii
	//--------------------------------------------------------------------------
	//
	case YAxisCylinderType:
		x = (maxX + minX) * 0.5f;
		z = (maxZ + minZ) * 0.5f;
		radius = maxX - x;

		their_x = (solid->maxX + solid->minX) * 0.5f;
		their_z = (solid->maxZ + solid->minZ) * 0.5f;
		their_radius = solid->maxX - their_x;

		radius += their_radius;
		x -= their_x;
		z -= their_z;

		return radius*radius >= x*x + z*z;

	//
	//------------------------------------------------------------------------
	// When colliding with a cone, we will base all the calculations on the
	// assumption that the collision can be detected identically by a cylinder
	// created by the the intersection of the bottom plane of this volume with
	// the cone
	//------------------------------------------------------------------------
	//
	case ConeType:
		x = (maxX + minX) * 0.5f;
		z = (maxZ + minZ) * 0.5f;
		radius = maxX - x;

		their_x = (solid->maxX + solid->minX) * 0.5f;
		their_z = (solid->maxZ + solid->minZ) * 0.5f;
		their_radius = solid->maxX - their_x;
		their_radius *= (solid->maxY - minY) / (solid->maxY - solid->minY);

		radius += their_radius;
		x -= their_x;
		z -= their_z;

		return radius*radius >= x*x + z*z;

	//
	//------------------------------
	// Fail on the unsupported types
	//------------------------------
	//
	default:
		#if defined(LAB_ONLY)
			DEBUG_STREAM << collision.GetTreeVolume()->solidType << endl << std::flush;
			Fail("Unsupported collision primative\n");
		#endif
		return False;
	}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::ProcessCollision(
		BoxedSolidCollision &collision,
		const Vector3D &velocity,
		BoxedSolidCollisionList *last_collisions,
		Normal *normal,
		Scalar *penetration
	)
{
	Check(this);
	Check(&collision);
	Check(&velocity);
	Check_Pointer(normal);
	Check_Pointer(penetration);

	Scalar
		result,
		len;
	int
		i,
		axes,
		axis,
		face,
		mask = 0,
		opp_face,
		temp;
	bool ok = false;
	BoxedSolid *solid = collision.GetTreeVolume();
	Vector3D
		position;

	//
	//-----------------------------------------------------
	// First verify that the collision really should happen
	//-----------------------------------------------------
	//
	if (Small_Enough(velocity.LengthSquared()))
	{
		return False;
	}
	if (!VerifyCollision(collision))
	{
		return False;
	}

	//
	//------------------------------------------------------------------------
	// Look at each of the three coordinates of motion, and make sure that the
	// collision slice is valid for at least one of them
	//------------------------------------------------------------------------
	//
	for (axis=X_Axis; axis<=Z_Axis; ++axis)
	{
		//
		//-------------------------------------------------------------------
		// Make sure we have actual velocity along this axis, then figure out
		// which face to test
		//-------------------------------------------------------------------
		//
		(*normal)[axis] = 0.0f;
		mask <<= 1;
		if (Small_Enough(velocity[axis]))
		{
			continue;
		}
		face = (axis<<1) + (velocity[axis]>0.0f);
		opp_face = face ^ 1;

		//
		//----------------------------------------------------------------------
		// Find out which faces of the cylinder's bounding box should be
		// considered for purposes of the collision.  Sides are only considered
		// valid if the leading face of the collision slice matches the leading
		// face of the disk.
		//----------------------------------------------------------------------
		//
		if (
			collision.collisionSlice[face] == (*this)[face] &&
			(
				(face&1) && (*this)[opp_face] < collision.collisionSlice[opp_face]
				 ||
					!(face&1)
					 && (*this)[opp_face] > collision.collisionSlice[opp_face]
			)
		)
			mask |= 1;
	}

	//
	//-----------------------------------------------------------------
	// Handle the actual collision based upon the type of model and the
	// components involved in the collision
	//-----------------------------------------------------------------
	//
	switch (solid->solidType)
	{

	//
	//--------------------------
	// Handle the Z facing Ramps
	//--------------------------
	//
	case RampFacingNegativeZType:
		normal->y = solid->maxZ - solid->minZ;
		normal->z = solid->maxY - solid->minY;
		position.y = collision.collisionSlice.minY - solid->maxY;
		position.z = collision.collisionSlice.minZ - solid->minZ;

		//
		//--------------------------------------------------------------------
		// The normal y and z values are assumed at this point to point in the
		// right direction although their lengths are screwed up
		//--------------------------------------------------------------------
		//
Compute_X:
		len = Sqrt(normal->y*normal->y + normal->z*normal->z);
		Verify(!Small_Enough(len));
		normal->y /= len;
		normal->z /= len;
		Check_Fpu();

		result = velocity.y*normal->y + velocity.z*normal->z;
		if (result > -SMALL)
		{
			if (mask&X_Axis_Bit)
			{
				goto X_Only;
			}
			else
			{
				break;
			}
		}
		*penetration =
			fabs((position.y*normal->y + position.z*normal->z) / result);
		Check_Fpu();

		//
		//-----------------------------------------------------------------------
		// If the side of the ramp was involved in the collision, pick either the
		// ramp or side to reflect off of based upon which opposes motion more
		//-----------------------------------------------------------------------
		//
		if (mask&X_Axis_Bit)
		{
			//
			//----------------------------------------------------------------
			// Compute plane equation and find the distance from the most
			// penetrating slice point to the plane.  If the ray comes out the
			// ramp, we hit that, otherwise we hit the side
			//----------------------------------------------------------------
			//
			result = fabs((solid->maxX - solid->minX) / velocity.x);
			Check_Fpu();
			if (result < *penetration)
			{
X_Only:
				normal->x = (velocity.x > 0.0f) ? -1.0f : 1.0f;
				normal->y = 0.0f;
				normal->z = 0.0f;
				*penetration = result;
			}
		}
		goto Reflect;

	case RampFacingPositiveZType:
		normal->y = solid->maxZ - solid->minZ;
		normal->z = solid->minY - solid->maxY;
		position.y = collision.collisionSlice.minY - solid->maxY;
		position.z = collision.collisionSlice.maxZ - solid->maxZ;
		goto Compute_X;

	case InvertedRampFacingNegativeZType:
		normal->y = solid->minZ - solid->maxZ;
		normal->z = solid->maxY - solid->minY;
		position.y = collision.collisionSlice.maxY - solid->maxY;
		position.z = collision.collisionSlice.minZ - solid->maxZ;
		goto Compute_X;

	case InvertedRampFacingPositiveZType:
		normal->y = solid->minZ - solid->maxZ;
		normal->z = solid->minY - solid->maxY;
		position.y = collision.collisionSlice.maxY - solid->maxY;
		position.z = collision.collisionSlice.maxZ - solid->minZ;
		goto Compute_X;

	//
	//--------------------------
	// Handle the X facing ramps
	//--------------------------
	//
	case RampFacingNegativeXType:
		normal->x = solid->maxY - solid->minY;
		normal->y = solid->maxX - solid->minX;
		position.y = collision.collisionSlice.minY - solid->maxY;
		position.x = collision.collisionSlice.minX - solid->minX;

		//
		//--------------------------------------------------------------------
		// The normal x and y values are assumed at this point to point in the
		// right direction although their lengths are screwed up
		//--------------------------------------------------------------------
		//
Compute_Z:
		len = Sqrt(normal->x*normal->x + normal->y*normal->y);
		Verify(!Small_Enough(len));
		normal->x /= len;
		normal->y /= len;
		Check_Fpu();

		result = velocity.y*normal->y + velocity.x*normal->x;
		if (result > -SMALL)
		{
			if (mask&Z_Axis_Bit)
			{
				goto Z_Only;
			}
			else
			{
				break;
			}
		}
		*penetration =
			fabs((position.y*normal->y + position.x*normal->x) / result);
		Check_Fpu();

		//
		//-----------------------------------------------------------------------
		// If the side of the ramp was involved in the collision, pick either the
		// ramp or side to reflect off of based upon which opposes motion more
		//-----------------------------------------------------------------------
		//
		if (mask&Z_Axis_Bit)
		{
			//
			//----------------------------------------------------------------
			// Compute plane equation and find the distance from the most
			// penetrating slice point to the plane.  If the ray comes out the
			// ramp, we hit that, otherwise we hit the side
			//----------------------------------------------------------------
			//
			result = fabs((solid->maxZ - solid->minZ) / velocity.z);
			Check_Fpu();
			if (result < *penetration)
			{
Z_Only:
				normal->x = 0.0f;
				normal->y = 0.0f;
				normal->z = (velocity.z > 0.0f) ? -1.0f : 1.0f;
				*penetration = result;
			}
		}
		goto Reflect;

	case RampFacingPositiveXType:
		normal->x = solid->minY - solid->maxY;
		normal->y = solid->maxX - solid->minX;
		position.y = collision.collisionSlice.minY - solid->maxY;
		position.x = collision.collisionSlice.maxX - solid->maxX;
		goto Compute_Z;

	case InvertedRampFacingNegativeXType:
		normal->x = solid->maxY - solid->minY;
		normal->y = solid->minX - solid->maxX;
		position.y = collision.collisionSlice.maxY - solid->maxY;
		position.x = collision.collisionSlice.minX - solid->maxX;
		goto Compute_Z;

	case InvertedRampFacingPositiveXType:
		normal->x = solid->minY - solid->maxY;
		normal->y = solid->minX - solid->maxX;
		position.y = collision.collisionSlice.maxY - solid->maxY;
		position.x = collision.collisionSlice.maxX - solid->minX;
		goto Compute_Z;

	//
	//------------------
	// Handle the wedges
	//------------------
	//
	case WedgeFacingNegativeZAndPositiveXType:
		normal->x = solid->minZ - solid->maxZ;
		normal->z = solid->maxX - solid->minX;
		position.x = collision.collisionSlice.maxX - solid->minX;
		position.z = collision.collisionSlice.minZ - solid->minZ;

		//
		//--------------------------------------------------------------------
		// The normal x and z values are assumed at this point to point in the
		// right direction although their lengths are screwed up
		//--------------------------------------------------------------------
		//
Compute_Y:
		len = Sqrt(normal->x*normal->x + normal->z*normal->z);
		Verify(!Small_Enough(len));
		normal->x /= len;
		normal->z /= len;
		Check_Fpu();

		result = velocity.x*normal->x + velocity.z*normal->z;
		if (result > -SMALL)
		{
			if (mask&Y_Axis_Bit)
			{
				goto Y_Only;
			}
			else
			{
				break;
			}
		}
		*penetration =
			fabs((position.x*normal->x + position.z*normal->z) / result);
		Check_Fpu();

		//
		//-----------------------------------------------------------------------
		// If the top of the wedge was involved in the collision, pick either the
		// side or top to reflect off of based upon which opposes motion more
		//-----------------------------------------------------------------------
		//
		if (mask&Y_Axis_Bit)
		{
			//
			//----------------------------------------------------------------
			// Compute plane equation and find the distance from the most
			// penetrating slice point to the plane.  If the ray comes out the
			// ramp, we hit that, otherwise we hit the side
			//----------------------------------------------------------------
			//
			result = fabs((solid->maxY - solid->minY) / velocity.y);
			Check_Fpu();
			if (result < *penetration)
			{
Y_Only:
				normal->x = 0.0f;
				normal->y = (velocity.y > 0.0f) ? -1.0f : 1.0f;
				normal->z = 0.0f;
				*penetration = result;
			}
		}
		goto Reflect;

	case WedgeFacingNegativeZAndNegativeXType:
		normal->x = solid->maxZ - solid->minZ;
		normal->z = solid->maxX - solid->minX;
		position.x = collision.collisionSlice.minX - solid->maxX;
		position.z = collision.collisionSlice.minZ - solid->minZ;
		goto Compute_Y;

	case WedgeFacingPositiveZAndNegativeXType:
		normal->x = solid->maxZ - solid->minZ;
		normal->z = solid->minX - solid->maxX;
		position.x = collision.collisionSlice.minX - solid->minX;
		position.z = collision.collisionSlice.maxZ - solid->minZ;
		goto Compute_Y;

	case WedgeFacingPositiveZAndPositiveXType:
		normal->x = solid->minZ - solid->maxZ;
		normal->z = solid->minX - solid->maxX;
		position.x = collision.collisionSlice.maxX - solid->maxX;
		position.z = collision.collisionSlice.maxZ - solid->minZ;
		goto Compute_Y;

	//
	//----------------------------------------------------------------------
	//----------------------------------------------------------------------
	//
	case XAxisCylinderType:
		normal->y = minY + maxY - solid->minY - solid->maxY;
		normal->z = minZ + maxZ - solid->minZ - solid->maxZ;
		len = normal->y*normal->y + normal->z*normal->z;
		if (Small_Enough(len))
		{
			goto Reverse;
		}
		len = Sqrt(len);
		normal->y /= len;
		normal->z /= len;
		Check_Fpu();
		*penetration = 1.0f;
		goto Reflect;

	case ZAxisCylinderType:
		normal->x = minX + maxX - solid->minX - solid->maxX;
		normal->y = minY + maxY - solid->minY - solid->maxY;
		len = normal->x*normal->x + normal->y*normal->y;
		if (Small_Enough(len))
		{
			goto Reverse;
		}
		len = Sqrt(len);
		normal->x /= len;
		normal->y /= len;
		Check_Fpu();
		*penetration = 1.0f;
		goto Reflect;

	case YAxisCylinderType:
		normal->x = minX + maxX - solid->minX - solid->maxX;
		normal->z = minZ + maxZ - solid->minZ - solid->maxZ;
		len = normal->x*normal->x + normal->z*normal->z;
		if (Small_Enough(len))
		{
			goto Reverse;
		}
		len = Sqrt(len);
		normal->x /= len;
		normal->z /= len;
		Check_Fpu();
		*penetration = 1.0f;

		if (mask & Y_Axis_Bit)
		{
			len = (velocity.z > 0.0f) ? -1.0f : 1.0f;
			if (len*velocity.y < velocity.x*normal->x + velocity.z*normal->z)
			{
				normal->x = 0.0f;
				normal->y = len;
				normal->z = 0.0f;
			}
		}
		goto Reflect;

	//
	//-------------
	// Handle cones
	//-------------
	//
	case ConeType:
		normal->x = minX + maxX - solid->minX - solid->maxX;
		normal->z = minZ + maxZ - solid->minZ - solid->maxZ;
		len = normal->x*normal->x + normal->z*normal->z;
		*penetration = 1.0f;
		if (Small_Enough(len))
		{
			goto Reverse;
		}
		normal->y =
			0.5f * (solid->maxX - solid->minX) * Sqrt(len)
			 / (solid->maxY - solid->minY);
		Check_Fpu();
		normal->Normalize(*normal);
		goto Reflect;

	//
	//---------------
	// Handle spheres
	//---------------
	//
	case SphereType:
		normal->x = minX + maxX - solid->minX - solid->maxX;
		normal->y = minY + maxY - solid->minY - solid->maxY;
		normal->z = minZ + maxZ - solid->minZ - solid->maxZ;
		len = normal->x*normal->x + normal->y*normal->y + normal->z*normal->z;
		if (Small_Enough(len))
		{
			goto Reverse;
		}
		len = Sqrt(len);
		normal->x /= len;
		normal->y /= len;
		normal->z /= len;
		Check_Fpu();
		*penetration = 1.0f;
		goto Reflect;

	//
	//---------------------------------
	// Handle reflecting off of a block
	//---------------------------------
	//
	case BlockType:
	case ReducibleBlockType:

		//
		//---------------------------------------------------
		// Count the number of axes involved in the collision
		//---------------------------------------------------
		//
		temp = mask;
		for (axes=0; temp; temp>>=1)
		{
			if (temp&1)
			{
				++axes;
			}
		}

		//
		//------------------------------------------------------
		// Handle as appropriate to the number of sides involved
		//------------------------------------------------------
		//
		switch (axes)
		{

		//
		//---------------------------------------------------------------------
		// If no sides were found to be involved in the collision, then we have
		// a weird case where the disk has completely penetrated through the
		// block
		//---------------------------------------------------------------------
		//
		case 0:

			//
			//-------------------------------------------------------------------
			// This is where we will check the current volume against our list of
			// volumes hit last time.  If we have not already hit this solid,
			// reverse the motion as we have penetrated the object
			//-------------------------------------------------------------------
			//
			i = 0;
			if (last_collisions)
			{
				for (i=0; i<last_collisions->GetCollisionCount(); ++i)
				{
					if ((*last_collisions)[i].GetTreeVolume() == solid)
					{
						break;
					}
				}
			}
			if (!last_collisions || i == last_collisions->GetCollisionCount())
			{
Reverse:
				*penetration = 1.0f;
				result = -1.0f / velocity.Length();
				Check_Fpu();
				normal->x = velocity.x*result;
				normal->y = velocity.y*result;
				normal->z = velocity.z*result;
				ok = true;
				break;
			}

			//
			//--------------------------------------------------------------------
			// Since this is not a fresh hit, check to see if any part of the
			// colliding disk sticks out of the slice.  If so, then we are heading
			// out of the solid in that direction, so let the motion continue.
			// WHAT IS THE VELOCITY TRIGGER FOR????
			//--------------------------------------------------------------------
			//
			position.x = (maxX + minX) * 0.5f;
			position.y = (maxY + minY) * 0.5f;
			position.z = (maxZ + minZ) * 0.5f;
			for (face=0; face<6; ++face)
			{
				i = face >> 1;
				if (face&1)
				{
					if
					(
						(*this)[face] > collision.collisionSlice[face]
						 && position[i] > collision.collisionSlice[face]
						 && velocity[i] < 5.0f
					)
						break;
				}
				else if
				(
					(*this)[face] < collision.collisionSlice[face]
					 && position[i] < collision.collisionSlice[face]
					 && velocity[i] > -5.0f
				)
					break;
			}
			if (face != 6)
			{
				break;
			}

			//
			//-------------------------------------------------------------------
			// Otherwise, generate a normal opposing motion towards the center of
			// the solid.  For each component, if it does not oppose motion, zero
			// it out.  This will give a better approximation of the brick
			//-------------------------------------------------------------------
			//
			normal->x = position.x - (solid->minX + solid->maxX) * 0.5f;
			normal->y = position.y - (solid->minY + solid->maxY) * 0.5f;
			normal->z = position.z - (solid->minZ + solid->maxZ) * 0.5f;
			for (i=0; i<3; ++i)
			{
				if ((*normal)[i]*velocity[i] >= 0.0)
				{
					(*normal)[i] = 0.0;
				}
			}

			//
			//----------------------------------------------------------------
			// Normalize the vector, making sure we don't get killed by a zero
			// case
			//----------------------------------------------------------------
			//
			result = normal->Length();
			if (!Small_Enough(result))
			{
				normal->x /= result;
				normal->y /= result;
				normal->z /= result;
				Check_Fpu();
			}
			else
			{
				break;
			}
			*penetration = 1.0f;
			goto Reflect;

		//
		//-----------------------------------------------------------------
		// If only one side is involved in the collision, reflect off of it
		//-----------------------------------------------------------------
		//
		case 1:
Test_Corner:
			position.x = 0.5f * (minX + maxX);
			position.z = 0.5f * (minZ + maxZ);
			switch (mask)
			{
			case X_Axis_Bit:
				face = X_Axis << 1;

				//
				//----------------------------------
				// Figure out the z bounce direction
				//----------------------------------
				//
				if (collision.collisionSlice.minZ > position.z)
				{
					normal->z = position.z - collision.collisionSlice.minZ;
				}
				else if (collision.collisionSlice.maxZ < position.z)
				{
					normal->z = position.z - collision.collisionSlice.maxZ;
				}

				//
				//----------------------------------
				// Figure out the x bounce direction
				//----------------------------------
				//
				len = maxX - position.x;
				normal->x = len*len - normal->z*normal->z;
				Verify(normal->x >= 0.0f);
				normal->x = Sqrt(normal->x);
				normal->x /= len;
				normal->z /= len;
				Check_Fpu();
				if (velocity.x > 0.0f)
				{
					normal->x = -normal->x;
				}
				break;

			case Y_Axis_Bit:
				face = Y_Axis << 1;
				*penetration =
					fabs(
						(
							collision.collisionSlice[face]
							 - collision.collisionSlice[face^1]
						) / velocity[face>>1]
					);
				Check_Fpu();
				goto Pick_Side;

			case Z_Axis_Bit:
				face = Z_Axis << 1;

				//
				//----------------------------------
				// Figure out the x bounce direction
				//----------------------------------
				//
				if (collision.collisionSlice.minX > position.x)
				{
					normal->x = position.x - collision.collisionSlice.minX;
				}
				else if (collision.collisionSlice.maxX < position.x)
				{
					normal->x = position.x - collision.collisionSlice.maxX;
				}

				//
				//----------------------------------
				// Figure out the z bounce direction
				//----------------------------------
				//
				len = maxX - position.x;
				normal->z = len*len - normal->x*normal->x;
				Verify(normal->z >= 0.0f);
				normal->z = Sqrt(normal->z);
				normal->x /= len;
				normal->z /= len;
				Check_Fpu();
				if (velocity.z > 0.0f)
				{
					normal->z = -normal->z;
				}
				break;
			}
			*penetration =
				fabs(
					(
						collision.collisionSlice[face]
						 - collision.collisionSlice[face^1]
					) / velocity[face>>1]
				);
			Check_Fpu();

			//
			//-------------------------------------------------------------------
			// Check to make sure that the normal opposes velocity.  If it
			// doesn't look to see if this is the first hit on the object, and if
			// so, reverse the velocity
			//-------------------------------------------------------------------
			//
Reflect:
			Check(normal);
			ok = (*normal * velocity < 0.0f);
			if (!ok)
			{
				i = 0;
				if (last_collisions)
				{
					for (i=0; i<last_collisions->GetCollisionCount(); ++i)
					{
						if ((*last_collisions)[i].GetTreeVolume() == solid)
						{
							break;
						}
					}
				}
				if (!last_collisions || i == last_collisions->GetCollisionCount())
				{
					goto Reverse;
				}
			}
			break;

		 //
		 //---------------------------------------------------------------------
		 // If two faces were involved, we hit an edge, so figure out which face
		 // we really hit and mask off the other face
		 //---------------------------------------------------------------------
		 //
		 case 2:
			switch (mask)
			{
			case X_Axis_Bit|Z_Axis_Bit:
				position.x =
					collision.collisionSlice.maxX - collision.collisionSlice.minX;
				position.z =
					collision.collisionSlice.maxZ - collision.collisionSlice.minZ;
				mask ^=
					(fabs(position.z*velocity.x) > fabs(position.x*velocity.z))
						? Z_Axis_Bit : X_Axis_Bit;
				break;

			case X_Axis_Bit|Y_Axis_Bit:
				position.x =
					collision.collisionSlice.maxX - collision.collisionSlice.minX;
				position.y =
					collision.collisionSlice.maxY - collision.collisionSlice.minY;
				mask ^=
					(fabs(position.y*velocity.x) > fabs(position.x*velocity.y))
						? Y_Axis_Bit : X_Axis_Bit;
				break;

			case Y_Axis_Bit|Z_Axis_Bit:
				position.y =
					collision.collisionSlice.maxY - collision.collisionSlice.minY;
				position.z =
					collision.collisionSlice.maxZ - collision.collisionSlice.minZ;
				mask ^=
					(fabs(position.y*velocity.z) > fabs(position.z*velocity.y))
						? Y_Axis_Bit : Z_Axis_Bit;
				break;
			}
			goto Test_Corner;


		//
		//-----------------------------------------------------------------------
		// Otherwise, project the intruding vertex back along the velocity vector
		// to try and find out which face to bounce off of
		//-----------------------------------------------------------------------
		//
		default:
			position.x =
				collision.collisionSlice.maxX - collision.collisionSlice.minX;
			position.y =
				collision.collisionSlice.maxY - collision.collisionSlice.minY;
			position.z =
				collision.collisionSlice.maxZ - collision.collisionSlice.minZ;

			position.x = fabs(position.x/velocity.x);
			position.y = fabs(position.y/velocity.y);
			position.z = fabs(position.z/velocity.z);

			if (position.x < position.y)
			{
				face = (position.x < position.z) ? (X_Axis<<1) : (Z_Axis<<1);
			}
			else
			{
				face = (position.y < position.z) ? (Y_Axis<<1) : (Z_Axis<<1);
			}
			*penetration = position[face>>1];

Pick_Side:
			if (velocity[face>>1] > 0.0f)
			{
				(*normal)[face>>1] = -1.0f;
			}
			else
			{
				(*normal)[face>>1] = 1.0f;
			}
			goto Reflect;
		}
		break;

	#if defined(LAB_ONLY)
		default:
			Fail("Unsupported collision type!...\n");
			break;
	#endif
	}

	if (ok)
	{
		Check(normal);
		if (*penetration < 0.0f)
		{
			return False;
		}
	}
	return ok;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::IntersectsBounded(const ExtentBox &extents)
{
	Check(this);
	Check(&extents);

	Verify(minX <= extents.minX);
	Verify(maxX >= extents.maxX);
	Verify(minY <= extents.minY);
	Verify(maxY >= extents.maxY);
	Verify(minZ <= extents.minZ);
	Verify(maxZ >= extents.maxZ);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxX - x;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	Scalar x2 = x;
	Scalar z2 = z;
	Clamp(x2, extents.minX, extents.maxX);
	Clamp(z2, extents.minZ, extents.maxZ);
	x2 -= x;
	z2 -= z;
	return x2*x2 + z2*z2 <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::ContainsBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);
	Verify(maxX >= point.x);
	Verify(minY <= point.y);
	Verify(maxY >= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxX - x;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	x -= point.x;
	z -= point.z;
	return x*x + z*z <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Scalar
	BoxedYAxisCylinder::FindDistanceBelowBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);	//z);   //GY
	Verify(maxX >= point.x);	//z);   //GY
	Verify(minY <= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxX - x;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	x -= point.x;
	z -= point.z;
	if (x*x + z*z > radius*radius)
	{
		return -1.0f;
	}

	Scalar height = point.y - maxY;
	return Abs(height);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::HitByBounded(
		Line *line,
		Scalar enters,
		Scalar leaves
	)
{
	Check(this);
	Check(line);

	Verify(enters <= leaves);
	Verify(leaves >= 0.0f);

	Scalar x = (minX + maxX) * 0.5f;
	Scalar z = (minZ + maxZ) * 0.5f;
	Scalar radius = maxX - x;

	Scalar a =
		line->direction.x*line->direction.x
		 + line->direction.z*line->direction.z;

	//
	//--------------------------------------------------------------------------
	// If the line is parallel to the cylinder, see if it hits the bottom/top of
	// the cylinder in the X-Z plane.  If it does, it hits the cylinder when it
	// hits the box, otherwise it misses altogether
	//--------------------------------------------------------------------------
	//
	if (Small_Enough(a))
	{
		x -= line->origin.x;
		z -= line->origin.z;
		if (x*x + z*z > radius*radius)
		{
			return False;
		}
		line->length = Max(enters, 0.0f);
		return True;
	}

	//
	//-----------------------------------------------------------------------
	// The line is not parallel to the cylinder, so solve the equation giving
	// the intersection of the line with the cylinder using the quadratic
	// equation
	//-----------------------------------------------------------------------
	//
	x = line->origin.x - x;
	z = line->origin.z - z;

	Scalar b = 2.0f * (line->direction.x*x + line->direction.z*z);
	Scalar c = x*x + z*z - radius*radius;
	Scalar i = b*b - 4.0f*a*c;

	if (i < SMALL)
	{
		return False;
	}

	//
	//-------------------------------------------------------------------------
	// If the line does hit the cylinder, update the enter/leaving distances to
	// take the cylinder surface into account
	//-------------------------------------------------------------------------
	//
	Verify(a > SMALL);
	Scalar ratio = Sqrt(a)/(a*a);
	i = Sqrt(i);
	b *= -ratio;
	i *= ratio;
	Scalar enter = (b - i) * 0.5f;
	if (enter > enters)
	{
		enters = enter;
	}

	Scalar leave = enter + i;
	if (leave < leaves)
	{
		leaves = leave;
	}

	//
	//------------------------------------------------------------------------
	// If we have pushed the entering distance after the leaving distance, the
	// cylinder was missed, otherwise it is hit at the new entering distance
	//------------------------------------------------------------------------
	//
	if (enters > leaves || enters > line->length || leaves < 0.0f)
	{
		return False;
	}

	line->length = Max(enters, 0.0f);
	return True;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedYAxisCylinder::TestInstance() const
{
	return solidType == YAxisCylinderType;
}

//#############################################################################
//##########################    BoxedZAxisCylinder    #########################
//#############################################################################

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedZAxisCylinder::BoxedZAxisCylinder(
	const ExtentBox &extents,
	BoxedSolid::Material material,
	Simulation *owner,
	BoxedSolid *next_solid
):
	BoxedSolid(extents, ZAxisCylinderType, material, owner, next_solid)
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
BoxedZAxisCylinder::~BoxedZAxisCylinder()
{
	Check_Pointer(this);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedZAxisCylinder::IntersectsBounded(const ExtentBox &extents)
{
	Check(this);
	Check(&extents);

	Verify(minX <= extents.minX);
	Verify(maxX >= extents.maxX);
	Verify(minY <= extents.minY);
	Verify(maxY >= extents.maxY);
	Verify(minZ <= extents.minZ);
	Verify(maxZ >= extents.maxZ);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar y = (minY + maxY) * 0.5f;
	Scalar radius = maxX - x;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	Scalar x2 = x;
	Scalar y2 = y;
	Clamp(x2, extents.minX, extents.maxX);
	Clamp(y2, extents.minY, extents.maxY);
	x2 -= x;
	y2 -= y;
	return x2*x2 + y2*y2 <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedZAxisCylinder::ContainsBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);
	Verify(maxX >= point.x);
	Verify(minY <= point.y);
	Verify(maxY >= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar y = (minY + maxY) * 0.5f;
	Scalar radius = maxX - x;

	//
	//-----------------------------------------------------------------------
	// Convert the point to the coordinates of the cylinder, putting the
	// center point of the cylinder at the origin.  Note that we are
	// subtracting the point from the center point as opposed to the normal
	// way.  This will result in both X and Y being negated, but since we are
	// squaring them, this will not matter
	//-----------------------------------------------------------------------
	//
	x -= point.x;
	y -= point.y;
	return x*x + y*y <= radius*radius;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Scalar
	BoxedZAxisCylinder::FindDistanceBelowBounded(const Point3D &point)
{
	Check(this);
	Check(&point);

	Verify(minX <= point.x);//z); //GY
	Verify(maxX >= point.x);//z); //GY
	Verify(minY <= point.y);
	Verify(minZ <= point.z);
	Verify(maxZ >= point.z);

	//
	//---------------------------------------------------------------------
	// Find the center point in the XZ plane of the cylinder, and find the
	// radius of the cylinder by measuring in the X axis.  The Y value will
	// automatically be within the cylinder if X & Z are
	//---------------------------------------------------------------------
	//
	Scalar x = (minX + maxX) * 0.5f;
	Scalar radius = maxX - x;

	//
	//--------------------------------------------------------------------------
	// figure out the thickness of cylinder slice where a plane perpendicular to
	// X drops through the cylinder and the test point
	//--------------------------------------------------------------------------
	//
	x = point.x - x;
	Scalar y = radius - Sqrt(radius*radius - x*x);
	if (point.y > maxY - y)
	{
		return point.y - maxY + y;
	}
	else if (point.y < minY + y)
	{
		return -1.0f;
	}
	else
	{
		return 0.0f;
	}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedZAxisCylinder::HitByBounded(
		Line *line,
		Scalar enters,
		Scalar leaves
	)
{
	Check(this);
	Check(line);

	Verify(enters <= leaves);
	Verify(leaves >= 0.0f);

	Scalar x = (minX + maxX) * 0.5f;
	Scalar y = (minY + maxY) * 0.5f;
	Scalar radius = maxX - x;

	Scalar a =
		line->direction.y*line->direction.y
		 + line->direction.x*line->direction.x;

	//
	//--------------------------------------------------------------------------
	// If the line is parallel to the cylinder, see if it hits the bottom/top of
	// the cylinder in the X-Z plane.  If it does, it hits the cylinder when it
	// hits the box, otherwise it misses altogether
	//--------------------------------------------------------------------------
	//
	if (Small_Enough(a))
	{
		y -= line->origin.y;
		x -= line->origin.x;
		if (y*y + x*x > radius*radius)
		{
			return False;
		}
		line->length = Max(enters, 0.0f);
		return True;
	}

	//
	//-----------------------------------------------------------------------
	// The line is not parallel to the cylinder, so solve the equation giving
	// the intersection of the line with the cylinder using the quadratic
	// equation
	//-----------------------------------------------------------------------
	//
	y = line->origin.y - y;
	x = line->origin.x - x;

	Scalar b = 2.0f * (line->direction.y*y + line->direction.x*x);
	Scalar c = y*y + x*x - radius*radius;
	Scalar i = b*b - 4.0f*a*c;

	if (i < SMALL)
	{
		return False;
	}

	//
	//-------------------------------------------------------------------------
	// If the line does hit the cylinder, update the enter/leaving distances to
	// take the cylinder surface into account
	//-------------------------------------------------------------------------
	//
	Verify(a > SMALL);
	Scalar ratio = Sqrt(a)/(a*a);
	i = Sqrt(i);
	b *= -ratio;
	i *= ratio;
	Scalar enter = (b - i) * 0.5f;
	if (enter > enters)
	{
		enters = enter;
	}

	Scalar leave = enter + i;
	if (leave < leaves)
	{
		leaves = leave;
	}

	//
	//------------------------------------------------------------------------
	// If we have pushed the entering distance after the leaving distance, the
	// cylinder was missed, otherwise it is hit at the new entering distance
	//------------------------------------------------------------------------
	//
	if (enters > leaves || enters > line->length || leaves < 0.0f)
	{
		return False;
	}

	line->length = Max(enters, 0.0f);
	return True;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
Logical
	BoxedZAxisCylinder::TestInstance() const
{
	return solidType == ZAxisCylinderType;
}