#ifndef SCHLICKPHOTONDISTRIBUTOR_HXX
 #define SCHLICKPHOTONDISTRIBUTOR_HXX
 
 #include "PhotonDistributor.hxx"
 
 class SchlickPhotonDistributor : public PhotonDistributor
 {
 private:
 	Randomizer randomizer;
 public:
 
 	Vec3f	   hit_point;           // the hit point on the surface
 	float      p12,                 // square of the isotropy factor for primary surface
 		   p22;                 // square of the isotropy factor for primary surface
 	float      mirror,              // Material based mirroring    percentage
 		   glossy_mirror,       // Material based glossy reflection percentage
 		   glossy_transparency, // Material based glossy refraction percentage
 		   transparency;        // Material based transparency percentage
 
 	SchlickPhotonDistributor(Scene *scene, Material *material)
 		: PhotonDistributor(scene, material)
 // Optionally fix the random seed for parallization
 #ifdef RANDOMSEED
 		, randomizer(RANDOMSEED)
 #endif
 	{
 		// Initialization of reflection and refraction
 		if (material->mirror > 1.0 || material->transparency > 1.0 || material->mirror + material->transparency > 1.0)
 			cerr << "WARNING: Material '" << material << "' has wrong reflection/refraction percentages (> 1.0)!" << endl;
 
 		mirror	      = (material->mirror == UNDEF) ? 0.0f : material->mirror;
 		glossy_mirror       = (material->glossy_mirror       == UNDEF) ? 0.0f : material->glossy_mirror;
 		glossy_transparency = (material->glossy_transparency == UNDEF) ? 0.0f : material->glossy_transparency;
 		transparency	= (material->transparency == UNDEF) ? 0.0f : material->transparency;
 
 		// Precomputation of p^2
 		p12 = material->p1 * material->p1;
 		p22 = material->p2 * material->p2;
 	};
 
 	virtual ~SchlickPhotonDistributor()
 	{};
 
 	// As this is a rather common method, it is only maintained once
 	#include "SchlickTerm.include"
 
 	virtual void FancyDistribution(TravellingPhoton &photon, vector<Photon*> &list, const Vec3f &normal, const bool &into, const float &what_to_do)
 	{
 		float m = mirror,
 		      t = transparency;
 
 		// Get the real distribution of the mirror coefficient first
 		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 		float cosa	     = 1.0,
 		      cosb	     = 1.0,
 		      cosb2,
 		      new_refrac_index = 1.0,
 		      c,
 		      n		= 1.0,
 		      schlick;
 
 		if (t)
 		{
 			// Cosine of the incident angle
 			cosa = Dot(normal, -photon.dir);
 
 			if (into)
 				new_refrac_index = photon.refrac_index + (material->refrac_index - 1.0f);
 			else
 				new_refrac_index = photon.refrac_index - (material->refrac_index - 1.0f);
 
 			n     = photon.refrac_index / new_refrac_index,
 			cosb2 = 1.0f - n * n * (1.0f - cosa*cosa);
 
 			if (cosb2 < 0.0f)
 				// Total reflection
 				m += t;
 			else
 			{
 				cosb = sqrtf(cosb2),
 				c    = (n - 1) / (n + 1);
 
 				c *= c;
 
 				// We use the angle in the less dense medium for Schlick's term
 				if (into)
 					schlick = c + (1 - c) * powf((1 - cosa), 5);
 				else
 					schlick = c + (1 - c) * powf((1 - cosb), 5);
 
 				schlick *= t;
 
 				m += schlick;
 				t -= schlick;
 			}
 		}
 
 		// Reset the photon's path
 		photon.org = photon.org + photon.t * photon.dir;
 		photon.t   = Infinity;
 		photon.hit = NULL;
 
 		// Evaluate against the random parameter
 		// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 		if (what_to_do + t > 1.0f)
 		{
 			// Transparency
 			// ~ ~ ~ ~ ~ ~
 			photon.dir = (n * cosa - cosb) * normal + n * photon.dir;
 
 			// Apply refraction index now as we are sure to have changed the medium
 			photon.refrac_index = new_refrac_index;
 
 			if (material->refrac_filter != UNDEF)
 				photon.energy = Product(photon.energy, material->refrac_filter       * material->color
 								     + (1 - material->refrac_filter) * Vec3f(1.0f));
 		}
 		else
 		{
 			// Mirror
 			// ~ ~ ~
 			// Project the photon direction onto the normal vector...
 			Vec3f projected = Dot(normal, -photon.dir) * normal;
 
 			// ...and mirror away
 			photon.dir += 2 * projected;
 			Normalize(photon.dir);
 
 			if (material->mirror_filter != UNDEF)
 				photon.energy = Product(photon.energy, material->mirror_filter       * material->color
 								     + (1 - material->mirror_filter) * Vec3f(1.0f));
 		}
 
 		scene->PhotonTrace(photon, list);
 	}
 
 	virtual void Distribute(TravellingPhoton &photon, vector<Photon*> &list)
 	{
 
 		Vec3f H,
 		      N = photon.hit->GetNormal(photon);
 		float t,		 // +
 		      u,		 //  |
 		      v,		 //   > cosines of various angles between normal vectors
 		      x,		 //  |
 		      w;		 // +
 		bool  into;
 
 		into = true;
 		if (Dot(N,photon.dir) > 0)
 		{
 			N = -N;
 			into = false;
 		}
 
 		hit_point = photon.org + photon.t * photon.dir;
 
 		for (vector<Procedural*>::iterator it = material->procedural.begin(); it != material->procedural.end(); ++it)
 			if ((*it)->ProvidesBumpMap())
 				(*it)->BumpMap(N, hit_point);
 
 		float what_to_do = randomizer.drand();
 
 		if (what_to_do < 1.0f - mirror - transparency)
 		{
 			// Here, we do some approximation for specular reflectance:
 			// - We choose a random reflection direction within both the upper hemisphere
 			//   and the incident plane
 			// - Then, we calculate the BRDF according to this point (i.e. the inverse BRDF!)
 			// - The photon is now split in two: One is reflected carrying to the specular energy,
 			//   while the other is absorbed by the surface keepig the remianing energy for the
 			//   photon map
 
 			// Get a new random direction
 			// ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 		       // Randomize the parametrized direction
 			float incident_part = randomizer.drand() * 2 - 1,
 			      normal_part   = randomizer.drand();
 
 			// Span the incident plane
 			Vec3f incident_direction,
 			      new_direction;
 
 			incident_direction = photon.dir + Dot(N, -photon.dir) * N;
 
 			if (Nonzero(incident_direction))
 			{
 				Normalize(incident_direction);
 
 				// Deparametrize this direction
 				new_direction = incident_part * incident_direction + normal_part * N;
 				Normalize(new_direction);
 			}
 			else
 				new_direction = N;
 
 			// Calculate the inverse BRDF and act accordingly
 			// We do this in ways of the Schlick model (Compare Schlick Shader for further details)
 			// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 			float specular;
 
 			// Get the halfway vector in between incoming and outgoing ray
 			H = new_direction - photon.dir;
 			if (Nonzero(H))
 			{
 				Normalize(H);
 
 				// Calculation of various scaling factors
 				t = fabs(Dot(H,	     N));
 				u = fabs(Dot(-photon.dir,   H));
 				v = fabs(Dot(-photon.dir,   N));
 				x = fabs(Dot(new_direction, N));
 				w = fabs(Dot(H - N,	 N));
 
 				// Precomputation of w^2
 				w = w * w;
 
 				specular = SchlickTerm(t, u, v, x, w);
 			}
 			else
 				specular = 0.0f;
 
 			// Now perform the split and make sure that the reflected photon carries color information etc.
 			Vec3f pixelColor = material->color;
 
 			for (vector<Procedural*>::iterator it = material->procedural.begin(); it != material->procedural.end(); ++it)
 				if ((*it)->ProvidesTexture())
 					pixelColor = Product(pixelColor, (*it)->GetColor(hit_point));
 
 			photon.org = hit_point;
 
 			if (glossy_mirror < 1.0)
 			{
 				// A new slim photon is created to be inserted into the photon map
 				Photon *absorbed_photon = new Photon(photon);
 				absorbed_photon->energy *= (1.0f - glossy_mirror) * (1.0f - specular);
 
 				list.push_back(absorbed_photon);
 			}
 
 			if (specular)
 			{
 				// The old photon is used to follow the further path
 				photon.energy = Product((glossy_mirror + specular - glossy_mirror * specular) * photon.energy, pixelColor);
 
 				if (photon.energy.x + photon.energy.y + photon.energy.z > MINIMAL_PHOTONENERGY)
 				{
 					photon.dir = new_direction;
 					photon.t   = Infinity;
 					photon.hit = NULL;
 
 					scene->PhotonTrace(photon, list);
 				}
 			}
 		}
 		else 
 			// Any other effect is performed
 			FancyDistribution(photon, list, N, into, what_to_do);
 
 	};
 };
 
 #endif