unit RayTr;
{ Everything-independent super ultra ray-tracing
  by Brion Vibber, 12-13-92, based on lots of stuff
  last revised 12-25-92, unless I forget to change this line.
  12-25-92: set for Floats
}
interface

{$N+}

uses Objects; { For TCollection, a super-useful thingy }

type
  Float = Double;
const
  MinReal = 5.0e-324;
  MaxReal = 1.7e308;

type
  PUniverse = ^TUniverse;
  PObject3d = ^TObject3d;
  PLightSource = ^TLightSource;
  PPoint3d = ^TPoint3d;
  PRayInfo = ^TRayInfo;
  {PPolygon3d = ^TPolygon3d;}
  {PPolygonList = ^TPolygonList;}

  TPoint3d = object
    x,y,z: Float;
    procedure SetTo(ax,ay,az: Float);
    procedure MakeUnit;
    procedure Cross(v: TPoint3d);
    function Dot(v: TPoint3d): Float;
    procedure Add(v: TPoint3d);
    procedure Sub(v: TPoint3d);
    procedure Scale(s: Float);
    function Length: Float;
    procedure RotateFrom(v: TPoint3d);
    procedure RotateTo(v: TPoint3d);
  end;
  TRay = record
    case Integer of
      0: (Origin, Direction: TPoint3d);
      1: (o,d: TPoint3d); { Shorthand versions }
      2: (a,b: TPoint3d); { "       " "      " }
  end;
  RgbTriple = record
    r,g,b: Float; { 0 - 1 }
  end;
  TRayInfo = record
    Hit: Boolean;
    Time: Float;
    GoingIn: Boolean;
    Pos,Normal: TPoint3d;
    Reflectivity,Opacity,Color: RgbTriple;
    Brightness: Float; { 0 - 1 }
    SpecularN,IndexRefraction: Float;
  end;

  TObject3d = object(TObject)
    Owner: PUniverse;
    procedure GetInfo(Ray: TRay; var Inf: TRayInfo); virtual;
    function SpecularW(Angle: Float): Float; virtual;
  end;
  TLightSource = object(TObject)
    Pos: TPoint3d;
    Color: RgbTriple;
    constructor Init(x,y,z,r,g,b: Float);
  end;

  TUniverse = object(TObject)
    Shading,Shadows,Specular,Transparency,Reflection,ScaleBase: Boolean;
    ScrWidth,ScrHeight: Integer;
    Eye,Gaze,Up,vX,vY,vU,vV,vH,vM: TPoint3d;
    DistScreen,AngleHoriz,AngleVert: Float; { Use degrees }
    BackColor,Ambient: RgbTriple;
    ReflectCount,MaxReflect: Integer;
    Items,Lights: TCollection;
    constructor Init(aScrWidth,aScrHeight,aMaxReflect: Integer;
      anEyeX,anEyeY,anEyeZ,
      aGazeX,aGazeY,aGazeZ,anUpX,anUpY,anUpZ,aDistScreen,
      anAngleHoriz,anAngleVert,anAmbientR,anAmbientG,anAmbientB,
      aBackR,aBackG,aBackB: Float);
    destructor Done; virtual;
    procedure TraceRay(Ray: TRay; var Col: RgbTriple);
    procedure TracePoint(x,y: Integer; var Col: RgbTriple);
    procedure Insert(Item: PObject3d);
    procedure Delete(Item: PObject3d);
    procedure InsertLight(Item: PLightSource);
    procedure DeleteLight(Item: PLightSource);
    procedure DeleteAllObjects;
    procedure DeleteAllLights;
    procedure DeleteAll;
    procedure FreeAllObjects;
    procedure FreeAllLights;
    procedure FreeAll;
  end;

  PSphere = ^TSphere;
  PPlane = ^TPlane;
  PPlanePts = ^TPlanePts;
  PCylinder = ^TCylinder;
  PCone = ^TCone;

  TSphere = object(TObject3d)
    Pos: TPoint3d;
    r,r2: Float;
    constructor Init(xa,ya,za,ra: Float);
    procedure GetInfo(Ray: TRay; var Inf: TRayInfo); virtual;
  end;
  TPlane = object(TObject3d)
    Pos,Normal: TPoint3d;
    constructor Init(xa,ya,za,vx,vy,vz: Float);
    procedure GetInfo(Ray: TRay; var Inf: TRayInfo); virtual;
  end;
  TPlanePts = object(TPlane)
    constructor Init(x1,y1,z1,x2,y2,z2,x3,y3,z3: Float);
  end;
  TCylinder = object(TObject3d)
    Pos,Normal: TPoint3d;
    r,r2,h: Float;
    Ends: ShortInt; { 0 = no, 1 = bottom, 2 = top, 3 = both }
    constructor Init(xa,ya,za,vx,vy,vz,ra,ha: Float; aends: ShortInt);
    procedure GetInfo(Ray: TRay; var Inf: TRayInfo); virtual;
  end;
  TCone = object(TObject3d)
    Pos,Normal: TPoint3d;
    ra,rb,ra2,rb2,h,th: Float;
    Ends: ShortInt; { 0 = no, 1 = bottom, 2 = top, 3 = both }
    constructor Init(xa,ya,za,vx,vy,vz,ra1,rb1,ha: Float; aends: ShortInt);
    procedure GetInfo(Ray: TRay; var Inf: TRayInfo); virtual;
  end;

