program h3 (input, output);

const
   pi = 3.1415927;
   numpolys = 360;
   numpoints = 390;

type
   vector = array [1..4] of real;
   m4x4 = array [1..4] of vector;
   polytype = array [1..4] of integer;

var
   last,cl:integer;
   xrot, yrot, zrot, Tper, Tscale, Trot, mtemp: m4x4;
   xa, ya, za, xtrans, ytrans, ztrans, d, scale: real;
   p, Nav: array [1..numpoints] of vector;
   Etot: array [1..numpoints] of real;
   polys: array [1..numpolys] of polytype;
   vis: array [1..numpolys] of boolean;
   minz: array [1..numpolys] of real;
   N: array [1..numpolys] of vector;
   i, j: integer;
   c: char;
   viewer, light: vector;
   Rn, Ws, Id, Ips, k, Rp: real;
   data, parms: text;

%include '786pas.inc'

procedure add (v1, v2: vector; var vr: vector);

var 
   i: integer;

begin
   for i := 1 to 4 do
      vr [i] := v1 [i] + v2 [i]
end;


procedure sub (v1, v2: vector; var vr: vector);

var
   i: integer;

begin
   for i := 1 to 4 do
      vr [i] := v1 [i] - v2 [i]
end;


function mag (v: vector): real;

begin
   mag := sqrt (v [1] * v [1] + v [2] * v [2] + v [3] * v [3])
end;


procedure MVmult (m: m4x4; v: vector; var vr: vector);

var
   i, j: integer;

begin
   for i := 1 to 4 do
   begin
      vr [i] := 0;
      for j := 1 to 4 do
         vr [i] := vr [i] + m [i, j] * v [j]
   end
end;


procedure MMmult (m1, m2: m4x4; var mr: m4x4);

var
   i, j, k: integer;

begin
   for i := 1 to 4 do
      for j := 1 to 4 do
      begin
         mr [i, j] := 0;
         for k := 1 to 4 do
            mr [i, j] := mr [i, j] + m1 [i, k] * m2 [k, j]
      end
end;


procedure mout (m: m4x4);

var
   i, j: integer;

begin
   for i := 1 to 4 do
   begin
      for j := 1 to 4 do
         write (m [i, j]: 10: 4);
      writeln
   end
end;


procedure vout (v: vector);

var
   j: integer;

begin
   for j := 1 to 4 do
   begin
      write (v [j]: 10: 4);
      writeln
   end
end;


procedure readfile;

var
   i, j, line, vert: integer;

begin
   reset (data, 'gpract1');
   readln (data);
   readln (data);
   for i := 1 to numpoints do
   begin
      read (data, line);
      read (data, c);
      for j := 1 to 3 do
         read (data, p [i, j]);
      p [i, 4] := 1;
      readln (data)
   end;

   for i := 1 to 3 do
      readln (data);

   for i := 1 to numpolys do
   begin
      read (data, c);
      read (data, vert);
      for j := 1 to 4 do
      begin
         read (data, c);
         read (data, c);
         read (data, polys [i, j])
      end;
      readln (data)
   end;

   for i := 1 to 3 do
      readln (data);

   read (data, xtrans);
   read (data, ytrans);
   read (data, d);

   read (data, xa);
   read (data, ya);
   read (data, za)

end;


procedure setline (var m: m4x4; l: integer; a, b, c, d: real);

begin
   m [l, 1] := a; m [l, 2] := b; m [l, 3] := c; m [l, 4] := d
end;


procedure norm (var v: vector);

var
   i: integer;

begin
   for i := 1 to 3 do
      v [i] := v [i] / v [4];
   v [4] := 1
end;


function dot_prod (v1, v2: vector): real;

begin
   dot_prod := v1 [1] * v2 [1] + v1 [2] * v2 [2] + v1 [3] * v2 [3]
end;


procedure cross_prod (v1, v2: vector; var cross: vector);

begin
   cross [1] := v1 [2] * v2 [3] - v1 [3] * v2 [2];
   cross [2] := v1 [3] * v2 [1] - v1 [1] * v2 [3];
   cross [3] := v1 [1] * v2 [2] - v1 [2] * v2 [1];
   cross [4] := 1
end;


procedure draw_poly (pol: integer);

var
   i, j, pmin, x, y: integer;
   min, xs, xf, xsinc, xfinc, Es, Ef, Esinc, Efinc, Einc, E: real;
   used: array [1..4] of boolean;
   xp, yp, Ep: array [1..4] of real;


procedure plottemp (x, y: integer; E: real);

begin
{
   writeln ('### ', x:10, y:10, E:10:4);
}
   if E > 15 then E := 15 else if E < 0 then E := 0;
   cl:=round(e);
   if last<>cl then defcolors(round(e),5);
   last:=cl;
   point(1,0);
end;


procedure find_incs (start: integer);

var
   i: integer;
   xfstep, xsstep, ydiff: real;

