//To render this file correctly, you have two choices:
// (1) Turn path_tracing off and use radiosity with an offical POV-Ray version
// (2) Make the following changes to lighting.cpp based on
// 	the POV-Ray 3.6.1 source:
//     (Don't worry, cut-and-paste.)
//
//  In function Reflect():
//  Line 2714:
//      VECTOR e2, e3;
//      DBL x1, x2, r1, r2, y1;
//
//  Line 2733:
//      e3[X]=0.12590125091205;
//      e3[Y]=-0.951025125012;
//      e3[Z]=-0.9192051285012;
//      VCross(e2,Jitter_Normal,e3);
//      VNormalize(e3,e2);
//      VCross(e2,Jitter_Normal,e3);
//      r1=sqrt(drand48());
//      r2=2.0*M_PI*drand48();
//      x1=r1*cos(r2);
//      x2=r1*sin(r2);
//      y1=sqrt(1.0-r1*r1);
//      VLinComb3(New_Ray.Direction, x1, e2, x2, e3, y1, Jitter_Normal);
//
//  Now just compile and you've got yourself a path-tracer!
//
//  (Actually, this just hacks in diffuse reflection in the place of
//   specular reflection, but the result is something like path-tracing
//   capability.)


#include "transforms.inc"


#declare path_tracing = on;
#declare focal_blur = on;

//Note:  If you don't want the focal blur, that's fine, but
//you should just set aperture to ~0 for the path tracing since
//it seems to be the simplest way to get good sampling.  I guess
//you could use +AM2 +R6 to get the samples, but it's probably slower.

// Scene setup:
//light_source{<25,15,-35> rgb 1} //(for debug)
camera{
    location <1.5,1.8,-3>
    look_at <0.5,0,0>
    #if(focal_blur=on)
	focal_point vnormalize(<1,2,-3>)*0.9
	aperture 0.3
	blur_samples 2000
	variance 1/10000 //optimize?
	confidence 0.99999
    #end
}
sky_sphere{
    pigment{
	gradient y
	color_map{
	    [0 rgb 0]
	    [1 rgb <0.9,0.9,1.0>*0.7]
	}
    }
}
plane{y,-1.1
    pigment{rgb 1} 
    #if(path_tracing=on)
	finish{ambient 0 diffuse 0 reflection{rgb <1.0,1.0,1.0>}}
    #else
	finish{ambient 0 diffuse 1.0 reflection 0}
    #end
    normal{granite 0.4 scale 0.7}
}
//BIG sphere off-camera to provide some smooth lighting.
sphere{<25,15,5>,10 pigment{rgb <1.0,0.8,0.6>} finish{ambient 3}}
global_settings{
    max_trace_level 8      //Seems low, but it works.  
                           //Equivalent of recursion_limit for path-tracing
    #if(path_tracing=off)
	//I'm SURE you can optimize this quite a bit.  I'm no
	//good at optimizing radiosity.  Here's proof:
	radiosity{ //(untested):
	    normal on
	    count 400
	    error_bound 0.03
	    recursion_limit 4
	    pretrace_start 0.08
	    pretrace_end 0.01
	}
    #end
}


//==================MACROS================//

//Hermite spline:
//  evaluate a spline
//  pv = point vector
//  vv = velocity vector
//  tv = parameterize the spline into segments
//       (should be equally spaced C1 continuity)
//  evaluated at tval
#macro hermite(pv,vv,tv,tval)
    #local l=dimension_size(pv,1);
    #local cntr=0;
    #while(cntr<l & tval>tv[cntr+1])
	#local cntr=cntr+1;
    #end
    #local tvaln=(tval-tv[cntr])/(tv[cntr+1]-tv[cntr]);
    #local p1 = pv[cntr];
    #local p2 = pv[cntr+1];
    #local v1 = vv[cntr];
    #local v2 = vv[cntr+1];
    #local asp=2.0*p1+v1-2.0*p2+v2;
    #local bsp=-3.0*p1-2.0*v1+3.0*p2-v2;
    (asp*tvaln*tvaln*tvaln+bsp*tvaln*tvaln+v1*tvaln+p1)
