#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)
		: a(a),b(b),c(c),d(d),e(e),f(f),g(g),h(h),i(i),j(j)
	{
		this->shader = shader;
	};

	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