procedure GetSphereInfo(Ray: TRay; Pos: TPoint3d; r: Float; var close,farther: Float);
procedure GetCylinderInfo(Ray: TRay; Pos,Normal: TPoint3d; r,h: Float; ends: ShortInt;
  var close,farther,end1,end2: Float);
procedure GetConeInfo(Ray: TRay; Pos,Normal: TPoint3d; ra,rb,h: Float; ends: ShortInt;
  var close,farther,end1,end2: Float);
procedure GetPlaneInfo(Ray: TRay; Pos,Normal: TPoint3d; var time: Float);

procedure MakeTrip(r,g,b: Float; var c: RgbTriple);
function Power(x,y: Float): Float;
function Log(x,y: Float): Float;
function ArcSin(x: Float): Float;
function ArcCos(x: Float): Float;
function Tan(x: Float): Float;

implementation

procedure MakeTrip(r,g,b: Float; var c: RgbTriple);
begin
  c.r := r; c.g := g; c.b := b;
end;

function Power(x,y: Float): Float;
var
  i: Integer;
  x1: Float;
begin
  {Power := Exp(y * Ln(x));}
  x1 := x;
  for i := 1 to Trunc(y)-1 do
    x1 := x1 * x;
  Power := x;
end;

function Log(x,y: Float): Float;
begin
  Log := Ln(x) / Ln(y);
end;

function ArcSin(x: Float): Float;
begin
  ArcSin := ArcTan(x / Sqrt(1 - (x*x)));
end;

function ArcCos(x: Float): Float;
begin
  ArcCos := ArcTan(Sqrt(1 - Sqr(x)) / x);
end;

function Tan(x: Float): Float;
begin
  Tan := Sin(x) / Cos(x);
end;

procedure TPoint3d.SetTo(ax,ay,az: Float);
begin
  x := ax; y := ay; z := az;
end;

procedure TPoint3d.MakeUnit;
var
  d: Float;
begin
  d := Sqrt(Sqr(x)+Sqr(y)+Sqr(z));
  if d <> 0 then begin
    x := x / d;
    y := y / d;
    z := z / d;
  end;
end;

procedure TPoint3d.Cross(v: TPoint3d);
var
  xa,ya: Float;
begin
  xa := y*v.z - z*v.y;
  ya := z*v.x - x*v.z;
  z := x*v.y - y*v.x;
  x := xa; y := ya;
end;

function TPoint3d.Dot(v: TPoint3d): Float;
begin
  Dot := x*v.x + y*v.y + z*v.z;
end;

procedure TPoint3d.Add(v: TPoint3d);
begin
  x := x + v.x;
  y := y + v.y;
  z := z + v.z;
end;

procedure TPoint3d.Sub(v: TPoint3d);
begin
  x := x - v.x;
  y := y - v.y;
  z := z - v.z;
end;

procedure TPoint3d.Scale(s: Float);
begin
  x := x * s;
  y := y * s;
  z := z * s;
end;

function TPoint3d.Length: Float;
begin
  Length := Sqrt(Sqr(x)+Sqr(y)+Sqr(z));
end;

procedure TPoint3d.RotateFrom(v: TPoint3d);
var
  d,xa,ya,za,cos1,sin1: Float;
begin
  v.MakeUnit;
  d := Sqrt(Sqr(v.y)+Sqr(v.z));
  cos1 := v.z/d; sin1 := v.y/d;
  ya := y*cos1-z*sin1; z := y*sin1+z*cos1; y := ya;
  cos1 := d; sin1 := -v.x;
  za := z*cos1-x*sin1; x := z*sin1+x*cos1; z := za;
end;

procedure TPoint3d.RotateTo(v: TPoint3d);
var
  d,xa,ya,za,cos1,sin1: Float;
begin
  v.MakeUnit;
  d := Sqrt(Sqr(v.y)+Sqr(v.z));
  cos1 := -v.z/d; sin1 := -v.y/d;
  ya := y*cos1-z*sin1; z := y*sin1+z*cos1; y := ya;
  cos1 := -d; sin1 := v.x;
  za := z*cos1-x*sin1; x := z*sin1+x*cos1; z := za;
end;