#end

//Hermite spline velocity vector:
//  Same equations, just return a velocity instead.
#macro hermite_vel(pv,vv,tv,tval)
    #local l=dimension_size(pv,1);
    #local cntr=0;
    #while(cntr<l & tval>tv[cntr+1])
	#local cntr=cntr+1;
    #end
    #local tvaln=(tval-tv[cntr])/(tv[cntr+1]-tv[cntr]);
    #local p1 = pv[cntr];
    #local p2 = pv[cntr+1];
    #local v1 = vv[cntr];
    #local v2 = vv[cntr+1];
    #local asp=2.0*p1+v1-2.0*p2+v2;
    #local bsp=-3.0*p1-2.0*v1+3.0*p2-v2;
    (3.0*asp*tvaln*tvaln+2*bsp*tvaln+v1)
#end

//Smooth quad based on four points, four normals
#macro squad(p1,n1,p2,n2,p3,n3,p4,n4)
    smooth_triangle{p1,n1,p2,n2,p3,n3}
    smooth_triangle{p1,n1,p3,n3,p4,n4}
#end
//Flat quad with four points
#macro quad(p1,p2,p3,p4)
    triangle{p1,p2,p3}
    triangle{p1,p3,p4}
#end

//Create a transformation matrix from a coordinate system:
//  (a,b,c) at point d
#macro v2m(a,b,c,d)
    transform{matrix <a.x,a.y,a.z,b.x,b.y,b.z,c.x,c.y,c.z,d.x,d.y,d.z>}
#end


//==================SPHERE================//

//One big long declaration:
#declare orb = union{

#declare ss=0.12;                     //square size
#declare ts=0.17;                     //triangle size

#declare srot = -40;                  //square rotation (degrees)
#declare trot = 60;                   //triangle rotation

#declare s1 = vrotate(x,y*srot);      //square coord system:
#declare s2 = y;
#declare s3 = vrotate(z,y*srot);

#declare t2 = vnormalize(<1,1,1>);    //triangle coord system:
#declare t1 = vaxis_rotate(vnormalize(<-1,0,1>),t2,trot);
#declare t3 = vnormalize(vcross(t1,t2));

#declare tp = vnormalize(<1,1,1>)*1.2;//triangle center point
#declare sp = 1.0*y;                 //square center point

#declare wid=0.08;                    //object thickness

//Just the top and bottom facets of the square part
#declare square = 
    union{
	triangle{<-ss,0,-ss>,<ss,0,-ss>,<ss,0,ss>}
	triangle{<-ss,0,-ss>,<-ss,0,ss>,<ss,0,ss>}
	triangle{<-ss,-wid,-ss>,<ss,-wid,-ss>,<ss,-wid,ss>}
	triangle{<-ss,-wid,-ss>,<-ss,-wid,ss>,<ss,-wid,ss>}
    }

//Corner points of the triangle:
//     tp3
//    /   \
//   /     \
// tp1-----tp2
#declare tp1 = <-ts*sqrt(3)*0.5,0,-ts*0.5>;
#declare tp2 = < ts*sqrt(3)*0.5,0,-ts*0.5>;
#declare tp3 = <0,0,ts>;

//Scalar multiple for the spline tangent (velocity) vectors (bigger=more curvature):
#declare vs = 1.02;

