/*
	name:	intersct.c

	Intersection-routines
	---------------------

	These routines will calculate the t-parameter of a line,
	which corresponds to the point where the given line intersects with
	the given primitive. The routines return zero (0) if no intersection
	occured.


    This source-code is part of the RayLab 1.1 package, and it is provided
    for your compiling pleasure.  You may use it, change it, re-compile it
    etc., as long as nobody else but you receive the changes/compilations
    that you have made!

    You may not use any part(s) of this source-code for your own productions
    without the permission from the author of RayLab. Please read the legal
    information found in the users documentation for RayLab for more details.

*/


#include  "defs.h"
#include  "extern.h"




/**********************************************************************
 *
 *	Try all objects for intersection, and return object number
 *	and intersection information.
 *
 **********************************************************************/

long Intersect_All_Objects(INTERSECTION *Intersct, LINE *RayLine)
{
	INTERSECTION TmpIntersct;
	double	t, tnew;
	long	i, ObjNum;

	ObjNum=-1L;

	t=-1.0;
	for(i=0;i<NumObjects;i++) {
	    tnew=Intersect_Object(&TmpIntersct, RayLine, ObjectArray[i]);
	    if((tnew>EPSILON)&&((tnew<t)||(t<=EPSILON))) {
		ObjNum=i;
		t=tnew;
		*Intersct=TmpIntersct;
	    }
	}

	return(ObjNum);
}



/**********************************************************************
 *
 *	Get (possible) intersection with a given object
 *
 **********************************************************************/


