unit quadrics;

{$I init.inc}

{*Unit Quadrics**********************************************************

 Quadric rountines!

 Created            : 24/05/93
 Last Change        : 13/06/93
 Revisions          : none

 ************************************************************************}

interface

uses globals;

type
  PQuadric = ^TQuadric;
  TQuadric = Object(TPrim)
    QA,QB,QC,
    QD,QE,QF,
    QG,QH,QI,QJ : Flt;
    NB, NC, ND, NF, NG, NI  : Flt;
    constructor Init(A,B,C,D,E,F,G,H,I,J : Flt);
    function Intersect(Ray : TRay; var Hitlist : PHitlist) : byte; virtual;
    procedure Normal(P : Point; var N : Vec); virtual;
    destructor Done; virtual;
  end;

implementation

{$I vectors.inc }

{*Quadric Routines**********************************************************}

constructor TQuadric.Init;

begin
  TPrim.Init;

 QA := A; QB := B; QC := C;
 QD := D; QE := E; QF := F;
 QG := G; QH := H; QI := I; QJ := J;

  { Pre-computed constants **** optimizations! **** }

  NB := 2.0 * QB;
  NC := 2.0 * QC;
  ND := 2.0 * QD;
  NF := 2.0 * QF;
  NG := 2.0 * QG;
  NI := 2.0 * QI;

  {$IFDEF Debug}
    Writeln(con,'Initialised Quadric ... ');
    Writeln(con);
  {$ENDIF}

end;

function TQuadric.Intersect;

{************************************************************************

 Source: Essential Ray Tracing Algorithms, Eric Haines.
     in: An Introduction To Raytracing, p 72.

 ************************************************************************}

var
 Aq,
 NBq,
 Cq,
 Ka,
 Kb,
 KA_KB,
 SQRT_KA_KB,
 t0,
 t1  : Flt;
 roots : Byte;


begin

  roots := 0;

  {$IFDEF Debug}
     Writeln(con,'Intersecting Quadric...');
     writeln(con);
  {$ENDIF}

  stats.qdrintersects := stats.qdrintersects + 1;

  with Ray do
  begin

    { P[] = Ray.P[] (X0)    D[] = Ray.D[] (Xd) }

    Aq := D[0] * ( QA * D[0] + NB * D[1] + NC * D[2]) +
                      D[1] * ( QE * D[1] + NF * D[2]) +
                                       QH * Sqr(D[2]);

    NBq := D[0] * ( QA * P[0] + QB * P[1] + QC * P[2] + QD) +
           D[1] * ( QB * P[0] + QE * P[1] + QF * P[2] + QG) +
           D[2] * ( QC * P[0] + QF * P[1] + QH * P[2] + QI);

    Cq := (P[0] * ( (QA * P[0]) + (NB * P[1]) + (NC * P[2]) + ND)) +
                        (P[1] * ( (QE * P[1]) + (NF * P[2]) + NG)) +
                                      (P[2] * ( (QH * P[2]) + NI)) + QJ;

  end;

  if Aq = 0 then
  begin
    t0 := -(Cq/(NBq*2));
    Hitlist^.Add(New(PHitItem,Init(t0,0,@self)));
    {$IFDEF Debug}
        writeln(con, 'Found intersection t1 ',t1:7:3); writeln(con);
    {$ENDIF Debug}
    stats.qdrintersectssucceeded := stats.qdrintersectssucceeded + 1;
    Intersect := 1;
    exit;
  end;

  Ka := -(NBq/Aq);
  Kb := Cq/Aq;

  KA_KB := sqr(Ka) - Kb;

  if KA_KB < 0 then begin intersect := 0; exit; end;

  SQRT_KA_KB := sqrt(KA_KB);

  t1 := Ka + SQRT_KA_KB;
  if t1>rayeps then
  begin
    Hitlist^.Add(New(PHititem,Init(t1,0,@self)));
    {$IFDEF Debug}
      writeln(con, 'Found intersection t1 ',t1:7:3); writeln(con);
    {$ENDIF Debug}
    stats.qdrintersectssucceeded := stats.qdrintersectssucceeded + 1;
    inc(roots);
  end;

  t0 := Ka - SQRT_KA_KB;
  if t0>rayeps then
  begin
    Hitlist^.Add(New(PHititem,Init(t0,0,@self)));
    {$IFDEF Debug}
      writeln(con, 'Found intersection t0 ',t0:7:3); writeln(con);
    {$ENDIF Debug}
    inc(roots);
  end;

  intersect:=roots;

end;

procedure TQuadric.Normal;

begin
  N[0] := QA * P[0] + QB * P[1] + QC * P[2] + QD;
  N[1] := QB * P[0] + QE * P[1] + QF * P[2] + QG;
  N[2] := QC * P[0] + QF * P[1] + QH * P[2] + QI;
  VecUnit(N);
end;

destructor TQuadric.Done;

begin
end;

end.