```///////////////////////////////////////
//
//  QJuliaFragment.cg
//  4/17/2004
//
//
//  Intersects a ray with the qj set w/ parameter mu and returns
//  the color of the phong shaded surface (estimate)
//  (Surface also colored by normal direction)
//
//     Keenan Crane (kcrane@uiuc.edu)
//
//
//

// Some constants used in the ray tracing process.  (These constants
// were determined through trial and error and are not by any means
// optimal.)

#define BOUNDING_RADIUS_2  3.0      // radius of a bounding sphere for the set used to accelerate intersection

#define ESCAPE_THRESHOLD   1e1      // any series whose points' magnitude exceed this threshold are considered
// divergent

#define DEL                1e-4     // delta is used in the finite difference approximation of the gradient
// (to determine normals)

// --------- quaternion representation -----------------------------------------
//
// Each quaternion can be specified by four scalars q = A + Bi + Cj + Dk, so are
// stored as a float4.  I've tried a struct containing a separate scalar and
// 3-vector to avoid a lot of swizzling, but the float4 representation ends up
// using fewer instructions.  A matrix representation is also possible.
//

// -------- quatMult() ----------------------------------------------------------
//
// Returns the product of quaternions q1 and q2.
// Note that quaternion multiplication is NOT commutative (i.e., q1 ** q2 != q2 ** q1 ).
//

float4 quatMult( float4 q1, float4 q2 )
{
float4 r;

r.x   = q1.x*q2.x - dot( q1.yzw, q2.yzw );
r.yzw = q1.x*q2.yzw + q2.x*q1.yzw + cross( q1.yzw, q2.yzw );

return r;
}

// ------- quatSq() --------------------------------------------------------------
//
// Returns the square of quaternion q.  This function is a special (optimized)
// case of quatMult().
//

float4 quatSq( float4 q )
{
float4 r;

r.x   = q.x*q.x - dot( q.yzw, q.yzw );
r.yzw = 2*q.x*q.yzw;

return r;
}

// ------- iterateIntersect() -----------------------------------------------------
//
// Iterates the quaternion q for the purposes of intersection.  This function also
// produces an estimate of the derivative at q, which is required for the distance
// estimate.  The quaternion c is the parameter specifying the Julia set, and the
// integer maxIterations is the maximum number of iterations used to determine
// whether a point is in the set or not.
//
// To estimate membership in the set, we recursively evaluate
//
// q = q*q + c
//
// until q has a magnitude greater than the threshold value (i.e., it probably
// diverges) or we've reached the maximum number of allowable iterations (i.e.,
// it probably converges).  More iterations reveal greater detail in the set.
//
// To estimate the derivative at q, we recursively evaluate
//
// q' = 2*q*q'
//
// concurrently with the evaluation of q.
//

void iterateIntersect( inout float4 q, inout float4 qp, float4 c, int maxIterations )
{
for( int i=0; i<maxIterations; i++ )
{
qp = 2.0 * quatMult(q, qp);
q = quatSq(q) + c;

if( dot( q, q ) > ESCAPE_THRESHOLD )
{
break;
}
}
}

// ----------- normEstimate() -------------------------------------------------------
//
// Create a shading normal for the current point.  We use an approximate normal of
// the isosurface of the potential function, though there are other ways to
// generate a normal (e.g., from an isosurface of the potential function).
//

float3 normEstimate(float3 p, float4 c, int maxIterations )
{
float3 N;
float4 qP = float4( p, 0 );

float4 gx1 = qP - float4( DEL, 0, 0, 0 );
float4 gx2 = qP + float4( DEL, 0, 0, 0 );
float4 gy1 = qP - float4( 0, DEL, 0, 0 );
float4 gy2 = qP + float4( 0, DEL, 0, 0 );
float4 gz1 = qP - float4( 0, 0, DEL, 0 );
float4 gz2 = qP + float4( 0, 0, DEL, 0 );

for( int i=0; i<maxIterations; i++ )
{
gx1 = quatSq( gx1 ) + c;
gx2 = quatSq( gx2 ) + c;
gy1 = quatSq( gy1 ) + c;
gy2 = quatSq( gy2 ) + c;
gz1 = quatSq( gz1 ) + c;
gz2 = quatSq( gz2 ) + c;
}

return N;
}

// ---------- intersectQJulia() ------------------------------------------
//
// Finds the intersection of a ray with origin rO and direction rD with the
// quaternion Julia set specified by quaternion constant c.  The intersection
// is found using iterative sphere tracing, which takes a conservative step
// along the ray at each iteration by estimating the minimum distance between
// the current ray origin and the closest point in the Julia set.  The
// parameter maxIterations is passed on to iterateIntersect() which determines
// whether the current ray origin is in (or near) the set.
//

float intersectQJulia( inout float3 rO, inout float3 rD, float4 c, int maxIterations, float epsilon )
{
float dist; // the (approximate) distance between the first point along the ray within
// epsilon of some point in the Julia set, or the last point to be tested if
// there was no intersection.

while( 1 )
{

float4 z = float4( rO, 0 );         // iterate on the point at the current ray origin.  We
// want to know if this point belongs to the set.

float4 zp = float4( 1, 0, 0, 0 );   // start the derivative at real 1.  The derivative is
// needed to get a lower bound on the distance to the set.

// iterate this point until we can guess if the sequence diverges or converges.
iterateIntersect( z, zp, c, maxIterations );

// find a lower bound on the distance to the Julia set and step this far along the ray.
float normZ = length( z );
dist = 0.5 * normZ * log( normZ ) / length( zp );  //lower bound on distance to surface

rO += rD * dist; // (step)

// Intersection testing finishes if we're close enough to the surface
// (i.e., we're inside the epsilon isosurface of the distance estimator
// function) or have left the bounding sphere.
if( dist < epsilon || dot(rO, rO) > BOUNDING_RADIUS_2 )
break;
}

// return the distance for this ray
return dist;
}

// ----------- Phong() --------------------------------------------------
//
// Computes the direct illumination for point pt with normal N due to
// a point light at light and a viewer at eye.
//

float3 Phong( float3 light, float3 eye, float3 pt, float3 N )
{
float3 diffuse = float3( 1.00, 0.45, 0.25 ); // base color of shading
const int specularExponent = 10;             // shininess of shading
const float specularity = 0.45;              // amplitude of specular highlight

float3 L     = normalize( light - pt );  // find the vector to the light
float3 E     = normalize( eye   - pt );  // find the vector to the eye
float  NdotL = dot( N, L );              // find the cosine of the angle between light and normal
float3 R     = L - 2 * NdotL * N;        // find the reflected vector

diffuse += abs( N )*0.3;  // add some of the normal to the
// color to make it more interesting

// compute the illumination using the Phong equation
return diffuse * max( NdotL, 0 ) + specularity*pow( max(dot(E,R),0), specularExponent );
}

// ---------- intersectSphere() ---------------------------------------
//
// Finds the intersection of a ray with a sphere with statically
// sphere serves as a bounding volume for the Julia set.

float3 intersectSphere( float3 rO, float3 rD )
{
float B, C, d, t0, t1, t;

B = 2 * dot( rO, rD );
C = dot( rO, rO ) - BOUNDING_RADIUS_2;
d = sqrt( B*B - 4*C );
t0 = ( -B + d ) * 0.5;
t1 = ( -B - d ) * 0.5;
t = min( t0, t1 );
rO += t * rD;

return rO;
}

// ------------ main() -------------------------------------------------
//
//  Each fragment performs the intersection of a single ray with
//  the quaternion Julia set.  In the current implementation
//  the ray's origin and direction are passed in on texture
//  coordinates, but could also be looked up in a texture for a
//  more general set of rays.
//
//  The overall procedure for intersection performed in main() is:
//
//  -move the ray origin forward onto a bounding sphere surrounding the Julia set
//  -test the new ray for the nearest intersection with the Julia set
//  -if the ray does include a point in the set:
//      -estimate the gradient of the potential function to get a "normal"
//      -use the normal and other information to perform Phong shading
//      -cast a shadow ray from the point of intersection to the light
//  -return the shaded color if there was a hit and the background color otherwise
//

float4 main( float3 rO : TEXCOORD0,                // ray origin
float3 rD : TEXCOORD1,                // ray direction (unit length)

uniform float4 mu,                    // quaternion constant specifying the particular set
uniform float epsilon,                // specifies precision of intersection
uniform float3 eye,                   // location of the viewer
uniform float3 light,                 // location of a single point light
uniform int maxIterations )           // maximum number of iterations used to test convergence
: COLOR
{
const float4 backgroundColor = float4( 0.3, 0.3, 0.3, 0 );  //define the background color of the image

float4 color;  // This color is the final output of our program.

// Initially set the output color to the background color.  It will stay
// this way unless we find an intersection with the Julia set.

color = backgroundColor;

// First, intersect the original ray with a sphere bounding the set, and
// move the origin to the point of intersection.  This prevents an
// unnecessarily large number of steps from being taken when looking for
// intersection with the Julia set.

rD = normalize( rD );  //the ray direction is interpolated and may need to be normalized
rO = intersectSphere( rO, rD );

// Next, try to find a point along the ray which intersects the Julia set.
// (More details are given in the routine itself.)

float dist = intersectQJulia( rO, rD, mu, maxIterations, epsilon );

// We say that we found an intersection if our estimate of the distance to
// the set is smaller than some small value epsilon.  In this case we want
// to do some shading / coloring.

if( dist < epsilon )
{
// Determine a "surface normal" which we'll use for lighting calculations.
float3 N = normEstimate( rO, mu, maxIterations );

// Compute the Phong illumination at the point of intersection.
color.rgb = Phong( light, rD, rO, N );
color.a = 1;  // (make this fragment opaque)

// If the shadow flag is on, determine if this point is in shadow
{
// The shadow ray will start at the intersection point and go
// towards the point light.  We initially move the ray origin
// a little bit along this direction so that we don't mistakenly
// find an intersection with the same point again.

float3 L = normalize( light - rO );
rO += N*epsilon*2.0;
dist = intersectQJulia( rO, L, mu, maxIterations, epsilon );

// Again, if our estimate of the distance to the set is small, we say
// that there was a hit.  In this case it means that the point is in