procedure TObject3d.GetInfo(Ray: TRay; var Inf: TRayInfo);
begin
  Inf.Hit := False;
end;

function TObject3d.SpecularW(Angle: Float): Float;
begin
  SpecularW := Angle / 90;
end;


constructor TLightSource.Init(x,y,z,r,g,b: Float);
begin
  TObject.Init;
  Pos.x := x; Pos.y := y; Pos.z := z;
  Color.r := r; Color.g := g; Color.b := b;
end;

constructor TUniverse.Init(aScrWidth,aScrHeight,aMaxReflect: Integer;
  anEyeX,anEyeY,anEyeZ,
  aGazeX,aGazeY,aGazeZ,anUpX,anUpY,anUpZ,aDistScreen,
  anAngleHoriz,anAngleVert,anAmbientR,anAmbientG,anAmbientB,
  aBackR,aBackG,aBackB: Float);
begin
  TObject.Init;
  Items.Init(5,5);
  Lights.Init(5,5);
  Shading := True; Shadows := True; Specular := True;
  Transparency := True; Reflection := True; ScaleBase := False;
  ScrWidth := aScrWidth; ScrHeight := aScrHeight;
  MaxReflect := aMaxReflect;
  Eye.x := anEyeX; Eye.y := anEyeY; Eye.z := anEyeZ;
  Gaze.x := aGazeX; Gaze.y := aGazeY; Gaze.z := aGazeZ;
  Up.x := anUpX; Up.y := anUpY; Up.z := anUpZ;
  DistScreen := aDistScreen;
  AngleHoriz := anAngleHoriz; AngleVert := anAngleVert;
  Ambient.r := anAmbientR; Ambient.g := anAmbientG; Ambient.b := anAmbientB;
  MakeTrip(aBackR,aBackG,aBackB,BackColor);
  Gaze.MakeUnit;
  Up.MakeUnit;
  vX := Gaze; vX.Cross(Up); vX.MakeUnit;
  vY := vX; vY.Cross(Gaze); vY.MakeUnit;
  vM := Gaze; vM.Scale(DistScreen); vM.Add(Eye);
  vH := vX; vH.Scale(DistScreen*Tan(AngleHoriz*pi/180));
  vV := vY; vV.Scale(DistScreen*Tan(AngleVert*pi/180));
end;

destructor TUniverse.Done;
begin
  Items.Done;
  Lights.Done;
  TObject.Done;
end;

procedure TUniverse.TraceRay(Ray: TRay; var Col: RgbTriple);
type
  InfAr = array[0..(65520 div SizeOf(TRayInfo))-1] of TRayInfo;
var
  Inf: ^InfAr;
  Info: TRayInfo;
  i,j,k,l: Integer;
  c,lc: RgbTriple;
  ra,Min,cosTheta,ir,ig,ib,d: Float;
  ndotn, idotn, idotn_div_ndotn_x2: Float;
  li,v1,v2: TPoint3d;
  r: TRay;