double Intersect_Object(INTERSECTION *Intrsct, LINE *Line, OBJECT *Object)
{
	LINE	TransformedLine;
	register double	t,a,b,c;
	register long	i;
	VECTOR	tempv;

	if(CheckForBounding(Line, Object)==TRUE) {			/* Check for boundings, if any */
	  TransformedLine.Origin=Line->Origin;
	  TransformedLine.Direction=Line->Direction;

	  if(Object->Transform.NumTransforms>0) {
	    for(i=Object->Transform.NumTransforms-1L;i>=0L;i--) {
		switch(Object->Transform.Entry[i].Type) {
		    case TRANSFORM_SCALE:
			tempv=Object->Transform.Entry[i].Values;
			a=1.0/tempv.x; b=1.0/tempv.y; c=1.0/tempv.z;
			TransformedLine.Origin.x*=a; TransformedLine.Origin.y*=b; TransformedLine.Origin.z*=c;
			TransformedLine.Direction.x*=a; TransformedLine.Direction.y*=b; TransformedLine.Direction.z*=c;
			break;
		    case TRANSFORM_MOVE:
			SubVector((VECTOR *)&TransformedLine.Origin,(VECTOR *)&TransformedLine.Origin,&Object->Transform.Entry[i].Values);
			break;
		    case TRANSFORM_ROTATE:
			NegVector(&tempv,&Object->Transform.Entry[i].Values);
			RevRotatePoint(&TransformedLine.Origin,&tempv);
			RevRotateVector(&TransformedLine.Direction,&tempv);
			break;
		    case TRANSFORM_NONE:
		    default:
			break;
		}
	    }
	  }

	  switch(Object->ShapeType) {
	    case SHAPE_PLANE:
		t=Intersect_Plane(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_SPHERE:
		t=Intersect_Sphere(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_ELLIPSOID:
		t=Intersect_Ellipsoid(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_TRIANGLE:
		t=Intersect_Triangle(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_TRIANGLELIST:
		t=Intersect_TriangleList(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_BOX:
		t=Intersect_Box(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_DISC:
		t=Intersect_Disc(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_CYLINDER:
		t=Intersect_Cylinder(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_CONE:
		t=Intersect_Cone(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_TORUS:
		t=Intersect_Torus(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_PEAK:
		t=Intersect_Peak(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_DEOFLASKA:
		t=Intersect_Deoflaska(Intrsct, &TransformedLine, Object);
		break;
	    case SHAPE_CSG:
		t=Intersect_CSG(Intrsct, &TransformedLine, Object);
		break;
	    default:
		t=0.0;
		break;
	  }

	  if(t>EPSILON) {
	    if(Object->Transform.NumTransforms>0) {	/* Transform intersection to its real place/orientation in space */
		for(i=0L;i<Object->Transform.NumTransforms;i++) {
		    tempv=Object->Transform.Entry[i].Values;
		    switch(Object->Transform.Entry[i].Type) {
			case TRANSFORM_SCALE:
				a=tempv.x; b=tempv.y; c=tempv.z;
				Intrsct->IPoint.x*=a; Intrsct->IPoint.y*=b; Intrsct->IPoint.z*=c;
				Intrsct->Normal.x*=1.0/a; Intrsct->Normal.y*=1.0/b; Intrsct->Normal.z*=1.0/c;
				break;
			case TRANSFORM_MOVE:
				Intrsct->IPoint.x+=tempv.x;
				Intrsct->IPoint.y+=tempv.y;
				Intrsct->IPoint.z+=tempv.z;
				break;
			case TRANSFORM_ROTATE:
				RotatePoint(&Intrsct->IPoint,&tempv);
				RotateVector(&Intrsct->Normal,&tempv);
				break;
			case TRANSFORM_NONE:
			default:
				break;
		    }
		}
	    }
	    if(Object->Invert==INVERT_TRUE) {
		NegVector(&Intrsct->Normal,&Intrsct->Normal);	/* "Reverse" surface */
	    }
	  }
	}

	return(t);
}



/**********************************************************************
 *
 *	These are the actual intersection calculation routines
 *
 **********************************************************************/


double Intersect_Plane(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	t,pnx,pny,pnz,a,tmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	PLANE	*Plane;

	PlaneIntersectionAttempts++;
	Plane=(PLANE *)Obj->Shape;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Plane->Normal.x; pny=Plane->Normal.y; pnz=Plane->Normal.z; a=Plane->a;

	tmp=pnx*lvx+pny*lvy+pnz*lvz;
	if(tmp!=0.0) {
	    t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
	    if(t>EPSILON) {
		PlaneIntersections++;
		Intrsct->IPoint.x=lx0+lvx*t; Intrsct->IPoint.y=ly0+lvy*t; Intrsct->IPoint.z=lz0+lvz*t;
		Intrsct->Normal.x=pnx; Intrsct->Normal.y=pny; Intrsct->Normal.z=pnz;
		Intrsct->Texture=&Obj->Texture;
		Intrsct->Shadow=Obj->Shadow;
	    }
	    else t=0.0;
	}
	else t=0.0;

	return(t);
}



double Intersect_Sphere(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	double	vx,vy,vz,lx0,ly0,lz0,sx0,sy0,sz0,r;

	register double	s,sroot,vxdx,vydy,vzdz,vxsqr,vysqr,vzsqr;
	double	t;
	double	dx,dy,dz;
	double	dxsqr,dysqr,dzsqr;
	double	vsum,t1,t2;
	SPHERE	*Sphere;

	SphereIntersectionAttempts++;
	Sphere=(SPHERE *)Obj->Shape;

	vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	sx0=Sphere->Centre.x; sy0=Sphere->Centre.y; sz0=Sphere->Centre.z;
	r=Sphere->r;

						/* These are used more than once, */
	dx=lx0-sx0; dy=ly0-sy0; dz=lz0-sz0;	/* thus much time is saved by not */
	vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;	/* performing the same operation  */
	dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;	/* twice (or even more).          */
	vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;

	vsum=(vxsqr+vysqr+vzsqr);

	sroot=2*(vxdx*(vydy+vzdz)+vydy*vzdz)+r*r*vsum;
	sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);

	t=0.0;
	if(sroot>=0) {
		sroot=sqrt(sroot);

		s=-(vxdx+vydy+vzdz);
		t1=(s+sroot)/vsum;		/* Note: vsum!=0, always */
		t2=(s-sroot)/vsum;		/* Also: t2<=t1, as sroot>=0 */

		if(t2>EPSILON) t=t2;
		else if(t1>EPSILON) t=t1;
		if(t>EPSILON) {
			SphereIntersections++;
			Intrsct->IPoint.x=lx0+vx*t; Intrsct->IPoint.y=ly0+vy*t; Intrsct->IPoint.z=lz0+vz*t;
			Intrsct->Normal.x=Intrsct->IPoint.x-sx0;
			Intrsct->Normal.y=Intrsct->IPoint.y-sy0;
			Intrsct->Normal.z=Intrsct->IPoint.z-sz0;
			Intrsct->Texture=&Obj->Texture;
			Intrsct->Shadow=Obj->Shadow;
		}
	}

	return(t);
}



double Intersect_Ellipsoid(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	double	vx,vy,vz,lx0,ly0,lz0,ex0,ey0,ez0,a,b,c;

	register double	sroot,asqr,bsqr,csqr,vxdx,vydy,vzdz;
	double	t;
	double	dx,dy,dz;
	double	vxsqr,vysqr,vzsqr,dxsqr,dysqr,dzsqr;
	double	r,s,t1,t2;
	ELLIPSOID *Ellipsoid;

	EllipsoidIntersectionAttempts++;
	Ellipsoid=(ELLIPSOID *)Obj->Shape;

	vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	ex0=Ellipsoid->Centre.x; ey0=Ellipsoid->Centre.y; ez0=Ellipsoid->Centre.z;
	a=Ellipsoid->Radius.x; b=Ellipsoid->Radius.y; c=Ellipsoid->Radius.z;

	asqr=a*a; bsqr=b*b; csqr=c*c;		/* These are used more than once, */
	dx=lx0-ex0; dy=ly0-ey0; dz=lz0-ez0;	/* thus much time is saved by not */
	vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz;	/* performing the same operation  */
	dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;	/* twice (or even more).          */
	vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;

	s=-((vxdx/asqr)+(vydy/bsqr)+(vzdz/csqr));

	sroot=csqr*(vxsqr*bsqr+vysqr*asqr)+vzsqr*asqr*bsqr;
	sroot+=2*(vxdx*(csqr*vydy+bsqr*vzdz)+asqr*vydy*vzdz);
	sroot-=vxsqr*(csqr*dysqr+bsqr*dzsqr);
	sroot-=vysqr*(csqr*dxsqr+asqr*dzsqr);
	sroot-=vzsqr*(bsqr*dxsqr+asqr*dysqr);

	t=0.0;
	if(sroot>=0) {
		sroot=sqrt(sroot)/(a*b*c);

		r=((vxsqr/asqr)+(vysqr/bsqr)+(vzsqr/csqr));

		t1=(s+sroot)/r;
		t2=(s-sroot)/r;

		if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
			t=t1;
			EllipsoidIntersections++;
		}
		else if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
			t=t2;
			EllipsoidIntersections++;
		}
		if(t>EPSILON) {
			Intrsct->IPoint.x=lx0+vx*t; Intrsct->IPoint.y=ly0+vy*t; Intrsct->IPoint.z=lz0+vz*t;
			Intrsct->Normal.x=(Intrsct->IPoint.x-ex0)/(a*a);
			Intrsct->Normal.y=(Intrsct->IPoint.y-ey0)/(b*b);
			Intrsct->Normal.z=(Intrsct->IPoint.z-ez0)/(c*c);
			Intrsct->Texture=&Obj->Texture;
			Intrsct->Shadow=Obj->Shadow;
		}
	}

	return(t);
}



double Intersect_Triangle(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	double	t,ttmp,pnx,pny,pnz;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	double	ui,vi,u[4],v[4];
	register double un,vn,um,vm;
	short	sign,count,n,m;
	TRIANGLE *Triangle;

	TriangleIntersectionAttempts++;
	Triangle=(TRIANGLE *)Obj->Shape;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Triangle->Normal.x; pny=Triangle->Normal.y; pnz=Triangle->Normal.z;

	t=0.0;
	ttmp=pnx*lvx+pny*lvy+pnz*lvz;
	if(ttmp!=0.0) {
	    ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-Triangle->Offset)/ttmp;
	    if(ttmp>EPSILON) {
		switch(Triangle->Projection) {	/* Transform to 2D coordinates */
		    case PROJ_YZ:
			ui=ly0+ttmp*lvy;
			vi=lz0+ttmp*lvz;
			break;
		    case PROJ_XZ:
			ui=lx0+ttmp*lvx;
			vi=lz0+ttmp*lvz;
			break;
		    case PROJ_XY:
			ui=lx0+ttmp*lvx;
			vi=ly0+ttmp*lvy;
			break;
		}
		u[0]=Triangle->Corners[0].u-ui;
		u[1]=Triangle->Corners[1].u-ui;
		u[2]=Triangle->Corners[2].u-ui;
		v[0]=Triangle->Corners[0].v-vi;
		v[1]=Triangle->Corners[1].v-vi;
		v[2]=Triangle->Corners[2].v-vi;

		u[3]=u[0]; v[3]=v[0];
		count=0;			/* Determine amount of sides crossed */
		for(n=0;n<3;n++) {
		    m=n+1;
		    vn=v[n]; vm=v[m];
		    if(vn>=0.0) sign=1;  else sign=-1;
		    if(vm>=0.0) sign+=1; else sign-=1;
		    if(sign==0) {
			un=u[n]; um=u[m];
			if((un>=0.0)&&(um>=0.0)) count+=1;
			else if((un>=0.0)||(um>=0.0)) {
			    if((um-vm*(um-un)/(vm-vn))>0.0) count+=1;
			}
		    }
		}
		if((count&0x01)==1) {	/* If odd amount of sides were crossed... */
		    TriangleIntersections++;
		    t=ttmp;
		    Intrsct->Normal.x=pnx;
		    Intrsct->Normal.y=pny;
		    Intrsct->Normal.z=pnz;
		    Intrsct->IPoint.x=lx0+lvx*t;
		    Intrsct->IPoint.y=ly0+lvy*t;
		    Intrsct->IPoint.z=lz0+lvz*t;
		    Intrsct->Texture=&Obj->Texture;
		    Intrsct->Shadow=Obj->Shadow;
		}
	    }
	}

	return(t);
}



double Intersect_TriangleList(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	double	t,ttmp,pnx,pny,pnz;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	double	ui,vi,u[4],v[4];
	register double un,vn,um,vm;
	short	sign,count,n,m;
	TRIANGLELIST *Triangle;

	Triangle=(TRIANGLELIST *)Obj->Shape;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;

	t=0.0;

	while(Triangle!=NULL) {
	    TriangleIntersectionAttempts++;
	    pnx=Triangle->Normal.x; pny=Triangle->Normal.y; pnz=Triangle->Normal.z;
	    ttmp=pnx*lvx+pny*lvy+pnz*lvz;
	    if(ttmp!=0.0) {
		ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-Triangle->Offset)/ttmp;
		if(ttmp>EPSILON) {
		    switch(Triangle->Projection) {	/* Transform to 2D coordinates */
			case PROJ_YZ:
			    ui=ly0+ttmp*lvy;
			    vi=lz0+ttmp*lvz;
			    break;
			case PROJ_XZ:
			    ui=lx0+ttmp*lvx;
			    vi=lz0+ttmp*lvz;
			    break;
			case PROJ_XY:
			    ui=lx0+ttmp*lvx;
			    vi=ly0+ttmp*lvy;
			    break;
		    }
		    u[0]=Triangle->Corners[0].u-ui;
		    u[1]=Triangle->Corners[1].u-ui;
		    u[2]=Triangle->Corners[2].u-ui;
		    v[0]=Triangle->Corners[0].v-vi;
		    v[1]=Triangle->Corners[1].v-vi;
		    v[2]=Triangle->Corners[2].v-vi;

		    u[3]=u[0]; v[3]=v[0];
		    count=0;			/* Determine amount of sides crossed */
		    for(n=0;n<3;n++) {
			m=n+1;
			vn=v[n]; vm=v[m];
			if(vn>=0.0) sign=1;  else sign=-1;
			if(vm>=0.0) sign+=1; else sign-=1;
			if(sign==0) {
			    un=u[n]; um=u[m];
			    if((un>=0.0)&&(um>=0.0)) count+=1;
			    else if((un>=0.0)||(um>=0.0)) {
				if((um-vm*(um-un)/(vm-vn))>0.0) count+=1;
			    }
			}
		    }
		}

		if((count&0x01)==1) {	/* If odd amount of sides were crossed... */
		    if((ttmp>EPSILON)&&((t<EPSILON)||(ttmp<t))) {
			TriangleIntersections++;
			t=ttmp;
			Intrsct->Normal.x=pnx;
			Intrsct->Normal.y=pny;
			Intrsct->Normal.z=pnz;
			Intrsct->IPoint.x=lx0+lvx*t;
			Intrsct->IPoint.y=ly0+lvy*t;
			Intrsct->IPoint.z=lz0+lvz*t;
		    }
		}
	    }

	    Triangle=Triangle->NextTriangle;
	}

	if(t>EPSILON) {
		Intrsct->Texture=&Obj->Texture;
		Intrsct->Shadow=Obj->Shadow;
	}

	return(t);
}


double Intersect_Box(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	t,told,ttmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	register double	ipx,ipy,ipz;
	POINT	ip;
	VECTOR	norm;
	BOX	*Box;

	BoxIntersectionAttempts++;
	Box=(BOX *)Obj->Shape;

	told=-1.0;
	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;

	if(lvz!=0.0) {
	    ttmp=(Box->Corners[0].z-lz0)/lvz;	/* xy-plane, z-min */
	    if(ttmp>EPSILON) {
		ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)) {
			told=ttmp; ip.x=ipx; ip.y=ipy; ip.z=Box->Corners[0].z; norm.x=norm.y=0.0; norm.z=-1.0;
		}
	    }
	    ttmp=(Box->Corners[1].z-lz0)/lvz;	/* xy-plane, z-max */
	    if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
		ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)) {
			told=ttmp; ip.x=ipx; ip.y=ipy; ip.z=Box->Corners[1].z; norm.x=norm.y=0.0; norm.z=1.0;
		}
	    }
	}
	if(lvy!=0.0) {
	    ttmp=(Box->Corners[0].y-ly0)/lvy;	/* xz-plane, y-min */
	    if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
		ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
			told=ttmp; ip.x=ipx; ip.z=ipz; ip.y=Box->Corners[0].y; norm.x=norm.z=0.0; norm.y=-1.0;
		}
	    }
	    ttmp=(Box->Corners[1].y-ly0)/lvy;	/* xz-plane, y-max */
	    if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
		ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
			told=ttmp; ip.x=ipx; ip.z=ipz; ip.y=Box->Corners[1].y; norm.x=norm.z=0.0; norm.y=1.0;
		}
	    }
	}
	if(lvy!=0.0) {
	    ttmp=(Box->Corners[0].x-lx0)/lvx;	/* yz-plane, x-min */
	    if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
		ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
		if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
			told=ttmp; ip.y=ipy; ip.z=ipz; ip.x=Box->Corners[0].x; norm.y=norm.z=0.0; norm.x=-1.0;
		}
	    }
	    ttmp=(Box->Corners[1].x-lx0)/lvx;	/* yz-plane, x-max */
	    if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
		ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
		if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z)) {
			told=ttmp; ip.y=ipy; ip.z=ipz; ip.x=Box->Corners[1].x; norm.y=norm.z=0.0; norm.x=1.0;
		}
	    }
	}

	if(told>EPSILON) {
		t=told;
		BoxIntersections++;
		Intrsct->IPoint=ip;
		Intrsct->Normal=norm;
		Intrsct->Texture=&Obj->Texture;
		Intrsct->Shadow=Obj->Shadow;
	}
	else t=0.0;

	return(t);
}