begin
   xsinc := 1000;
   xfinc := -1000;
   for i := start to 4 do
   begin
      ydiff := yp [i] - y;
      xfstep := (xp [i] - xf) / ydiff;
      if xfstep > xfinc then
      begin
         xfinc := xfstep;
         Efinc := (Ep [i] - Ef) / ydiff
      end;
      xsstep := (xp [i] - xs) / ydiff;
      if xsstep < xsinc then
      begin
         xsinc := xsstep;
         Esinc := (Ep [i] - Es) / ydiff
      end
   end
end;


begin
{
   moveabs (round (p [polys [pol, 4], 1]), round (p [polys [pol, 4], 2]));
   for i := 1 to 4 do
      lineabs (round (p [polys [pol, i], 1]), round (p [polys [pol, i], 2]))
}


   for i := 1 to 4 do
      used [i] := false;
   for i := 1 to 4 do
   begin
      min := 1000;
      for j := 1 to 4 do
         if (p [polys [pol, j], 2] < min) and not used [j] then
         begin
            min := p [polys [pol, j], 2];
            pmin := j
         end;
      xp [i] := p [polys [pol, pmin], 1];
      yp [i] := p [polys [pol, pmin], 2];
      if i > 1 then if yp [i] = yp [i - 1] then yp [i] := yp [i] + 0.01;
{
writeln ('***xy ', xp [i], yp [i]);
}
      Ep [i] := Etot [polys [pol, pmin]];
      used [pmin] := true
   end;

   xs := xp [1];
   xf := xs;
   Es := Ep [1];
   Ef := Es;

   y := round (yp [1]);
   find_incs (2);
  opengcmbmain;
   for y := round (yp [1]) to round (yp [2]) - 1 do
   begin
      if xs <> xf then Einc := (Ef - Es) / (xf - xs) else Einc := 0;
      E := Es;
      absmov(round(xs+xtrans),y+round(ytrans));
      for x := round (xs) to round (xf) do
      begin
         plottemp (x, y, E);
         E := E + Einc
      end;
      Es := Es + Esinc;
      Ef := Ef + Efinc;
      xs := xs + xsinc;
      xf := xf + xfinc
   end;

   if round (yp [1]) = round (yp [2]) then
   begin
      if xp [1] < xp [2] then
      begin
         xf := xp [2];
         Ef := Ep [2]
      end
      else
      begin
         xs := xp [2];
         Es := Ep [2]
      end
   end
   else
   begin
      if abs (xp [2] - xs) < abs (xp [2] - xf) then
         xs := xp [2] else xf := xp [2]
   end;

   find_incs (3);
   for y := round (yp [2]) to round (yp [3]) - 1 do
   begin
      if xf <> xs then Einc := (Ef - Es) / (xf - xs) else Einc := 0;
      E := Es;
      absmov(round(xs+xtrans),y+round(ytrans));
      for x := round (xs) to round (xf) do
      begin
         plottemp (x, y, E);
         E := E + Einc
      end;
      Es := Es + Esinc;
      Ef := Ef + Efinc;
      xs := xs + xsinc;
      xf := xf + xfinc
   end;

   if round (yp [1]) = round (yp [3]) then
   begin
      if xp [1] < xp [3] then
      begin
         xf := xp [3];
         Ef := Ep [3]
      end
      else
      begin
         xs := xp [3];
         Es := Ep [3]
      end
   end
   else
   begin
      if abs (xp [3] - xs) < abs (xp [3] - xf) then
         xs := xp [3] else xf := xp [3]
   end;

   find_incs (4);
   for y := round (yp [3]) to round (yp [4]) do
   begin
      if xf <> xs then Einc := (Ef - Es) / (xf - xs) else Einc := 0;
      E := Es;
      absmov(round(xs+xtrans),y+round(ytrans));
      for x := round (xs) to round (xf) do
      begin
         plottemp (x, y, E);
         E := E + Einc
      end;
      Es := Es + Esinc;
      Ef := Ef + Efinc;
      xs := xs + xsinc;
      xf := xf + xfinc
   end;
     closegcmbmain;
end;


procedure find_surface_normals;

var
   i: integer;
   edge1, edge2: vector;

begin
   for i := 1 to numpolys do
   begin
      sub (p [polys [i, 2]], p [polys [i, 1]], edge1);
      sub (p [polys [i, 2]], p [polys [i, 3]], edge2);
      cross_prod (edge2, edge1, N [i])
   end
end;


procedure find_point_intensities;

var
   i, j: integer;
   dotLN, dotSV, Epd, Eps: real;
   L, V, S: vector;