begin
  GetMem(Inf,Items.Count * SizeOf(TRayInfo));
  Ray.d.MakeUnit;
  for i := 0 to Items.Count-1 do
    PObject3d(Items.At(i))^.GetInfo(Ray,Inf^[i]);
  Min := MaxReal; j := -1;
  for i := 0 to Items.Count-1 do
    with Inf^[i] do
      if Hit and (Time < Min) and (Time > 0.001) then begin
        Min := Time; j := i;
      end;
  if j = -1 then begin
    Col := BackColor
  end else begin
    Info := Inf^[j];
    with Info do begin
      Normal.MakeUnit;
      c := Color;
      if Shading then begin
        ir := Brightness * Ambient.r;
        ig := Brightness * Ambient.g;
        ib := Brightness * Ambient.b;
        for i := 0 to Lights.Count-1 do begin
          li := PLightSource(Lights.At(i))^.Pos;
          lc := PLightSource(Lights.At(i))^.Color;
          Li.Sub(Pos);
          with li do
            d := Sqrt(Sqr(x)+Sqr(y)+Sqr(z));
          li.MakeUnit;
          costheta := (Normal.x*Li.x + Normal.y*Li.y + Normal.z*Li.z);
          if Shadows then begin
            r.Origin := Pos;
            r.Direction := li;
            repeat
              for k := 0 to Items.Count-1 do
                PObject3d(Items.At(k))^.GetInfo(r,Inf^[k]);
              Min := MaxReal; l := -1;
              for k := 0 to Items.Count-1 do
                with Inf^[k] do
                  if Hit and (Time < Min) and (Time < d) and
                    (Time > 0.001) and ((Opacity.r+Opacity.g+Opacity.b) > 0)
                    then begin
                    Min := Time; l := k;
                  end;
              if l <> -1 then
                with Inf^[l] do begin
                  lc.r := lc.r * (1-Opacity.r);
                  lc.g := lc.g * (1-Opacity.g);
                  lc.b := lc.b * (1-Opacity.b);
                  r.o := Pos;
                end;
            until l = -1;
          end;
          if Specular and (costheta > 0) then begin
            { Don't bother if light is behind }
            v1 := Ray.d;
            ndotn := (normal.x * normal.x +
	              normal.y * normal.y +
	              normal.z * normal.z);
            idotn := (normal.x * li.x +
	              normal.y * li.y +
	              normal.z * li.z);
            idotn_div_ndotn_x2 := (2.0 * (idotn) / ndotn);
            v2.x := li.x - idotn_div_ndotn_x2 * Normal.x;
            v2.y := li.y - idotn_div_ndotn_x2 * Normal.y;
            v2.z := li.z - idotn_div_ndotn_x2 * Normal.z;
            v2.MakeUnit;
            ra := v1.Dot(v2);
            if ra > 0 then costheta := costheta + (
              PObject3d(Items.At(j))^.SpecularW(ArcCos(costheta)*180/pi) *
              Power(ra,SpecularN));
          end;
          if (costheta < 0) then begin
            costheta := -costheta;
            lc.r := lc.r * (1-Opacity.r);
            lc.g := lc.g * (1-Opacity.g);
            lc.b := lc.b * (1-Opacity.b);
          end;
          Ir := Ir + ((Brightness * lc.r) / (d + 0.001)) * costheta;
          Ig := Ig + ((Brightness * lc.g) / (d + 0.001)) * costheta;
          Ib := Ib + ((Brightness * lc.b) / (d + 0.001)) * costheta;
        end;
        c.r := c.r - 1 + Ir;
        c.g := c.g - 1 + Ig;
        c.b := c.b - 1 + Ib;
      end;
      if Reflection and (ReflectCount < MaxReflect) and
        ((Reflectivity.r+Reflectivity.g+Reflectivity.b) > 0) then begin
        ndotn := (normal.x * normal.x +
	          normal.y * normal.y +
	          normal.z * normal.z);
        idotn := (normal.x * ray.d.x +
	          normal.y * ray.d.y +
	          normal.z * ray.d.z);
        idotn_div_ndotn_x2 := (2.0 * (idotn) / ndotn);
        r.Origin := Pos;
        r.d.x := Ray.d.x - idotn_div_ndotn_x2 * Normal.x;
        r.d.y := Ray.d.y - idotn_div_ndotn_x2 * Normal.y;
        r.d.z := Ray.d.z - idotn_div_ndotn_x2 * Normal.z;
        Inc(ReflectCount);
        TraceRay(r,lc);
        Dec(ReflectCount);
        if ScaleBase then begin
          c.r := (c.r*(1-Reflectivity.r)) + (Reflectivity.r * lc.r);
          c.g := (c.g*(1-Reflectivity.g)) + (Reflectivity.g * lc.g);
          c.b := (c.b*(1-Reflectivity.b)) + (Reflectivity.b * lc.b);
        end else begin
          c.r := c.r + (Reflectivity.r * lc.r);
          c.g := c.g + (Reflectivity.g * lc.g);
          c.b := c.b + (Reflectivity.b * lc.b);
        end;
      end;
      if Transparency and (ReflectCount < MaxReflect) and
        (Opacity.r < 1) and (Opacity.g < 1) and (Opacity.b < 1) then begin
        if GoingIn then
          ra := 1 / IndexRefraction
        else
          ra := IndexRefraction;
        idotn := -(Normal.x*Ray.d.x +
                   Normal.y*Ray.d.y +
                   Normal.z*Ray.d.z);
        d := 1 + Sqr(ra)*(Sqr(idotn)-1);
        r.o := Pos;
        if d >= 0 then begin
          Min := ra*idotn - Sqrt(d);
          r.d.x := Ray.d.x*ra + Normal.x*Min;
          r.d.y := Ray.d.y*ra + Normal.y*Min;
          r.d.z := Ray.d.z*ra + Normal.z*Min;
        end else begin
          ndotn := (normal.x * normal.x +
                    normal.y * normal.y +
                    normal.z * normal.z);
          idotn := (normal.x * ray.d.x +
                    normal.y * ray.d.y +
                    normal.z * ray.d.z);
          idotn_div_ndotn_x2 := (2.0 * (idotn) / ndotn);
          r.d.x := Ray.d.x - idotn_div_ndotn_x2 * Normal.x;
          r.d.y := Ray.d.y - idotn_div_ndotn_x2 * Normal.y;
          r.d.z := Ray.d.z - idotn_div_ndotn_x2 * Normal.z;
        end;
        Inc(ReflectCount);
        TraceRay(r,lc);
        Dec(ReflectCount);
        if ScaleBase then begin
          c.r := (c.r*Opacity.r) + ((1-Opacity.r) * lc.r);
          c.g := (c.g*Opacity.g) + ((1-Opacity.g) * lc.g);
          c.b := (c.b*Opacity.b) + ((1-Opacity.b) * lc.b);
        end else begin
          c.r := c.r + ((1-Opacity.r) * lc.r);
          c.g := c.g + ((1-Opacity.g) * lc.g);
          c.b := c.b + ((1-Opacity.b) * lc.b);
        end;
      end;
      Col := c;
    end;
  end;
  FreeMem(Inf,Items.Count * SizeOf(TRayInfo));
end;

procedure TUniverse.TracePoint(x,y: Integer; var Col: RgbTriple);
var
  Ray: TRay;
  v,h: TPoint3d;
begin
  Ray.Origin := Eye;
  Ray.d := vM; Ray.d.Sub(Eye);
  h := vH; v := vV;
  h.Scale(-(2*(x/ScrWidth)-1));
  v.Scale(2*(y/ScrHeight)-1);
  Ray.d.Add(h);
  Ray.d.Add(v);
  TraceRay(Ray,Col);
end;

procedure TUniverse.Insert(Item: PObject3d);
begin
  Items.Insert(Item);
  Item^.Owner := @Self;
end;

procedure TUniverse.Delete(Item: PObject3d);
begin
  Items.Delete(Item);
  Item^.Owner := nil;
end;

procedure TUniverse.InsertLight(Item: PLightSource);
begin
  Lights.Insert(Item);
end;

procedure TUniverse.DeleteLight(Item: PLightSource);
begin
  Lights.Delete(Item);
end;

procedure TUniverse.DeleteAllObjects;
var i: Integer;
begin
  Items.DeleteAll;
end;

procedure TUniverse.DeleteAllLights;
begin
  Lights.DeleteAll;
end;

procedure TUniverse.DeleteAll;
begin
  Items.DeleteAll;
  Lights.DeleteAll;
end;

procedure TUniverse.FreeAllObjects;
begin
  Items.FreeAll;
end;

procedure TUniverse.FreeAllLights;
begin
  Lights.FreeAll;
end;

procedure TUniverse.FreeAll;
begin
  Items.FreeAll;
  Lights.FreeAll;
end;

constructor TSphere.Init(xa,ya,za,ra: Float);
begin
  TObject3d.Init;
  Pos.SetTo(xa,ya,za);
  r := ra; r2 := Sqr(ra);
end;

procedure TSphere.GetInfo(Ray: TRay; var Inf: TRayInfo);
var
  a, b, c, t1, t2, t3, close, farther: Float;
begin
  a := ray.d.x * ray.d.x + ray.d.y * ray.d.y + ray.d.z * ray.d.z;
  close := -1.0;
  farther := -1.0;
  if a <> 0 then begin
    b := 2.0 * ((ray.o.x - Pos.x) * ray.d.x
        + (ray.o.y - Pos.y) * ray.d.y
        + (ray.o.z - Pos.z) * ray.d.z);
    c := (ray.o.x - Pos.x) * (ray.o.x - Pos.x)
        + (ray.o.y - Pos.y) * (ray.o.y - Pos.y)
        + (ray.o.z - Pos.z) * (ray.o.z - Pos.z) - r2;
    t1 := b * b - 4.0 * a * c;
    if t1 > 0 then begin
      t2 := sqrt(t1);
      t3 := 2.0 * a;
      close := -(b + t2) / t3;
      farther := -(b - t2) / t3;
    end;
  end;
  if (close <= 0.001) and (farther > 0.001) then begin
    Inf.Time := farther; Inf.GoingIn := False;
  end else
    if (close > 0.001) and (farther <= 0.001) then begin
      Inf.Time := close; Inf.GoingIn := False;
    end else begin
      Inf.GoingIn := True;
      if close < farther then
        Inf.Time := close
      else
        Inf.Time := farther;
    end;
  Inf.Hit := (Inf.Time > 0.001);
  if not Inf.Hit then Exit;
  Inf.Pos.SetTo(
    ray.o.x+ray.d.x*Inf.Time,
    ray.o.y+ray.d.y*Inf.Time,
    ray.o.z+ray.d.z*Inf.Time);
  Inf.Normal := Inf.Pos;
  Inf.Normal.Sub(Pos);
  if not Inf.GoingIn then Inf.Normal.Scale(-1);
  MakeTrip(0.9,0.8,0.7,Inf.Color);
  Inf.Brightness := 1;
  MakeTrip(0.3,0.3,0.3,Inf.Reflectivity);
  MakeTrip(1,1,1,Inf.Opacity);
  Inf.SpecularN := 100;
  Inf.IndexRefraction := 0.95;
end;

constructor TPlane.Init(xa,ya,za,vx,vy,vz: Float);
begin
  TObject3d.Init;
  Pos.SetTo(xa,ya,za);
  Normal.SetTo(vx,vy,vz);
  Normal.MakeUnit;
end;

procedure TPlane.GetInfo(Ray: TRay; var Inf: TRayInfo);
var
  p1, p2, p3: Float;
begin
  p1 := Pos.x * Normal.x + Pos.y * Normal.y + Pos.z * Normal.z;
  p2 := ray.o.x * Normal.x + ray.o.y * Normal.y + ray.o.z * Normal.z;
  p3 := ray.d.x * Normal.x + ray.d.y * Normal.y + ray.d.z * Normal.z;
  if p3 = 0 then p3 := 0.001;
  Inf.Time := (p1-p2)/p3;
  Inf.Hit := (Inf.Time > 0.001);
  if not Inf.Hit then Exit;
  Inf.Normal := Normal;
  if Ray.d.Dot(Inf.Normal) > 0 then Inf.Normal.Scale(-1);
  Inf.Pos.SetTo(
    ray.o.x+ray.d.x*Inf.Time,
    ray.o.y+ray.d.y*Inf.Time,
    ray.o.z+ray.d.z*Inf.Time);
  MakeTrip(0.75,0.75,0.75,Inf.Color);
  Inf.GoingIn := True;
  Inf.Brightness := 1;
  MakeTrip(0.1,0.1,0.1,Inf.Reflectivity);
  MakeTrip(1,1,1,Inf.Opacity);
  Inf.SpecularN := 50;
  Inf.IndexRefraction := 0.99;
end;

constructor TPlanePts.Init(x1,y1,z1,x2,y2,z2,x3,y3,z3: Float);
begin
  TPlane.Init(x1,y1,z1,
    y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2),
    z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2),
    x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