double Intersect_Disc(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	t;

	DiscIntersectionAttempts++;

	t=Intersect_DiscShape(Intrsct, Line, (DISC *)Obj->Shape);
	if(t>EPSILON) {
		Intrsct->Texture=&Obj->Texture;
		Intrsct->Shadow=Obj->Shadow;
	}

	return(t);
}


double Intersect_DiscShape(INTERSECTION *Intrsct, LINE *Line, DISC *Disc)
{
	register double	t,pnx,pny,pnz,a,tmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0,dx,dy,dz;
	POINT	ip;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
	pnx=Disc->Plane.Normal.x; pny=Disc->Plane.Normal.y; pnz=Disc->Plane.Normal.z; a=Disc->Plane.a;

	t=0.0;
	tmp=pnx*lvx+pny*lvy+pnz*lvz;
	if(tmp!=0.0) {
		tmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
		if(tmp>EPSILON) {
			ip.x=lx0+lvx*tmp; ip.y=ly0+lvy*tmp; ip.z=lz0+lvz*tmp;
			dx=ip.x-Disc->Centre.x;
			dy=ip.y-Disc->Centre.y;
			dz=ip.z-Disc->Centre.z;
			if((dx*dx+dy*dy+dz*dz)<(Disc->r*Disc->r)) {
				t=tmp;
				DiscIntersections++;
				Intrsct->IPoint=ip;
				Intrsct->Normal.x=pnx; Intrsct->Normal.y=pny; Intrsct->Normal.z=pnz; 
			}
		}
	}

	return(t);
}


