#include "utils.h"

// For point P, gather illumination from the light sources and
// compute the smoke density at that point.  Do not count lights
// with a "__foglight" parameter that's set to 0.
void
smokedensity (point Pcur, 
	      output color Lscatter, output float smoke)
{
    Lscatter = 0;
    lights (Pcur) {
		float foglight = 1;
		getmessage("light", "__foglight",foglight);
        Lscatter += foglight * Cl;
    }
	smoke = 1;
}


// Return a component-by-component exp() of a color
color colorexp (color C)
{
    return color (exp(C[0]), exp(C[1]), exp(C[2]));
}

void swap (output float x, output float y)
{
	float tmp = x; x = y; y = tmp;
}


float quadratic (float a, float b, float c, output float roots[2])
{
    if (a == 0.0) {
        if (b == 0.0)
			return 0;
        else {        
            roots[0] = -c/b;
            return 1;
        }
    } 
	else {
        //  Compute the discriminant d = b^2 - 4ac
        float d = b*b - 4*a*c;
        if (d < 0.0)
			return 0;
		else if (d == 0.0) {
			roots[0] = -b/(2*a);
			return 1;
        } 
		else {
            // Two real roots
			d = sqrt (d);
			roots[0] = - (b - d)/(2*a);
			roots[1] = - (b + d)/(2*a);
			if (roots[0] > roots[1]) // Reorder if needed
				swap (roots[0], roots[1]);
            return 2;
        }
    }
}

volume mango_volume(
	float	opacdensity = 0,
	float	lightdensity = 0.005,
	color	scatter = 1,
	float	integstart = 0,
	float	integend = 10000,
	float	stepsize = 0.1,
	float	minsteps = 3,
	float	maxsteps = 100,  
	float	background_object = 0,
	string	vol_type[] = { "" },
	string	vol_name[] = { "" },
	string	vol_coord[] = { "" },
	float	vol_density[] = { 0 },
	color	vol_color[] = { 0 }
	)
{
    point Worigin = P - I;    // Origin of volume ray

    vector WIN = normalize(I);

	color Cv = 0, Ov = 0;   // color & opacity of volume that we accumulate

	for (int i = 0; i < arraylength (vol_name); i += 1) {
		float raymarchingstart, raymarchingend;
		// compute the intersection of the volume ray with the canonical volume
		point O = transform ("common", vol_coord[i], Worigin);
		vector D = transform ("common", vol_coord[i], WIN);
		if (vol_type[i] == "lightcone") { // intersection with a canonical cone (apex at origin, pointing downward the z axis, height = 1; radius = 1)
			// setup a quadratic equation
			float a = sqr (D[0]) + sqr (D[1]) - sqr (D[2]);
			float b = 2 * (D[0] * O[0] + D[1] * O[1] - D[2] * O[2]);
			float c = sqr (O[0]) + sqr (O[1]) - sqr (O[2]);
			float roots[2];
			if (quadratic (a, b, c, roots) != 2)
				continue;
			// clamp against plane z = 0 && z = -1;
			if (D[2] != 0.0)
			{
				float t0 = - O[2] / D[2];
				float t1 = (-1.0 - O[2]) / D[2];
				if (t0 > t1) // Reorder if needed
					swap (t0, t1);
				roots[0] = max (t0, roots[0]);
				roots[1] = min (t1, roots[1]);
				if (roots[0] >= roots[1])
					continue;
			}
			raymarchingstart = max (integstart, roots[0]);
			raymarchingend = min (integend, roots[1]);
		}
		else
			continue; // unknown volume type

		// Integrate forwards from the start point
		float d = raymarchingstart - random()*stepsize;
		// Calculate a reasonable step size
		float end = min (length (I), raymarchingend) - 0.0001;
		float minss = (end-d) / maxsteps;
		float maxss = (end-d) / minsteps;
		float ss = min (max (minss, stepsize), maxss);

	    float dtau, last_dtau;
		color li, last_li;

		// Get the in-scattered light and the local smoke density for the
		// beginning of the ray
		smokedensity (Worigin + d*WIN, last_li, last_dtau);


		while (d < end) {
			// Take a step and get the local scattered light and smoke density
			ss = clamp (ss, 0.005, end-d);
			d += ss;
			smokedensity (Worigin + d*WIN, li, dtau);
	
			// Find the blocking and light scattering contribution of the
			// portion of the volume covered by this step.
			float tau = opacdensity * (ss/2) * (dtau + last_dtau);
			color lighttau = vol_color[i] * lightdensity * vol_density[i] * ss/2 * (li*dtau + last_li*last_dtau);
	
			// Composite with exponential extinction of background light
			Cv += (1-Ov) * lighttau * scatter;
			Ov += (1-Ov) * (1 - colorexp (-tau*scatter));
			last_dtau = dtau;
			last_li = li;
		}
	}
    
    // C & opacity are the color and opacity of the background element.
    // Now Cv is the light contributed by the volume itself, and Ov is the
    // opacity of the volume, i.e. (1-Ov)*C is the light from the background
    // which makes it through the volume.
    if ( raylevel() > 0 || background_object != 0 ) {
        // Standard compositing model:  volume over surface
        C = Cv + (1-Ov)*C;
        opacity = Ov + (1-Ov)*opacity;
    } else {
        // Currently in Gelato this shader always calculates the volume
        // contribution between P and the camera.  But a closer, non-opaque
        // object also accounts for part of that ray.  This alternate
        // compositing model actually gets the math right in that case!
        // BUT... this only works if the farthest object is either opaque
        // or uses the standard compositing.  So if you use smoke, be sure
        // to put a "back-plane" on the scene (it can be transparent black,
        // but MUST use standard_compositing = 1).
        C = Cv*opacity + (1-Ov)*C;
    }
}	