end;

constructor TCylinder.Init(xa,ya,za,vx,vy,vz,ra,ha: Float; aends: ShortInt);
begin
  TObject3d.Init;
  Pos.SetTo(xa,ya,za);
  Normal.SetTo(vx,vy,vz);
  Normal.MakeUnit;
  r := ra; r2 := Sqr(ra);
  h := ha;
  ends := aends;
end;

procedure TCylinder.GetInfo(Ray: TRay; var Inf: TRayInfo);
var
  r1: TRay;
  a,b,c,t1,t2,t3,close1,farther1,close,farther,end1,end2: Float;
  p1,p2: TPoint3d;
begin
  Inf.Hit := False;
  r1 := Ray;
  r1.o.Sub(Pos);
  r1.o.RotateFrom(Normal);
  r1.d.RotateFrom(Normal);
  if (ends and 1) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end1 := (-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end1)+
      Sqr(r1.o.y+r1.d.y*end1)+
      Sqr(r1.o.z+r1.d.z*end1)) > r2 then end1 := -1;
  end else
    end1 := -1;
  if (ends and 2) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end2 := (h-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end2)+
      Sqr(r1.o.y+r1.d.y*end2)+
      Sqr(r1.o.z+r1.d.z*end2-h)) > r2 then end2 := -1;
  end else
    end2 := -1;
  a := r1.d.x * r1.d.x + r1.d.y * r1.d.y;
  close := -1.0;
  farther := -1.0;
  if a <> 0 then begin
    b := 2.0 * (r1.o.x * r1.d.x
        + r1.o.y * r1.d.y);
    c := r1.o.x * r1.o.x
        + r1.o.y * r1.o.y - r2;
    t1 := b * b - 4.0 * a * c;
    if t1 > 0 then begin
      t2 := sqrt(t1);
      t3 := 2.0 * a;
      close1 := -(b + t2) / t3;
      farther1 := -(b - t2) / t3;
      if close1 > 0 then begin
        p1 := r1.d;
        p1.Scale(close1);
        close := p1.Length;
        p1.Add(r1.o);
        if (p1.z < 0) or (p1.z > h) then close := -1;
      end;
      if farther1 > 0 then begin
        p2 := r1.d;
        p2.Scale(farther1);
        farther := p2.Length;
        p2.Add(r1.o);
        if (p2.z < 0) or (p2.z > h) then farther := -1;
      end;
    end;
  end;
  Inf.Time := MaxReal;
  if (close < Inf.Time) and (close > 0.001) then Inf.Time := close;
  if (farther < Inf.Time) and (farther > 0.001) then Inf.Time := farther;
  if (end1 < Inf.Time) and (end1 > 0.001) then Inf.Time := end1;
  if (end2 < Inf.Time) and (end2 > 0.001) then Inf.Time := end2;
  Inf.Hit := Inf.Time < MaxReal;
  if not Inf.Hit then Exit;
  Inf.Pos := Ray.d;
  Inf.Pos.Scale(Inf.Time);
  Inf.Pos.Add(Ray.o);
  Inf.Normal := Normal;
  if (Inf.Time = close) or (Inf.Time = farther) then begin
    if Inf.Time = close then
      Inf.Normal.Scale(p1.z)
    else
      Inf.Normal.Scale(p2.z);
    Inf.Normal.Add(Pos);
    Inf.Normal.Sub(Inf.Pos);
  end else begin
    if Inf.Time = end2 then Inf.Normal.Scale(-1);
  end;
  Inf.GoingIn := Ray.d.Dot(Inf.Normal) > 0;
  if Inf.GoingIn then Inf.Normal.Scale(-1);
  MakeTrip(0.7,0.8,0.9,Inf.Color);
  Inf.Brightness := 1;
  MakeTrip(0,0,0,Inf.Reflectivity);
  MakeTrip(1,1,1,Inf.Opacity);
  Inf.SpecularN := 100;
  Inf.IndexRefraction := 0.95;