//mesh the bridge between triangle and square faces.  Really, just create one for a triangle.
//Then copy it three times for that triangle.  Then copy the triangle, with three bridges
//attached, eight times.  Nothing more.
#declare bridge = 
    mesh{
	//Stay with me.  To get the square corner points in the triangle coordinate system,
	//we must first transform the point in square coordinates to absolute coordinates.
	//Next, inverse-transform the point in absolute coordinates to get its representation
	//in triangle coordinates.  Then if you draw a point in the triangle's coordinate
	//system, it will show up at the corner of a square.  Just basic matrix math.

	//two relevant corners of the square (rotated AND translated):
	#declare s1p = vinv_transform( vtransform(<ss,0,ss>,v2m(s1,s2,s3,y)),v2m(t1,t2,t3,tp));
	#declare s2p = vinv_transform( vtransform(<-ss,0,ss>,v2m(s1,s2,s3,y)),v2m(t1,t2,t3,tp));

	//Tangent vectors corresponding to those corners (note that directional vectors
	//are only rotated, NOT translated):
	#declare sv1p = vinv_transform( vtransform(z,v2m(s1,s2,s3,<0,0,0>)),v2m(t1,t2,t3,<0,0,0>));
	#declare sv2p = vinv_transform( vtransform(z,v2m(s1,s2,s3,<0,0,0>)),v2m(t1,t2,t3,<0,0,0>));

	//Normal vector for the square:
	#declare svn = vinv_transform( vtransform(y,v2m(s1,s2,s3,<0,0,0>)),v2m(t1,t2,t3,<0,0,0>));

	//tangent vectors at the two corners of the triangle from which the bridge is drawn:
	#declare tvec1 = <sqrt(3)/2,0,0.5>;
	#declare tvec2 = <sqrt(3)/2,0,0.5>;

	//Calculate four corner points for the square, setting two in by the feature width:
	#declare sqp1 = s1p;
	#declare sqp2 = s2p;
	#declare sqp3 = s1p-svn*wid;
	#declare sqp4 = s2p-svn*wid;
	//Pick some nice curvature values for the four edges of the bridge:
	#declare sqv1 = -sv1p*2.0*vs;
	#declare sqv2 = -sv2p*3.1*vs;
	#declare sqv3 = -sv1p*2.0*vs;
	#declare sqv4 = -sv2p*3.1*vs;
	//Permute how the edges match up to make it a little more interesting:
	#declare c=0;
	#while(c<1)
	    #declare tempp=sqp1;
	    #declare sqp1=sqp2;
	    #declare sqp2=sqp4;
	    #declare sqp4=sqp3;
	    #declare sqp3=tempp;
	    #declare tempv=sqv1;
	    #declare sqv1=sqv2;
	    #declare sqv2=sqv4;
	    #declare sqv4=sqv3;
	    #declare sqv3=tempv;
	    #declare c=c+1;
	#end
	//Parameterize the spline.  Simple, just 0 to 1 for a two-point spline:
	#declare tv = array[2]{0.0,1.0}

	//Create four point arrays to describe the four splines:
	#declare pv1 = array[2]{tp3,sqp1}
	#declare pv2 = array[2]{tp2,sqp2}
	#declare pv3 = array[2]{tp3-y*wid,sqp3}
	#declare pv4 = array[2]{tp2-y*wid,sqp4}

	//Create four tangent vector arrays:
	#declare vv1 = array[2]{tvec1*2.0*vs,sqv1} //inner
	#declare vv3 = array[2]{tvec1*1.8*vs,sqv3} //inner
	#declare vv2 = array[2]{tvec2*2.8*vs,sqv2} //outer
	#declare vv4 = array[2]{tvec2*3.0*vs,sqv4} //outer
	
	//Split the bridge up into (np) segments:
	#declare np=30;
	#declare c=0;
	#while(c<np)
	    #declare tval=c/(np-1);
	    //Get the spline values at this point 
	    // (the four should form a rectangle
	    // around the perimeter of the bridge):
	    #declare p1 = hermite(pv1,vv1,tv,tval);
	    #declare p2 = hermite(pv2,vv2,tv,tval);
	    #declare p3 = hermite(pv3,vv3,tv,tval);
	    #declare p4 = hermite(pv4,vv4,tv,tval);

	    //Get the tangent vectors:
	    #declare v1 = hermite_vel(pv1,vv1,tv,tval);
	    #declare v2 = hermite_vel(pv2,vv2,tv,tval);
	    #declare v3 = hermite_vel(pv3,vv3,tv,tval);
	    #declare v4 = hermite_vel(pv4,vv4,tv,tval);

	    //Create normal vectors for the four sides:
	    #declare nt = vnormalize(vcross(p1-p2,v1));
	    #declare nb = vnormalize(vcross(p3-p4,v3));
	    #declare nl = vnormalize(vcross(p1-p3,v1));
	    #declare nr = vnormalize(vcross(p2-p4,v2));

	    //As long as it's not the 0th segment, then
	    //we have some previous values from which to
	    //join points from this station.  So create
	    //the smooth quades based on points and normal
	    //vectors.
	    #if(c>0)
		squad(p10,nt0,p20,nt0,p2,nt,p1,nt) //top
		squad(p30,nb0,p40,nb0,p4,nb,p3,nb) //bot
		squad(p10,nl0,p1,nl,p3,nl,p30,nl0) //left
		squad(p20,nr0,p2,nr,p4,nr,p40,nr0) //right
	    #end

	    //Store points from this station for reference
	    //on the next iteration of the loop.
	    #declare p10 = p1;
	    #declare p20 = p2;
	    #declare p30 = p3;
	    #declare p40 = p4;
	    #declare nt0 = nt;
	    #declare nb0 = nb;
	    #declare nl0 = nl;
	    #declare nr0 = nr;
	    #declare c=c+1;
	#end
    }