begin
   for i := 1 to numpoints do
      for j := 1 to 4 do
         Nav [i, j] := 0;
   for i := 1 to numpolys do
      for j := 1 to 4 do
         add (Nav [polys [i, j]], N [i], Nav [polys [i, j]]);
{
writeln ('*** Nav follow ...');
for i := 1 to numpoints do vout (Nav [i]);
writeln;
}
   Epd := Id * Rp;
   for i := 1 to numpoints do
   begin
      sub (light, p [i], L);
{
writeln ('*** L ...');
vout (L);
writeln;
}
      dotLN := dot_prod (L, Nav [i]) / (mag (L) * mag (Nav [i]));
{
writeln (dot_prod (L, Nav [i]):10:4, mag (L):10:4, mag (Nav [i]):10:4);
writeln (dotLN);
}
      if dotLN <= 0.001 then Eps := 0
      else
      begin
         for j := 1 to 3 do
            S [j] := 2 * Nav [i, j] * dotLN / mag (Nav [i]) - L [j];
         sub (viewer, p [i], V);
{
writeln ('*** S & V ...');
vout (S); vout (V);
writeln;
}
         dotSV := dot_prod (S, V) / (mag (S) * mag (V));
         if dotSV <= 0 then Eps := 0
         else Eps := (Ips / (mag (V) + k)) * (Rp * dotLN +
            Ws * exp (Rn * ln (dotSV)));
         if Eps < 0 then Eps := 0
      end;
      Etot [i] := Epd + Eps;
{
writeln ('*** Point intensity (Etot) ... ', Etot [i]);
writeln;
}
   end
end;


procedure kill_face_aways;

var
   i: integer;
   pos_vec2: vector;

begin
   for i := 1 to numpolys do
   begin
      sub (p [polys [i, 2]], viewer, pos_vec2);
      vis [i] := (dot_prod (N [i], pos_vec2) > 0)
   end
end;


procedure find_min_zs;

var
   i, j: integer;

begin
   for i := 1 to numpolys do
   begin
      minz [i] := 1000;
      for j := 1 to 4 do
         if p [polys [i, j], 3] < minz [i] then minz [i] := p [polys [i, j], 3]
   end
end;


procedure draw_polys;

var
   i, pmin: integer;
   min: real;
   active: boolean;

begin
   repeat
      active := false;
      min := 1000;
      for i := 1 to numpolys do
         if (minz [i] < min) and vis [i] then
         begin
            min := minz [i];
            pmin := i;
            active := true
         end;
{
      if active then
         for i := 1 to numpolys do
         begin
            if minz [i] = min then
            begin
               draw_poly (i);
               vis [i] := false
            end
         end
}
      if active then
      begin
         draw_poly (pmin);
         vis [pmin] := false
      end
   until not active
end;


begin
   write ('Scale');
   readln (scale);

   writeln ('Light coords (x, y, z):');
   readln (light [1]);
   readln (light [2]);
   readln (light [3]);
   light [4] := 1;

{
   writeln ('Enter: Rn, Ws, Id, Ips, k, Rp');
}
   reset (parms, 'parms');
   readln (parms, Rn);
   readln (parms, Ws);
   readln (parms, Id);
   readln (parms, Ips);
   readln (parms, k);
   readln (parms, Rp);

   writeln ('Reading data ...');
   readfile;
   writeln ('Finished data.');

   viewer [1] := 0;
   viewer [2] := 0;
   viewer [3] := -d;
   viewer [4] := 1;

   xa := xa/360*2*pi;
   ya := ya/360*2*pi;
   za := za/360*2*pi;

   setline (xrot, 1, 1, 0, 0, 0);
   setline (xrot, 2, 0, cos (xa), -sin (xa), 0);
   setline (xrot, 3, 0, sin (xa), cos (xa), 0);
   setline (xrot, 4, 0, 0, 0, 1);

   setline (yrot, 1, cos (ya), 0, sin (ya), 0);
   setline (yrot, 2, 0, 1, 0, 0);
   setline (yrot, 3, -sin (ya), 0, cos (ya), 0);
   setline (yrot, 4, 0, 0, 0, 1);

   setline (zrot, 1, cos (za), -sin (za), 0, 0);
   setline (zrot, 2, sin (za), cos (za), 0, 0);
   setline (zrot, 3, 0, 0, 1, 0);
   setline (zrot, 4, 0, 0, 0, 1);

   setline (Tper, 1, 1, 0, 0, 0);
   setline (Tper, 2, 0, 1, 0, 0);
   setline (Tper, 3, 0, 0, 0, 0);
   setline (Tper, 4, 0, 0, -1/d, 1);

   setline (Tscale, 1, scale, 0, 0, 0);
   setline (Tscale, 2, 0, scale, 0, 0);
   setline (Tscale, 3, 0, 0, scale, 0);
   setline (Tscale, 4, 0, 0, 0, 1);

   MMmult (yrot, xrot, mtemp);
   MMmult (zrot, mtemp, Trot);

   writeln ('Rotating points ...');
   for i := 1 to numpoints do
      MVmult (Trot, p [i], p [i]);

   writeln ('Finding minimum z coords for each poly ...');
   find_min_zs;

   writeln ('Finding surface normals ...');
   find_surface_normals;

{
   writeln ('Removing hidden surfaces ...');
   kill_face_aways;
}

   writeln ('Finding point intensities ...');
   find_point_intensities;

   writeln ('Transforming points ...');
   for i := 1 to numpoints do
   begin
      MVmult (Tper, p [i], p [i]);
      norm (p [i]);
      MVmult (Tscale, p [i], p [i])
   end;
   
   writeln ('Drawing polygons ...');
   gpinit;
   last:=0;
   draw_polys;

   writeln ('Finished.')

end.