double Intersect_Cylinder(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{ 
	register double lvx,lvy,l0x,l0y,lvx2plvy2;
	double	rsqr,a,b,c,t,t1,t2,z1,z2;
	int	sidehit;
	CYLINDER *Cylinder;

	CylinderIntersectionAttempts++;
	Cylinder=(CYLINDER *)Obj->Shape;

	lvx=Line->Direction.x; lvy=Line->Direction.y;
	l0x=Line->Origin.x; l0y=Line->Origin.y;
	rsqr=Cylinder->r; rsqr*=rsqr;
	lvx2plvy2=1.0/(lvx*lvx+lvy*lvy);

	a=(rsqr-l0x*l0x-l0y*l0y)*lvx2plvy2;
	b=(lvx*l0x+lvy*l0y)*lvx2plvy2;
	c=a+b*b;
	t=0.0;
	if(c>=0.0) {
		c=sqrt(c);
		t1=-b-c;	/*  t1 <= t2  */
		t2=-b+c;
		z1=z2=Line->Origin.z;
		z1+=Line->Direction.z*t1;
		z2+=Line->Direction.z*t2;
		a=Cylinder->h;		/* a = height of cylinder */
		sidehit=FALSE;
		if(t1>EPSILON) {
		    if((z1>=0.0)&&(z1<=a)) {
			t=t1;
			Intrsct->IPoint.z=z1;
			sidehit=TRUE;
		    }
		}
		else if(t2>EPSILON) {
		    if((z2>=0.0)&&(z2<=a)) {
			t=t2;
			Intrsct->IPoint.z=z2;
			sidehit=TRUE;
		    }
		}
		b=Line->Origin.z; c=1.0/Line->Direction.z;
		t1=(0.0-b)*c;			/* t1 <=> bottom hit */
		t2=(a-b)*c;			/* t2 <=> top hit */
		if(t1>EPSILON) {
		    b=l0x+lvx*t1; c=l0y+lvy*t1;
		    if((b*b+c*c)>rsqr) t1=0.0;		/* Outside of the cyliunder? */
		}
		if(t2>EPSILON) {
		    if((t1>EPSILON)&&(t1<t2)) t2=0.0;
		    else {
			b=l0x+lvx*t2; c=l0y+lvy*t2;
			if((b*b+c*c)>rsqr) t2=0.0;
			else t1=0.0;
		    }
		}
		if((t1>EPSILON)&&((t1<t)||(t<EPSILON))) {
		    t=t1;
		    CreatePoint(&Intrsct->IPoint,l0x+lvx*t1,l0y+lvy*t1,0.0);
		    CreateVector(&Intrsct->Normal,0.0,0.0,-1.0);
		    sidehit=FALSE;
		}
		else if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) {
		    t=t2;
		    CreatePoint(&Intrsct->IPoint,b,c,a);
		    CreateVector(&Intrsct->Normal,0.0,0.0,1.0);
		    sidehit=FALSE;
		}
		if(sidehit==TRUE) {
			Intrsct->Normal.x=Intrsct->IPoint.x=l0x+lvx*t;
			Intrsct->Normal.y=Intrsct->IPoint.y=l0y+lvy*t;
			Intrsct->Normal.z=0.0;
		}
		if(t>=EPSILON) {
			CylinderIntersections++;
			Intrsct->Texture=&Obj->Texture;
			Intrsct->Shadow=Obj->Shadow;
		}
	}

	return(t);
}