end;

constructor TCone.Init(xa,ya,za,vx,vy,vz,ra1,rb1,ha: Float; aends: ShortInt);
begin
  TObject3d.Init;
  Pos.SetTo(xa,ya,za);
  Normal.SetTo(vx,vy,vz);
  Normal.MakeUnit;
  ra := ra1; ra2 := Sqr(ra);
  rb := rb1; rb2 := Sqr(rb);
  h := ha;
  th := h+((rb*h)/(ra-rb));
  ends := aends;
end;

procedure TCone.GetInfo(Ray: TRay; var Inf: TRayInfo);
var
  r1: TRay;
  v1,v2,a,b,c,t1,t2,t3,close1,farther1,close,farther,end1,end2: Float;
  p1,p2: TPoint3d;
begin
  Inf.Hit := False;
  r1 := Ray;
  r1.o.Sub(Pos);
  r1.o.RotateFrom(Normal);
  r1.d.RotateFrom(Normal);
  if (ends and 1) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end1 := (-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end1)+
      Sqr(r1.o.y+r1.d.y*end1)+
      Sqr(r1.o.z+r1.d.z*end1)) > ra2 then end1 := -1;
  end else
    end1 := -1;
  if (ends and 2) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end2 := (h-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end2)+
      Sqr(r1.o.y+r1.d.y*end2)+
      Sqr(r1.o.z+r1.d.z*end2-h)) > rb2 then end2 := -1;
  end else
    end2 := -1;
  a := r1.d.x * r1.d.x + r1.d.y * r1.d.y;
  (*
    xý+yý=(1-(z/th))rý
    take (1-(z/th))*ra2 and use as r2 in cyl stuff
  *)
  { get p1.z, p2.z }
  v1 := (1-(p1.z/th))*ra2;
  v2 := (1-(p2.z/th))*ra2;
  a := r1.d.x * r1.d.x + r1.d.y * r1.d.y;
  close := -1.0;
  farther := -1.0;
  if a <> 0 then begin
    b := 2.0 * (r1.o.x * r1.d.x
        + r1.o.y * r1.d.y);
    c := r1.o.x * r1.o.x
        + r1.o.y * r1.o.y - ra2;
    t1 := b * b - 4.0 * a * c;
    if t1 > 0 then begin
      t2 := sqrt(t1);
      t3 := 2.0 * a;
      close1 := -(b + t2) / t3;
      farther1 := -(b - t2) / t3;
      if close1 > 0 then begin
        p1 := r1.d;
        p1.Scale(close1);
        close := p1.Length;
        p1.Add(r1.o);
        if (p1.z < 0) or (p1.z > h) then close := -1;
      end;
      if farther1 > 0 then begin
        p2 := r1.d;
        p2.Scale(farther1);
        farther := p2.Length;
        p2.Add(r1.o);
        if (p2.z < 0) or (p2.z > h) then farther := -1;
      end;
    end;
  end;
  Inf.Time := MaxReal;
  if (close < Inf.Time) and (close > 0.001) then Inf.Time := close;
  if (farther < Inf.Time) and (farther > 0.001) then Inf.Time := farther;
  if (end1 < Inf.Time) and (end1 > 0.001) then Inf.Time := end1;
  if (end2 < Inf.Time) and (end2 > 0.001) then Inf.Time := end2;
  Inf.Hit := Inf.Time < MaxReal;
  if not Inf.Hit then Exit;
  Inf.Pos := Ray.d;
  Inf.Pos.Scale(Inf.Time);
  Inf.Pos.Add(Ray.o);
  Inf.Normal := Normal;
  if (Inf.Time = close) or (Inf.Time = farther) then begin
    if Inf.Time = close then
      Inf.Normal.Scale(p1.z)
    else
      Inf.Normal.Scale(p2.z);
    Inf.Normal.Add(Pos);
    Inf.Normal.Sub(Inf.Pos);
  end else begin
    if Inf.Time = end2 then Inf.Normal.Scale(-1);
  end;
  Inf.GoingIn := Ray.d.Dot(Inf.Normal) > 0;
  if Inf.GoingIn then Inf.Normal.Scale(-1);
  MakeTrip(0.7,0.8,0.9,Inf.Color);
  Inf.Brightness := 1;
  MakeTrip(0,0,0,Inf.Reflectivity);
  MakeTrip(1,1,1,Inf.Opacity);
  Inf.SpecularN := 100;
  Inf.IndexRefraction := 0.95;
