program raybench;
uses Windows, Math;
const
  Spheres: array [0..23] of Single =
             (0.0,-6.95149,-0.17298,6.932,
             -1.64317,-0.408679,-3.07338,1,
             1.66359,-0.256469,-5.76149,1,
             1.62878,0.840152,-5.60949,1,
             1.59027,1.74296,-5.4756,1,
             -3.78282,1.27461,-11.7239,1);
  Col: array [0..17] of Single =
             (0.2, 0.4, 0.1,
             0.5, 0.2, 0.6,
             0.7, 0.4, 0.4,
             0.2, 0.6, 0.7,
             0.1, 0.8, 0.5,
             0.5, 0.3, 0.7);
  Width = 640;
  Height = 480;
  HWidth = Width div 2;
  HHeight = Height div 2;
  SubGrid = 16;
type
  TSphere = array [0..3] of Single;
  PSphere = ^TSphere;
  RGB = record
    R, G, B: Byte;
  end;
  Vec = record
    X, Y, Z: Single;
  end;
var
  Pixels: array [0..WIDTH*HEIGHT - 1] of RGB;
function VLength(const V: Vec): Single;
begin
  Result:=Sqrt(V.X*V.X + V.Y*V.Y + V.Z*V.z);
end;
procedure Normalize(var V: Vec);
var
  VL: Single;
begin
  VL:=VLength(V);
  V.X:=V.X/VL;
  V.Y:=V.Y/VL;
  V.Z:=V.Z/VL;
end;
function RayHitSphere(const O, D: Vec; S: PSphere; var IP: Vec): Boolean;
var
  SRX, SRY, SRZ, b, C, DD, E, T0, T1, T: Single;
begin
  SRX:=O.X - S[0];
  SRY:=O.Y - S[1];
  SRZ:=O.Z - S[2];
  B:=SRX*D.X + SRY*D.Y + SRZ*D.Z;
  C:=(SRX*SRX + SRY*SRY + SRZ*SRZ) - S[3]*S[3];
  DD:=B*B - C;
  if DD > 0 then begin
    E:=Sqrt(DD);
    T0:=-B-E; T1:=-B+E;
    if T1 < T0 then T:=T1 else T:=T0;
    if T < 0.0001 then begin Result:=False; Exit; end;
    IP.X:=O.X + D.X*T;
    IP.Y:=O.Y + D.Y*T;
    IP.Z:=O.Z + D.Z*T;
    Result:=True;
  end else Result:=False;
end;
function RayHitSpheres(const O, D: Vec; var IP: Vec): Integer;
var
  Dist, CDist: Single;
  RIP: Vec;
  I: Integer;
begin
  CDist:=$FFFFFFFF;
  Result:=-1;
  I:=0;
  while I < High(Spheres) do begin
    if RayHitSphere(O, D, PSphere(@Spheres[I]), RIP) then begin
      Dist:=VLength(RIP);
      if (Result=-1) or (CDist > Dist) then begin
        CDist:=Dist;
        Result:=I div 4;
        IP:=RIP;
      end;
    end;
    Inc(I, 4);
  end;
end;
procedure CalcSubPixel(var RGB: Vec; X, Y, SX, SY: Single);
var
  DL: Single;
  SIdx: Integer;
  O, D, N, IP: Vec;
begin
  O.X:=0; O.Y:=0; O.Z:=0;
  D.X:=X/HWidth - 1 + SX;
  D.Y:=1 - Y/HHeight/(HWidth/HHeight) + SY;
  D.Z:=-1;
  Normalize(D);
  SIdx:=RayHitSpheres(O, D, IP);
  if SIdx=-1 then begin
    RGB.X:=RGB.X + 0.062;
    RGB.Y:=RGB.Y + 0.125;
    RGB.Z:=RGB.Z + 0.392;
  end else begin
    N.X:=IP.X - Spheres[SIdx*4];
    N.Y:=IP.Y - Spheres[SIdx*4 + 1];
    N.Z:=IP.Z - Spheres[SIdx*4 + 2];
    Normalize(D);
    DL:=0.2 + N.Z;
    if DL < 0.2 then DL:=0.2;
    if DL > 1 then DL:=1;
    RGB.X:=RGB.X + Col[SIdx*3]*DL;
    RGB.Y:=RGB.Y + Col[SIdx*3 + 1]*DL;
    RGB.Z:=RGB.Z + Col[SIdx*3 + 2]*DL;
  end;
end;
procedure CalcPixel(var RGB: RGB; X, Y: Integer);
var
  RGBF: Vec;
  GX, GY: Integer;
begin
  RGBF.X:=0; RGBF.Y:=0; RGBF.Z:=0;
  for GY:=0 to SubGrid - 1 do
    for GX:=0 to SubGrid - 1 do
      CalcSubPixel(RGBF, X, Y,
        GX/Width/SubGrid, GY/Width/SubGrid);
  RGB.R:=Trunc((RGBF.X/(SubGrid*SubGrid))*255.0);
  RGB.G:=Trunc((RGBF.Y/(SubGrid*SubGrid))*255.0);
  RGB.B:=Trunc((RGBF.Z/(SubGrid*SubGrid))*255.0);
end;
procedure CalcImage;
var
  X, Y: Integer;
begin
  for Y:=0 to Height - 1 do
    for X:=0 to Width - 1 do
      CalcPixel(Pixels[Y*Width + X], X, Y);
end;
procedure DumpPixels;
var
  X, Y: Integer;
begin
  Writeln('P3');
  Writeln(Width, ' ', Height);
  Writeln('255');
  for Y:=0 to Height - 1 do begin
    for X:=0 to Width - 1 do begin
      if X > 0 then Write(' ');
      Write(Pixels[Y*Width + X].R, ' ', Pixels[Y*Width + X].G, ' ', Pixels[Y*Width + X].B);
    end;
    Writeln;
  end;
end;
var
  Start: DWORD;
begin
  if ParamStr(1) <> 'dump' then begin
    Writeln('Warmup pass');
    CalcImage;
  end;
  Start:=GetTickCount;
  CalcImage;
  if ParamStr(1)='dump' then DumpPixels
  else Writeln('Elapsed time = ', ((GetTickCount - Start)/1000):0:3);
end.