#ifndef QUADRIC_HXX
 #define QUADRIC_HXX
 
 #include "Primitive.hxx"
 #include "Box.hxx"
 
 class Quadric : public Primitive
 {
 	float a, b, c, d, e, f, g, h, i, j,
 	      alpha, beta, gamma, p, radicand, t[2];
 	int   num;
 public:
 	Quadric(float a, float b, float c, float d, float e, float f, float g, 
 	        float h, float i, float j, Shader *shader, PhotonDistributor *photon_distributor)
 		: Primitive(shader, photon_distributor),a(a),b(b),c(c),d(d),e(e),f(f),g(g),h(h),i(i),j(j)
 	{};
 
 	virtual ~Quadric(){};
 
 	virtual bool Intersect(Ray &ray)
 	{
 		alpha =       a * ray.dir.x * ray.dir.x
 		      +       b * ray.dir.y * ray.dir.y
 		      +       c * ray.dir.z * ray.dir.z
 		      +       d * ray.dir.x * ray.dir.y
 		      +       e * ray.dir.x * ray.dir.z
 		      +       f * ray.dir.y * ray.dir.z;
 
 		beta  = 2.0 * a * ray.org.x * ray.dir.x
 		      + 2.0 * b * ray.org.y * ray.dir.y
 		      + 2.0 * c * ray.org.z * ray.dir.z
 		      +       d * ray.org.x * ray.dir.y
 		      +       d * ray.org.y * ray.dir.x
 		      +       e * ray.org.x * ray.dir.z
 		      +       e * ray.org.z * ray.dir.x
 		      +       f * ray.org.y * ray.dir.z
 		      +       f * ray.org.z * ray.dir.y
 		      +       g *             ray.dir.x
 		      +       h *             ray.dir.y
 		      +       i *             ray.dir.z;
 
 		gamma =       a * ray.org.x * ray.org.x
 		      +       b * ray.org.y * ray.org.y
 		      +       c * ray.org.z * ray.org.z
 		      +       d * ray.org.x * ray.org.y
 		      +       e * ray.org.x * ray.org.z
 		      +       f * ray.org.y * ray.org.z
 		      +       g * ray.org.x
 		      +       h * ray.org.y
 		      +       i * ray.org.z
  		      +       j;
 
 		// Apply the p-q-Formula:
 		p        = beta / (2.0 * alpha);
 		radicand = p * p - gamma / alpha;
 
 		if (radicand >= 0.0)
 		{
 			// We have two solutions. The Sqrt is only calculated once:
 			radicand = sqrt(radicand);
 			t[0]     = -p + radicand;
 			t[1]     = -p - radicand;
 
 			// Check the closer intersection point and take it if necessary
 			num = t[1] < t[0];
 			if (t[num] < ray.t && t[num] > Epsilon)
 			{
 				ray.t = t[num];
 				ray.hit = this;
 				return true;
 			}
 
 			// Check the farer intersection point and take it if necessary
 			num = !num;
 			if (t[num] < ray.t && t[num] > Epsilon)
 			{
 				ray.t = t[num];
 				ray.hit = this;
 				return true;
 			}
 		}
 
 		// We could not find any intersection point
 		return false;
 	};
 
 	virtual Vec3f GetNormal(Ray &ray)
 	{
 		Vec3f p      = ray.org + ray.t * ray.dir;
 		Vec3f normal = Vec3f(2 * a * p.x + d * p.y + e * p.z + g,
 		                     2 * b * p.y + d * p.x + f * p.z + h,
 		                     2 * c * p.z + e * p.x + f * p.y + i);
 		Normalize(normal);
 		return normal;
 	};
 
 	virtual Box CalcBounds()
 	{
 		// Return real bounding boxes later on (TODO)
 		return Box();
 	}
 };
 
 #endif