end;

procedure GetSphereInfo(Ray: TRay; Pos: TPoint3d; r: Float; var close,farther: Float);
var
  r2, a, b, c, t1, t2, t3: Float;
begin
  r2 := Sqr(r);
  a := ray.d.x * ray.d.x + ray.d.y * ray.d.y + ray.d.z * ray.d.z;
  close := -1.0;
  farther := -1.0;
  if a <> 0 then begin
    b := 2.0 * ((ray.o.x - Pos.x) * ray.d.x
        + (ray.o.y - Pos.y) * ray.d.y
        + (ray.o.z - Pos.z) * ray.d.z);
    c := (ray.o.x - Pos.x) * (ray.o.x - Pos.x)
        + (ray.o.y - Pos.y) * (ray.o.y - Pos.y)
        + (ray.o.z - Pos.z) * (ray.o.z - Pos.z) - r2;
    t1 := b * b - 4.0 * a * c;
    if t1 > 0 then begin
      t2 := sqrt(t1);
      t3 := 2.0 * a;
      close := -(b + t2) / t3;
      farther := -(b - t2) / t3;
    end;
  end;
end;

procedure GetCylinderInfo(Ray: TRay; Pos,Normal: TPoint3d; r,h: Float; ends: ShortInt;
  var close,farther,end1,end2: Float);