double Intersect_Cone(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{ 
	double	lvx,lvy,lvz,l0x,l0y,l0z;
	double	r1,r2,h,A,B,C,a,b,c,t,t1,t2,x,y,z;
	int	sidehit;
	CONE	*Cone;

	ConeIntersectionAttempts++;
	Cone=(CONE *)Obj->Shape;

	t = 0.0;

	r1 = Cone->r1; r2 = Cone->r2; h = Cone->h;
	lvx = Line->Direction.x; lvy = Line->Direction.y; lvz = Line->Direction.z;
	l0x = Line->Origin.x; l0y = Line->Origin.y; l0z = Line->Origin.z;

	c = (r2-r1)/h;
	A = (lvx*lvx + lvy*lvy - c*c*lvz*lvz);
	if(A!=0.0) {
	    A = 1.0 / A;
	    B = lvx*l0x + lvy*l0y - lvz*c*(l0z*c + r1);
	    C = l0x*l0x + l0y*l0y - l0z*c*(l0z*c + 2.0*r1) - r1*r1;

	    a = A*B;
	    b = a*a - A*C;
	    if(b>=0) {
		b = sqrt(b);		/* b > 0 */
		t1 = -a - b;		/* t1 < t2 */
		t2 = -a + b;

		sidehit = FALSE;
	/* Intersection with sides */
		if(t2 > EPSILON) {
		    z = l0z + t2*lvz;
		    if((z>0.0)&&(z<h)) {
			t = t2; sidehit = TRUE;
		    }
		    if(t1 > EPSILON) {
			z = l0z + t1*lvz;
			if((z>0.0)&&(z<h)) {
			    t = t1; sidehit = TRUE;
			}
		    }
		}
	/* Intersection with top/bottom */
		if(lvz!=0.0) {
		    t1 = -l0z/lvz;  t2 = (h-l0z)/lvz;
		    if((t1>EPSILON)&&((t<EPSILON)||(t1<t))) {
			x = l0x+lvx*t1; y = l0y+lvy*t1;
			if((x*x + y*y) <= (r1*r1)) {
			    t = t1; sidehit = FALSE;
			    Intrsct->IPoint.x = x; Intrsct->IPoint.y = y;
			    Intrsct->Normal.x = Intrsct->Normal.y = 0.0;
			    Intrsct->Normal.z = -1.0;
			}
		    }
		    if((t2>EPSILON)&&((t<EPSILON)||(t2<t))) {
			x = l0x+lvx*t2; y = l0y+lvy*t2;
			if((x*x + y*y) <= (r2*r2)) {
			    t = t2; sidehit = FALSE;
			    Intrsct->IPoint.x = x; Intrsct->IPoint.y = y;
			    Intrsct->Normal.x = Intrsct->Normal.y = 0.0;
			    Intrsct->Normal.z = 1.0;
			}
		    }
		}
		if(t>EPSILON) {
		    ConeIntersections++;
		    Intrsct->IPoint.z = l0z+lvz*t;
		    if(sidehit == TRUE) {
			x = Intrsct->IPoint.x = Intrsct->Normal.x = l0x+lvx*t;
			y = Intrsct->IPoint.y = Intrsct->Normal.y = l0y+lvy*t;
			Intrsct->Normal.z = -sqrt(x*x + y*y)*c;
		    }
		    Intrsct->Texture = &Obj->Texture;
		    Intrsct->Shadow = Obj->Shadow;
		}
	    }
	}

	return(t);
}


double Torus_Function(double t, double r0, double r1, LINE *Line)
{
	register double	a,b,c,f;

	a=Line->Origin.x+Line->Direction.x*t;
	b=Line->Origin.y+Line->Direction.y*t;
	c=Line->Origin.z+Line->Direction.z*t;

	f = sqrt(a*a + b*b) - r0;
	f = f*f + c*c - r1*r1;

	return(f);
}


double Torus_Derivata(double t, double r0, LINE *Line)
{
	register double	a,b,c,df;

	a=Line->Origin.x+Line->Direction.x*t;
	b=Line->Origin.y+Line->Direction.y*t;
	c=Line->Origin.z+Line->Direction.z*t;

	df = (1.0 - r0/sqrt(a*a + b*b))*(a + b);
	df *= (Line->Direction.x + Line->Direction.y);
	df += c*Line->Direction.z;
	df += df;

	return(df);
}

double Intersect_Torus(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	r0,r1,tn,fn,tb,fb,diff;
	double	tnp1,fnp1;
	double	dx,dy,dz,x0,y0,z0;
	double	cylr,tcyl1o,tcyl2o,tcyl1i,tcyl2i,tpln1,tpln2,z1,z2,pr1,pr2;
	double	dx2pdy2,a,b,c,ts[6],t;
	int	i,tamount,sortchange,FoundHit;
	TORUS	*Torus;

	TorusIntersectionAttempts++;
	Torus=(TORUS *)Obj->Shape;

	dx=Line->Direction.x; dy=Line->Direction.y; dz=Line->Direction.z;
	x0=Line->Origin.x; y0=Line->Origin.y; z0=Line->Origin.z;
	r0=Torus->r0; r1=Torus->r1;
	dx2pdy2=1.0/(dx*dx+dy*dy);

	cylr=(r0+r1);		/* Check for outer cylinder-bounding */
	a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
	b=(dx*x0+dy*y0)*dx2pdy2;
	c=a+b*b;

	t=0.0;
	if(c>=0.0) {		/* If found cylinder, try the rest */
	    c=sqrt(c);
	    tcyl1o=-b-c;	/*  tcyl1o < tcyl2o */
	    tcyl2o=-b+c;
	    if(tcyl2o>EPSILON) {  /* Don't do anything if BOTH hits behind */
		z1=(z0+dz*tcyl1o);
		z2=(z0+dz*tcyl2o);
		if(!(((z1>r1)&&(z2>r1))||((z1<-r1)&&(z2<-r1)))) {
		    if(fabs(z1)>r1) tcyl1o=0.0;
		    if(fabs(z2)>r1) tcyl2o=0.0;
		    cylr=(r0-r1);	  /* Check for inner cylinder */
		    a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
		    b=(dx*x0+dy*y0)*dx2pdy2;
		    c=a+b*b;
		    if(c>=0.0) {
			c=sqrt(c);
			tcyl1i=-b-c;	/*  tcyl1i < tcyl2i */
			tcyl2i=-b+c;
			z1=(z0+dz*tcyl1i);
			z2=(z0+dz*tcyl2i);
			if(fabs(z1)>r1) tcyl1i=0.0;
			if(fabs(z2)>r1) tcyl2i=0.0;
		    }
		    else {
			tcyl1i=tcyl2i=0.0;
		    }
		    tpln1=(r1-z0)/dz;	/* Check for plane intersection */
		    tpln2=(-r1-z0)/dz;
		    a=x0+dx*tpln1; b=y0+dy*tpln1; pr1=sqrt(a*a+b*b);
		    a=x0+dx*tpln2; b=y0+dy*tpln2; pr2=sqrt(a*a+b*b);
		    if((pr1>(r0+r1))||(pr1<(r0-r1))) tpln1=0.0;
		    if((pr2>(r0+r1))||(pr2<(r0-r1))) tpln2=0.0;
		    if( (tcyl1o>EPSILON)||(tcyl2o>EPSILON)||
			(tcyl1i>EPSILON)||(tcyl2i>EPSILON)||
			(tpln1>EPSILON)||(tpln2>EPSILON) ) {
			tamount=0;
			if(fabs(tcyl1o)>EPSILON) ts[tamount++]=tcyl1o;
			if(fabs(tcyl2o)>EPSILON) ts[tamount++]=tcyl2o;
			if(fabs(tcyl1i)>EPSILON) ts[tamount++]=tcyl1i;
			if(fabs(tcyl2i)>EPSILON) ts[tamount++]=tcyl2i;
			if(fabs(tpln1)>EPSILON)  ts[tamount++]=tpln1;
			if(fabs(tpln2)>EPSILON)  ts[tamount++]=tpln2;
			if(((tamount&1)!=0)||(tamount>4)) {
			    tamount=-1;
			}
			else {
			    sortchange=1;
			    while(sortchange>0) {  /* Sort bound-hits */
				sortchange=0;	   /* (closest last)  */
				for(i=0;i<tamount-1;i++) {
				    if(ts[i+1]>ts[i]) {
					a=ts[i]; ts[i]=ts[i+1]; ts[i+1]=a;
					sortchange++;
				    }
				}
			    }
			    a=ts[1];
			    if(tamount>2) {
				 b=ts[2]; c=ts[3];
			    }
			    ts[1]=(ts[0]+a)*0.5; ts[2]=a;
			    tamount++;
			    if(tamount>3) {
				ts[3]=b; ts[4]=(b+c)*0.5; ts[5]=c;
				tamount++;
			    }
			}
			FoundHit=FALSE;
			while((tamount>=2)&&(FoundHit==FALSE)) {
			    tb=ts[tamount-2];	/* tb is further away than... */
			    tn=ts[tamount-1];	/* ...tn  */
			    fb=Torus_Function(tb,r0,r1,Line);
			    fn=Torus_Function(tn,r0,r1,Line);
			    if(((fn<0.0)&&(fb>0.0))||((fn>0.0)&&(fb<0.0))) {
				if((tb>EPSILON)||(tn>EPSILON)) {
				    i=0; diff=1.0;	/* Trapets-method */
				    while((diff>NUMEPSILON)&&(i<10)) {
					tnp1=tn-fn*(tb-tn)/(fb-fn);
					diff=fabs(tnp1-tn);
					if(diff>NUMEPSILON) {
					    fnp1=Torus_Function(tnp1,r0,r1,Line);
					    if(((fnp1<0.0)&&(fb<0.0))||((fnp1>0.0)&&(fb>0.0))) {
						fb=fn; tb=tn;	/* Switch direction */
					    }
					    fn=fnp1;
					    i++;
					}
					tn=tnp1;
				    }
				    if(tn>NUMEPSILON) {
					FoundHit=TRUE;
				    }
				}
			    }
			    else {			/* Newton-Raphson */
				if((tamount==2)||(tamount==5)) {
				    tn=tb; fn=fb;
				}
				i=0; diff=1.0;
				while((diff>NUMEPSILON)&&(i<6)) {
				    diff=fn/Torus_Derivata(tn,r0,Line);
				    tn-=diff;
				    diff=fabs(diff);
				    if(diff>NUMEPSILON) {
					fn=Torus_Function(tn,r0,r1,Line);
					i++;
				    }
				}
				if((diff<=NUMEPSILON)&&(tn>NUMEPSILON)) {
				    FoundHit=TRUE;
				}
			    }
			    if(FoundHit==FALSE) {
				tamount--;
				if(tamount==4) tamount--;
			    }
			}
			if(FoundHit==TRUE) {
			    TorusIntersections++;
			    t=tn;
			    Intrsct->Texture=&Obj->Texture;
			    Intrsct->Shadow=Obj->Shadow;
			    Intrsct->IPoint.x=x0+dx*t;
			    Intrsct->IPoint.y=y0+dy*t;
			    Intrsct->IPoint.z=z0+dz*t;
			    a=Intrsct->IPoint.x; b=Intrsct->IPoint.y;
			    c=1.0-r0/sqrt(a*a+b*b);
			    Intrsct->Normal.x=a*c;
			    Intrsct->Normal.y=b*c;
			    Intrsct->Normal.z=Intrsct->IPoint.z;
			}
		    }
		}
	    }
	}

	return(t);
}


double Peak_Function(double t, LINE *Line)
{
	register double x,y,z,f;

	x = Line->Origin.x + Line->Direction.x*t;
	y = Line->Origin.y + Line->Direction.y*t;
	z = Line->Origin.z + Line->Direction.z*t;

	f = x*x + y*y;
	if(z>0.0) f -= 1.0/(z*z);
	else f += 1.0/(z*z);

	return( f );
}

double Peak_Derivata(double t, LINE *Line)
{
	register double x,y,z,f;

	x = Line->Origin.x + Line->Direction.x*t;
	y = Line->Origin.y + Line->Direction.y*t;
	z = Line->Origin.z + Line->Direction.z*t;

	f = x*Line->Direction.x + y*Line->Direction.y;
	if(z>0.0) f += Line->Direction.z/(z*z*z);
	else f -= Line->Direction.z/(z*z*z);

	return( 2.0*f );
}

double Intersect_Peak(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	double		lvx,lvy,lvz,lx0,ly0,lz0;
	double		t,t_a,t_b,f1,f2,tnp1,fnp1,fprim;
	register double	A,B,C,D,a,b;
	int		i;

	PeakIntersectionAttempts++;

	lvx = Line->Direction.x; lvy = Line->Direction.y; lvz = Line->Direction.z;
	lx0 = Line->Origin.x; ly0 = Line->Origin.y; lz0 = Line->Origin.z;

	a = 1.0/lvx; b = 1.0/lvy;
	t_a = -(ly0*a+lx0*b)/(lvy*a+lvx*b);	/* Closest point to z-axis */
	t_b = -lz0/lvz;				/* Intersection with z=0 */


  /* Use Newton-Raphson to find hit with "floor" */
	tnp1 = t_b; b = 1.0; i=0;
	tnp1 = 0.0;
	while((fabs(b)>NUMEPSILON)&&((i++)<6)) {
	    b = Peak_Function(tnp1,Line) / Peak_Derivata(tnp1,Line);
	    tnp1 -= b;
	}
	t_b = tnp1;				/* t_b = intersection with "floor" */


 /*	a = t_a;
	f1 = D + C*a;
	a *= t_a; f1 += B*a;
	a *= t_a; f1 += A*a; */			/* f1 = f(t_a) */

 /*	if(f1>0.0) !utanför! */

	t = t_b;

	if(t>(NUMEPSILON*4.0)) {
	    PeakIntersections++;
	    a = lx0 + lvx*t; b = ly0 + lvy*t;
	    Intrsct->IPoint.x=a; Intrsct->Normal.x=a;
	    Intrsct->IPoint.y=b; Intrsct->Normal.y=b;
	    a = lz0 + lvz*t;
	    Intrsct->IPoint.z=a; Intrsct->Normal.y=1.0/(a*a*a);
	    NormalizeVector(&Intrsct->Normal,&Intrsct->Normal); /* Can be _HUGE_ and _SMALL_! */
	    Intrsct->Texture = &Obj->Texture;
	    Intrsct->Shadow = Obj->Shadow;
	}
	else t = 0.0;

	return(t);
}


double Intersect_Deoflaska(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	t,r0,dr,h,tn,tnp1,t1,t2,diff;
	double	dx,dy,dz,x0,y0,z0,unitdz;
	double	tcyl1,tcyl2,cylr,dx2pdy2,a,b,c,z1,z2;
	int	i, MultiplePeriods, LastCheck, FoundHit;
	DEOFLASKA	*Deo;

	DeoflaskaIntersectionAttempts++;
	Deo=(DEOFLASKA *)Obj->Shape;

	dx=Line->Direction.x; dy=Line->Direction.y; dz=Line->Direction.z;
	x0=Line->Origin.x; y0=Line->Origin.y; z0=Line->Origin.z;
	r0=Deo->r0; dr=Deo->dr; h=Deo->h;

	cylr=(r0+dr)*1.03;	/* Check for cylinder-bounding */
	dx2pdy2=1.0/(dx*dx+dy*dy);

	a=(cylr*cylr-x0*x0-y0*y0)*dx2pdy2;
	b=(dx*x0+dy*y0)*dx2pdy2;
	c=a+b*b;

	if(c>=0.0) {		/* If found cylinder, try numerical root-finding */
	    c=sqrt(c);
	    tcyl1=-b-c;		/*  tcyl1 < tcyl2 */
	    tcyl2=-b+c;
	    if(tcyl2>EPSILON) {	/* Don't do anything if BOTH hits behind */
		FoundHit=FALSE;
		LastCheck=FALSE;
		if(((dz*dz)/(dx*dx+dy*dy))>0.25) MultiplePeriods=TRUE;
		else MultiplePeriods=FALSE;
		t1=t2=tcyl1;
		unitdz=fabs(dz)/sqrt(dx*dx+dy*dy+dz*dz);
		while((FoundHit==FALSE)&&(LastCheck==FALSE)) {
		    t1=t2;		/* Determine interval t1 <= t <= t2 containing max one root */
		    if(MultiplePeriods==TRUE) {
			z1=(z0+dz*t1);
			if(dz>0.0) {
			    z2=(floor((z1+PID2)/PI)*PI)+PID2;
/*			    z2=floor((z1+PI)/PI)*PI; */
			    if((z2-z1)<EPSILON) z2+=PI;
			} else {
			    z2=(floor((z1-PID2)/PI)*PI)+PID2;
/*			    z2=floor(z1/PI)*PI; */
			    if((z1-z2)<EPSILON) z2-=PI;
			}
			t2=(z2-z0)/dz;
			if(t2>=tcyl2) {
			    t2=tcyl2;
			    LastCheck=TRUE;
			}
			tn=(t1+t2)*0.5;	/* Start in the middle of the interval */
		    }
		    else {
			if(t2<=(tcyl1+NUMEPSILON)) {
			    t1=tcyl1;
			    t2=(tcyl1+tcyl2)*0.5;
			    tn=tcyl1;
			}
			else {
			    t1=(tcyl1+tcyl2)*0.5;
			    t2=tcyl2;
			    tn=tcyl2;
			    LastCheck=TRUE;
			}
		    }

		    if(t2>EPSILON) {
			i=0; diff=1.0;
			while((diff>NUMEPSILON)&&(i<6)) {
			    a=(x0+dx*tn); b=(y0+dy*tn); c=(r0-dr*sin(z0+dz*tn));
			    diff=(a*a+b*b-c*c)/(2*a*dx+2*b*dy+2*c*dr*cos(z0+dz*tn)*dz);
			    tnp1=tn-diff;
			    diff=fabs(diff);
			    tn=tnp1;
			    i++;
			}
			if((diff<NUMEPSILON)&&(tn>(4.0*NUMEPSILON))&&(tn>=t1)&&(tn<=t2)) FoundHit=TRUE;
			else FoundHit=FALSE;
		    }
		}
		if((FoundHit==TRUE)&&(tn>(4.0*NUMEPSILON))&&(tn>=tcyl1)&&(tn<=tcyl2)) {
			DeoflaskaIntersections++;
			t=tn;
			Intrsct->Texture=&Obj->Texture;
			Intrsct->Shadow=Obj->Shadow;
			Intrsct->IPoint.x=Intrsct->Normal.x=x0+dx*t;
			Intrsct->IPoint.y=Intrsct->Normal.y=y0+dy*t;
			Intrsct->IPoint.z=z0+dz*t;
			Intrsct->Normal.z=(r0-dr*sin(Intrsct->IPoint.z))*dr*cos(Intrsct->IPoint.z);
		}
		else t=0.0;
	    }
	    else t=0.0;
	}
	else t=0.0;

	return(t);
}


double Intersect_CSG(INTERSECTION *Intrsct, LINE *Line, OBJECT *Obj)
{
	register double	t, t1, t2, ttmp;
	short	operator, inside1, inside2, leftobj;
	INTERSECTION Intersct1, Intersct2;
	LINE	NewLine;
	CSG	*CsgObject;

	CsgIntersectionAttempts++;
	CsgObject=(CSG *)Obj->Shape;

	t=0.0; operator=CsgObject->Operator;

	switch(operator) {


	/* --- Logical AND (also used for MINUS) --- */

	    case CSG_OP_AND:
		t=t1=t2=0.0;
		inside1=inside2=CSG_OUTSIDE;

		leftobj=FALSE;
		NewLine.Direction=Line->Direction;
		NewLine.Origin=Line->Origin;
		while((leftobj==FALSE)&&(inside1==CSG_OUTSIDE)) {
		    ttmp=Intersect_Object(&Intersct1,&NewLine,CsgObject->Object1);
		    if(ttmp>EPSILON) {
			t1+=ttmp;
			inside1=InsideObject(CsgObject->Object2,&Intersct1.IPoint);
			NewLine.Origin=Intersct1.IPoint;
		    }
		    else leftobj=TRUE;
		}

		leftobj=FALSE;
		NewLine.Origin=Line->Origin;
		while((leftobj==FALSE)&&(inside2==CSG_OUTSIDE)) {
		    ttmp=Intersect_Object(&Intersct2,&NewLine,CsgObject->Object2);
		    if(ttmp>EPSILON) {
			t2+=ttmp;
			inside2=InsideObject(CsgObject->Object1,&Intersct2.IPoint);
			NewLine.Origin=Intersct2.IPoint;
		    }
		    else leftobj=TRUE;
		}

		if((inside1==CSG_INSIDE)&&(inside2==CSG_OUTSIDE)) {
		    t=t1;
		    CopyIntersection(Intrsct,&Intersct1);
		}
		else if((inside1==CSG_OUTSIDE)&&(inside2==CSG_INSIDE)) {
		    t=t2;
		    CopyIntersection(Intrsct,&Intersct2);
		}
		else if((inside1==CSG_INSIDE)&&(inside2==CSG_INSIDE)) {
		    if(t1<t2) {
			t=t1;
			CopyIntersection(Intrsct,&Intersct1);
		    } else {
			t=t2;
			CopyIntersection(Intrsct,&Intersct2);
		    }
		}

		break;


	/* --- Logical OR (single surface) --- */

	    case CSG_OP_OR:
		t=t1=t2=0.0;
		inside1=inside2=CSG_INSIDE;

		leftobj=FALSE;
		NewLine.Direction=Line->Direction;
		NewLine.Origin=Line->Origin;
		while((leftobj==FALSE)&&(inside1==CSG_INSIDE)) {
		    ttmp=Intersect_Object(&Intersct1,&NewLine,CsgObject->Object1);
		    if(ttmp>EPSILON) {
			t1+=ttmp;
			inside1=InsideObject(CsgObject->Object2,&Intersct1.IPoint);
			NewLine.Origin=Intersct1.IPoint;
		    }
		    else leftobj=TRUE;
		}

		leftobj=FALSE;
		NewLine.Origin=Line->Origin;
		while((leftobj==FALSE)&&(inside2==CSG_INSIDE)) {
		    ttmp=Intersect_Object(&Intersct2,&NewLine,CsgObject->Object2);
		    if(ttmp>EPSILON) {
			t2+=ttmp;
			inside2=InsideObject(CsgObject->Object1,&Intersct2.IPoint);
			NewLine.Origin=Intersct2.IPoint;
		    }
		    else leftobj=TRUE;
		}

		if((inside1==CSG_OUTSIDE)&&(inside2==CSG_INSIDE)) {
		    t=t1;
		    CopyIntersection(Intrsct,&Intersct1);
		}
		else if((inside1==CSG_INSIDE)&&(inside2==CSG_OUTSIDE)) {
		    t=t2;
		    CopyIntersection(Intrsct,&Intersct2);
		}
		else if((inside1==CSG_OUTSIDE)&&(inside2==CSG_OUTSIDE)) {
		    if(t1<t2) {
			t=t1;
			CopyIntersection(Intrsct,&Intersct1);
		    } else {
			t=t2;
			CopyIntersection(Intrsct,&Intersct2);
		    }
		}

		break;


	/* --- PLUS (closest hit, simple OR) --- */

	    case CSG_OP_PLUS:
		t1=Intersect_Object(&Intersct1,Line,CsgObject->Object1);
		t2=Intersect_Object(&Intersct2,Line,CsgObject->Object2);
		if(t1>EPSILON) {
			if(t1<t2) t=t1;
			else if(t2>EPSILON) t=t2;
			else t=t1;
		}
		else if(t2>EPSILON) t=t2;

		if(t>EPSILON) {
		    if(fabs(t-t1)<EPSILON) CopyIntersection(Intrsct,&Intersct1);
		    else      CopyIntersection(Intrsct,&Intersct2);
		}

		break;


	    default:
		break;

	}

	if(t>0) CsgIntersections++;

	return(t);
}



/**********************************************************************
 *
 *	Test for bounding intersections
 *
 **********************************************************************/

int CheckForBounding(LINE *Line, OBJECT *Obj)
{
	switch(Obj->BoundType) {
	   case BOUND_BOX:
		return(CheckForBoundingBox(Line,(BOX *)Obj->Bound));
		break;
	   case BOUND_SPHERE:
		return(CheckForBoundingSphere(Line,(SPHERE *)Obj->Bound));
		break;
	   case BOUND_NONE:
	   default:
		return(TRUE);
		break;
	}
}

int CheckForBoundingBox(LINE *Line, BOX *Box)
{
	register double	ttmp;
	double	lvx,lvy,lvz,lx0,ly0,lz0;
	register double	ipx,ipy,ipz;

	lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
	lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;

	if(lvz!=0.0) {
	    ttmp=(Box->Corners[0].z-lz0)/lvz;	/* xy-plane, z-min */
	    if(ttmp>EPSILON) {
		ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y))
			return(TRUE);
	    }
	    ttmp=(Box->Corners[1].z-lz0)/lvz;	/* xy-plane, z-max */
	    if(ttmp>EPSILON) {
		ipx=lx0+lvx*ttmp; ipy=ly0+lvy*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y))
			return(TRUE);
	    }
	}
	if(lvy!=0.0) {
	    ttmp=(Box->Corners[0].y-ly0)/lvy;	/* xz-plane, y-min */
	    if(ttmp>EPSILON) {
		ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
			return(TRUE);
	    }
	    ttmp=(Box->Corners[1].y-ly0)/lvy;	/* xz-plane, y-max */
	    if(ttmp>EPSILON) {
		ipx=lx0+lvx*ttmp; ipz=lz0+lvz*ttmp;
		if((ipx>=Box->Corners[0].x)&&(ipx<=Box->Corners[1].x)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
			return(TRUE);
	    }
	}
	if(lvx!=0.0) {
	    ttmp=(Box->Corners[0].x-lx0)/lvx;	/* yz-plane, x-min */
	    if(ttmp>EPSILON) {
		ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
		if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
			return(TRUE);
	    }
	    ttmp=(Box->Corners[1].x-lx0)/lvx;	/* yz-plane, x-max */
	    if(ttmp>EPSILON) {
		ipy=ly0+lvy*ttmp; ipz=lz0+lvz*ttmp;
		if((ipy>=Box->Corners[0].y)&&(ipy<=Box->Corners[1].y)&&(ipz>=Box->Corners[0].z)&&(ipz<=Box->Corners[1].z))
			return(TRUE);
	    }
	}

	return(FALSE);
}



int CheckForBoundingSphere(LINE *Line, SPHERE *Sphere)
{
	double	vx,vy,vz,r;
	register double	sroot,vydy,vzdz,vxsqr,vysqr,vzsqr;
	double	dx,dy,dz;
	double	dxsqr,dysqr,dzsqr;

	vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
	r=Sphere->r;

	dx=Line->Origin.x-Sphere->Centre.x;
	dy=Line->Origin.y-Sphere->Centre.y;
	dz=Line->Origin.z-Sphere->Centre.z;
	vydy=vy*dy; vzdz=vz*dz;
	dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz;
	vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;

	sroot=2*(vx*dx*(vydy+vzdz)+vydy*vzdz)+r*r*(vxsqr+vysqr+vzsqr);
	sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);

	if(sroot>=0.0) return(TRUE);
	else           return(FALSE);
}




/*  Misc. */

void CopyIntersection(INTERSECTION *i1, INTERSECTION *i2)
{
	i1->IPoint=i2->IPoint;
	i1->Normal=i2->Normal;
	i1->Texture=i2->Texture;
	i1->Shadow=i2->Shadow;
}