#ifndef ABSTRACTATMOSPHERE_HXX
#define ABSTRACTATMOSPHERE_HXX

#include "SimpleParticle.hxx"
#include "Rainbow.hxx"

/**
   An abstract superclass of any atmosphere.
   Supports volume tracing. The intensities along the ray are integrated.
   The intensity at a given point in space is determined by a trilinear
   interpolation between the 8 neighbouring points in a regular grid.
   Furthermore, rainbows can be generated.
*/
class AbstractAtmosphere
{
 public:
    virtual ~AbstractAtmosphere() {
	// delete all particles in the atmosphere
	for (int x=0; x<=resX; x++)
	    for (int y=0; y<=resY; y++)
		for (int z=0; z<=resZ; z++)
		    delete particles[x][y][z];

	// delete all rainbows
	for (unsigned int i=0; i<rainbows.size(); i++) delete rainbows[i];
    }

    /**
       Combines the given color of the scene with the color of the atmosphere
       that is obtained by an integration along the ray.
    */
    Vec3f rayTrace( const Vec3f &sceneColor, const Ray &ray ) const {
	// front to back approach
	// initialize the color and transparency
	Vec3f color(0,0,0);
	float transparency = 1.0;
	for (float t=0; t<maxT*deltaT; t+=deltaT) {
	    if (t >= ray.t) break; // break when intersection point is reached
	    Vec3f point = ray.org + t * ray.dir; // compute point in space
	    // get the transparency at this point (by a trilinear interpolation)
	    float transpAtPoint = transparencyAtPoint(point);
	    // get the color at this point (by a trilinear interpolation)
	    Vec3f colAtPoint = colorAtPoint(point, ray);
	    // update the color and transparency
	    color += transparency * colAtPoint * (1 - transpAtPoint);
	    transparency *= transpAtPoint;

	    // break if you can hardly see the following things
	    if (transparency < 0.01) break; 
	}
	
	// add the original color of the scene
	// (weighted by the remaining transparency at the intersection point)
	color += transparency * sceneColor;

	// add the colors of all rainbows
	for (unsigned int i=0; i<rainbows.size(); i++)
	    color = rainbows[i]->rayTrace(color,ray);
	
	color.assertInterval(0,1);
	return color;
    }
    
 protected:
    static const int resX = 100;
    static const int resY = 100;
    static const int resZ = 100;
    
    float deltaT; // step size of the volume tracing
    int maxT; // upper bound for the volume tracing

    Vec3f min; // min point of the volume
    Vec3f max; // max point of the volume

    /** the particles in the atmosphere (ordered in a regular grid) */
    AbstractParticle* particles[resX+1][resY+1][resZ+1];

    /** the rainbows (primary and secondary rainbows) */
    std::vector<Rainbow*> rainbows;
    
    /** Returns the transparency at the given point. */
    float transparencyAtPoint(const Vec3f &point) const {
	if ( ! ((point >= min) && (point < max))) return 1;
	float weights[8];
	// get the coordinates of the lower left neighbour
	// and compute the weights used for the trilinear interpolation
	Vec3f lowerLeft = interpolate(point, weights);
	int llX = (int)lowerLeft.x;
	int llY = (int)lowerLeft.y;
	int llZ = (int)lowerLeft.z;
	float result = 0;
	for (int x=0; x<=1; x++)
	    for (int y=0; y<=1; y++)
		for (int z=0; z<=1; z++)
		    result += weights[wi(x,y,z)] * particles[llX+x][llY+y][llZ+z]->transparency();
	return result;
    }

    /** Returns the color at the given point. */
    Vec3f colorAtPoint(const Vec3f &point, const Ray &ray) const {
	if ( ! ((point >= min) && (point < max))) return Vec3f(0,0,0);
	float weights[8];
	// get the coordinates of the lower left neighbour
	// and compute the weights used for the trilinear interpolation
	Vec3f lowerLeft = interpolate(point, weights);
	int llX = (int)lowerLeft.x;
	int llY = (int)lowerLeft.y;
	int llZ = (int)lowerLeft.z;
	Vec3f result(0,0,0);
	for (int x=0; x<=1; x++)
	    for (int y=0; y<=1; y++)
		for (int z=0; z<=1; z++)
		    result += weights[wi(x,y,z)] * particles[llX+x][llY+y][llZ+z]->color(ray);
	return result;
    }

    /**
       Performs a trilinear interpolation for the given point.
       The coordinates of the lower left neighbour are returned
       and the weights of the 8 neighbours are computed and set.
    */
    Vec3f interpolate(const Vec3f &point, float *weights) const {
	// determine lower left neighbour
	Vec3f posInArray = (point - min) / (max - min) * Vec3f(resX,resY,resZ);
	Vec3f lowerLeft;
	lowerLeft.x = (int)posInArray.x;
	lowerLeft.y = (int)posInArray.y;
	lowerLeft.z = (int)posInArray.z;
	Vec3f posRelative = posInArray - lowerLeft;
	// interpolate in x-direction
	for (int y=0; y<=1; y++)
	    for (int z=0; z<=1; z++) {
		weights[wi(0,y,z)] = 1 - posRelative.x;
		weights[wi(1,y,z)] = posRelative.x;
	    }
	// interpolate in y-direction
	for (int x=0; x<=1; x++)
	    for (int z=0; z<=1; z++) {
		weights[wi(x,0,z)] *= 1 - posRelative.y;
		weights[wi(x,1,z)] *= posRelative.y;
	    }
	// interpolate in z-direction
	for (int x=0; x<=1; x++)
	    for (int y=0; y<=1; y++) {
		weights[wi(x,y,0)] *= 1 - posRelative.z;
		weights[wi(x,y,1)] *= posRelative.z;
	    }
	return lowerLeft;
    }

    /** weights index (combines three indices into one) */
    int wi(int x, int y, int z) const {
	return x + 2*y + 4*z;
    }

    /** Initializes the particles in the atmosphere. */
    virtual void initParticles() = 0;
};


#endif // ABSTRACTATMOSPHERE_HXX