var
  r1: TRay;
  r2,a,b,c,t1,t2,t3,close1,farther1: Float;
  p1,p2: TPoint3d;
begin
  r2 := Sqr(r);
  r1 := Ray;
  r1.o.Sub(Pos);
  r1.o.RotateFrom(Normal);
  r1.d.RotateFrom(Normal);
  if (ends and 1) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end1 := (-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end1)+
      Sqr(r1.o.y+r1.d.y*end1)+
      Sqr(r1.o.z+r1.d.z*end1)) > r2 then end1 := -1;
  end else
    end1 := -1;
  if (ends and 2) <> 0 then begin
    t3 := r1.d.z;
    if t3 = 0 then t3 := 0.001;
    end2 := (h-r1.o.z)/t3;
    if (Sqr(r1.o.x+r1.d.x*end2)+
      Sqr(r1.o.y+r1.d.y*end2)+
      Sqr(r1.o.z+r1.d.z*end2-h)) > r2 then end2 := -1;
  end else
    end2 := -1;
  a := r1.d.x * r1.d.x + r1.d.y * r1.d.y;
  close := -1.0;
  farther := -1.0;
  if a <> 0 then begin
    b := 2.0 * (r1.o.x * r1.d.x
        + r1.o.y * r1.d.y);
    c := r1.o.x * r1.o.x
        + r1.o.y * r1.o.y - r2;
    t1 := b * b - 4.0 * a * c;
    if t1 > 0 then begin
      t2 := sqrt(t1);
      t3 := 2.0 * a;
      close1 := -(b + t2) / t3;
      farther1 := -(b - t2) / t3;
      if close1 > 0 then begin
        p1 := r1.d;
        p1.Scale(close1);
        close := p1.Length;
        p1.Add(r1.o);
        if (p1.z < 0) or (p1.z > h) then close := -1;
      end;
      if farther1 > 0 then begin
        p2 := r1.d;
        p2.Scale(farther1);
        farther := p2.Length;
        p2.Add(r1.o);
        if (p2.z < 0) or (p2.z > h) then farther := -1;
      end;
    end;
  end;
end;

procedure GetConeInfo(Ray: TRay; Pos,Normal: TPoint3d; ra,rb,h: Float; ends: ShortInt;
  var close,farther,end1,end2: Float);
begin
end;

procedure GetPlaneInfo(Ray: TRay; Pos,Normal: TPoint3d; var time: Float);
var
  p1, p2, p3: Float;
begin
  p1 := Pos.x * Normal.x + Pos.y * Normal.y + Pos.z * Normal.z;
  p2 := ray.o.x * Normal.x + ray.o.y * Normal.y + ray.o.z * Normal.z;
  p3 := ray.d.x * Normal.x + ray.d.y * Normal.y + ray.d.z * Normal.z;
  if p3 = 0 then p3 := 0.001;
  Time := (p1-p2)/p3;
end;

end.