//Now we have enough to declare the triangle part.
//It just has the basic top- and bottom-facet triangle
//along with three bridges, rotated 120 degrees apart:
#declare tri = 
    union{
	triangle{tp1,tp2,tp3}
	triangle{tp1-y*wid,tp2-y*wid,tp3-y*wid}
	object{bridge}
	object{bridge rotate y*120}
	object{bridge rotate y*240}

	#if(0) //draw the local coordinate system
	    cylinder{0,x/2,0.01 pigment{rgb x}}
	    cylinder{0,y/2,0.01 pigment{rgb y}}
	    cylinder{0,z/2,0.01 pigment{rgb z}}
	#end
    }

    union{
	//Main object declaration.  I'm glad this has symmetry
	//in the x- y- and z-planes, otherwise this would
	//sure be a lot uglier.

	//Draw six squares:
	object{square v2m(s1,s2,s3,sp)}
	object{square v2m(s1,s2,s3,sp) rotate x*90}
	object{square v2m(s1,s2,s3,sp) rotate -x*90}
	object{square v2m(s1,s2,s3,sp) rotate z*90}
	object{square v2m(s1,s2,s3,sp) rotate -z*90}
	object{square v2m(s1,s2,s3,sp) rotate x*180}
	//And eight triangles:
	object{tri v2m(t1,t2,t3,tp)}
	object{tri v2m(t1,t2,t3,tp) rotate y*90}
	object{tri v2m(t1,t2,t3,tp) rotate y*180}
	object{tri v2m(t1,t2,t3,tp) rotate y*270}
	object{tri v2m(t1,t2,t3,tp) scale <1,-1,-1>}
	object{tri v2m(t1,t2,t3,tp) rotate y*90 scale <1,-1,-1>}
	object{tri v2m(t1,t2,t3,tp) rotate y*180 scale <1,-1,-1>}
	object{tri v2m(t1,t2,t3,tp) rotate y*270 scale <1,-1,-1>}

	#if(path_tracing=on)
	    pigment{rgb 1} //doesn't matter, but at least it won't complain
	    finish{ambient 0 diffuse 0 reflection <1.0,0.9,0.8>}
	#else
	    pigment{rgb <1.0,0.9,0.8>}
	    finish{ambient 0 diffuse 1}
	#end
    }
}

//Place the spheres:
object{orb rotate -x*20 translate -x*5+z*4}
object{orb rotate z*99+y*35 translate x*2+z*3}
object{orb rotate y*105}

//Light the inside:
sphere{-x*5+z*4,0.6 pigment{rgb <1,0.5,0.5>} finish{ambient 6}}
sphere{x*2+z*3,0.6 pigment{rgb <0.4,0.7,0.9>} finish{ambient 6}}
//sphere{0,0.6 pigment{rgb <0.5,0.5,1>} finish{ambient 6}}