program FRACTAL;
{$E+}{$N+}
uses crt,grhacks,h3,r3,r4,r4x4,graph;

const  xmin=-100;
       xmax=100;
       ymin=-100;
       ymax=100;
       maxbrilho=255;
       esq=-40;
       dir=40;
       sup=40;
       inf=-40;
       alpha=0.75;
       nivelmax=7;

type  ponto_do_fractal=record
                          t:real;
                          ptg:h3.point;
                          cor:r3.T;
                        end;
      ponto_de_luz=record
                      coord:h3.point;
                     brilho:r3.T
                   end;

var n,L,v,contx,conty:word;
    lin,col:integer;
    x,y,RE,rz:real;
    cor:r3.T;
    anota:array[1..500] of ponto_do_fractal;
    PL:array[1..10] of ponto_de_luz;
    IM:array[1..4] of r4x4.T;
    id:R4X4.T;
    f,o,u,CE,cz:h3.point;
    mtp:h3.pmap;    {matriz de transf. de perspectiva}
    ch,nome:string;
    arq:file of char;
    cor_nivel:array[0..20] of r3.T;

(*******************************************************************)

function quad(x:double): double;
begin
  quad:=x*x
end;

(*******************************************************************)

function min(x,y:double): double;
begin
  if x<y then min:=x else min:=y
end;

(*******************************************************************)

procedure obtem_coord_cartesianas(p:h3.point);

var i:byte;

begin
  for i:=0 to 3 do p.c[i]:=p.c[i]/p.c[0]
end;

(*******************************************************************)

procedure cria_arquivo;

var c:char;

begin
  writeln ; write('Entre com nome do arquivo que sera gerado: ');
  readln(nome);
  assign(arq, nome);
  rewrite(arq);
  c:='P'; write(arq, c);
  c:='3'; write(arq, c);
  c:=' '; write(arq, c);
  c:='1'; write(arq, c);
  c:='9'; write(arq, c);
  c:='1'; write(arq, c);
  c:=' '; write(arq, c);
  c:='1'; write(arq, c);
  c:='9'; write(arq, c);
  c:='1'; write(arq, c);
  c:=' '; write(arq, c);
  c:='2'; write(arq, c);
  c:='5'; write(arq, c);
  c:='5'; write(arq, c);
  c:=' '; write(arq, c);
  c:=chr(13); write(arq,c);
  c:=chr(10); write(arq,c);
end;

(*******************************************************************)

procedure grava_o_byte(b:byte);

var b1,b2,b3,b4:byte; c,c1,c2,c3:char;

begin
  b4:=b; c:=' ';
  b3:=b4 mod 10; c3:=chr(b3+48); b4:=b4 div 10;
  b2:=b4 mod 10; c2:=chr(b2+48); b4:=b4 div 10;
  b1:=b4;        c1:=chr(b1+48);
  if b1<>0
  then write(arq, c1, c2, c3, c)
  else if b2<>0
       then write(arq, c2, c3, c)
       else write(arq, c3, c)
end;

(*******************************************************************)

procedure grava_cor_do_pixel(cor:r3.T);

var r,g,b:byte; c:char;

begin
  r:=round(maxbrilho*cor[0]); grava_o_byte(r);
  g:=round(maxbrilho*cor[1]); grava_o_byte(g);
  b:=round(maxbrilho*cor[2]); grava_o_byte(b);
  c:=chr(13); write(arq,c);
  c:=chr(10); write(arq,c)
end;

(*******************************************************************)

procedure constroi_IM(alpha, beta, teta, fi: real; var res:R4X4.T);

var INVM:array[1..5] of R4X4.T;

begin
{  M[1,0,0]:=1.0; M[1,0,1]:=0.0;        M[1,0,2]:=0.0;     M[1,0,3]:=0.0;
  M[1,1,0]:=0.0; M[1,1,1]:=alpha;      M[1,1,2]:=0.0;     M[1,1,3]:=0.0;
  M[1,2,0]:=0.0; M[1,2,1]:=0.0;        M[1,2,2]:=alpha;   M[1,2,3]:=0.0;
  M[1,3,0]:=0.0; M[1,3,1]:=0.0;        M[1,3,2]:=0.0;     M[1,3,3]:=alpha;

  M[2,0,0]:=1.0; M[2,0,1]:=0.0;        M[2,0,2]:=0.0;     M[2,0,3]:=rz*(1+beta);
  M[2,1,0]:=0.0; M[2,1,1]:=1.0;        M[2,1,2]:=0.0;     M[2,1,3]:=0.0;
  M[2,2,0]:=0.0; M[2,2,1]:=0.0;        M[2,2,2]:=1.0;     M[2,2,3]:=0.0;
  M[2,3,0]:=0.0; M[2,3,1]:=0.0;        M[2,3,2]:=0.0;     M[2,3,3]:=1.0;

  M[3,0,0]:=1.0; M[3,0,1]:=0.0;        M[3,0,2]:=0.0;       M[3,0,3]:=0.0;
  M[3,1,0]:=0.0; M[3,1,1]:=cos(teta);  M[3,1,2]:=sin(teta); M[3,1,3]:=0.0;
  M[3,2,0]:=0.0; M[3,2,1]:=-sin(teta); M[3,2,2]:=cos(teta); M[3,2,3]:=0.0;
  M[3,3,0]:=0.0; M[3,3,1]:=0.0;        M[3,3,2]:=0.0;       M[3,3,3]:=1.0;

  M[4,0,0]:=1.0; M[4,0,1]:=0.0;        M[4,0,2]:=0.0;       M[4,0,3]:=0.0;
  M[4,1,0]:=0.0; M[4,1,1]:=1.0;        M[4,1,2]:=0.0;       M[4,1,3]:=0.0;
  M[4,2,0]:=0.0; M[4,2,1]:=0.0;        M[4,2,2]:=cos(fi);   M[4,2,3]:=sin(fi);
  M[4,3,0]:=0.0; M[4,3,1]:=0.0;        M[4,3,2]:=-sin(fi);  M[4,3,3]:=cos(fi);
 }
  INVM[1,0,0]:=1.0; INVM[1,0,1]:=0.0;          INVM[1,0,2]:=0.0;       INVM[1,0,3]:=0.0;
  INVM[1,1,0]:=0.0; INVM[1,1,1]:=1/alpha;      INVM[1,1,2]:=0.0;       INVM[1,1,3]:=0.0;
  INVM[1,2,0]:=0.0; INVM[1,2,1]:=0.0;          INVM[1,2,2]:=1/alpha;   INVM[1,2,3]:=0.0;
  INVM[1,3,0]:=0.0; INVM[1,3,1]:=0.0;          INVM[1,3,2]:=0.0;       INVM[1,3,3]:=1/alpha;

  INVM[2,0,0]:=1.0; INVM[2,0,1]:=0.0;        INVM[2,0,2]:=0.0;     INVM[2,0,3]:=-rz*(1+beta);
  INVM[2,1,0]:=0.0; INVM[2,1,1]:=1.0;        INVM[2,1,2]:=0.0;     INVM[2,1,3]:=0.0;
  INVM[2,2,0]:=0.0; INVM[2,2,1]:=0.0;        INVM[2,2,2]:=1.0;     INVM[2,2,3]:=0.0;
  INVM[2,3,0]:=0.0; INVM[2,3,1]:=0.0;        INVM[2,3,2]:=0.0;     INVM[2,3,3]:=1.0;

  INVM[3,0,0]:=1.0; INVM[3,0,1]:=0.0;        INVM[3,0,2]:=0.0;       INVM[3,0,3]:=0.0;
  INVM[3,1,0]:=0.0; INVM[3,1,1]:=cos(teta);  INVM[3,1,2]:=-sin(teta); INVM[3,1,3]:=0.0;
  INVM[3,2,0]:=0.0; INVM[3,2,1]:=sin(teta);  INVM[3,2,2]:=cos(teta); INVM[3,2,3]:=0.0;
  INVM[3,3,0]:=0.0; INVM[3,3,1]:=0.0;        INVM[3,3,2]:=0.0;       INVM[3,3,3]:=1.0;

  INVM[4,0,0]:=1.0; INVM[4,0,1]:=0.0;        INVM[4,0,2]:=0.0;       INVM[4,0,3]:=0.0;
  INVM[4,1,0]:=0.0; INVM[4,1,1]:=1.0;        INVM[4,1,2]:=0.0;       INVM[4,1,3]:=0.0;
  INVM[4,2,0]:=0.0; INVM[4,2,1]:=0.0;        INVM[4,2,2]:=cos(fi);   INVM[4,2,3]:=-sin(fi);
  INVM[4,3,0]:=0.0; INVM[4,3,1]:=0.0;        INVM[4,3,2]:=sin(fi);  INVM[4,3,3]:=cos(fi);

  R4X4.mul(INVM[4], INVM[3], INVM[5]);
  R4X4.mul(INVM[5], INVM[2], INVM[5]);
  R4X4.mul(INVM[5], INVM[1], res)
end;

(*******************************************************************)

procedure le_os_dados;

const alpha=0.75;
      beta =0.0;
      teta =pi/2;
      fi   =pi/6;

var x1, y1, z1, x2, y2, z2, x3, y3, z3, r:double;
    i:word;

begin
  h3.mk_point(1.0,  0.0, 0.0, 0.0, CE); RE:=40;
  h3.mk_point(1.0,  0.0, 0.0, 0.0, cz); rz:=10;

  n:=2;

  constroi_IM(alpha, beta, teta,  fi, IM[1]);
  constroi_IM(alpha, beta, teta, -fi, IM[2]);


  id[0,0]:=1.0; id[0,1]:= 0.0;   id[0,2]:=0.0; id[0,3]:=0.0;
  id[1,0]:=0.0; id[1,1]:= 1.0;   id[1,2]:=0.0; id[1,3]:=0.0;
  id[2,0]:=0.0; id[2,1]:= 0.0;   id[2,2]:=1.0; id[2,3]:=0.0;
  id[3,0]:=0.0; id[3,1]:= 0.0;   id[3,2]:=0.0; id[3,3]:=1.0;


  L:=2;

  PL[1].coord.c[0]:=1.0; PL[1].coord.c[1]:=50;
  PL[1].coord.c[2]:=20;  PL[1].coord.c[3]:=50;
  PL[1].brilho[0]:=0.8;  PL[1].brilho[1]:=0.8; PL[1].brilho[2]:=0.8;

  PL[2].coord.c[0]:=1.0; PL[2].coord.c[1]:=50;
  PL[2].coord.c[2]:=-20;  PL[2].coord.c[3]:=50;
  PL[2].brilho[0]:=0.4;  PL[2].brilho[1]:=0.2; PL[2].brilho[2]:=0.3;

  for i:=0 to nivelmax do
  begin
    r:=i/nivelmax;
    cor_nivel[i][0]:=r;
    cor_nivel[i][1]:=1-r;
    cor_nivel[i][2]:=0.5
  end;
  writeln('Entre com obs, foco e u.');
  readln(x1, y1, z1, x2, y2, z2, x3, y3, z3);
  h3.mk_point(1.0, x1, y1, z1, o);
  h3.mk_point(1.0, x2, y2, z2, f);
  h3.mk_point(1.0, x3, y3, z3, u)

end;

(*******************************************************************)

function muito_proximos(p,q:h3.point) : boolean;

begin
  if (abs(p.c[1] - q.c[1]) < 0.01) and
     (abs(p.c[2] - q.c[2]) < 0.01) and
     (abs(p.c[3] - q.c[3]) < 0.01)
  then muito_proximos:=true
  else muito_proximos:=false
end;

(*******************************************************************)

procedure anota_ponto(x:real; ptg:h3.point; cor:r3.T);
begin
  v:=v+1; anota[v].t:=x;
  anota[v].ptg:=ptg;
  anota[v].cor:=cor;
end;

(*******************************************************************)

function minimo_da_lista : word;

var  i,m:word; d,d1:double;

begin
  if v=0
  then minimo_da_lista:=0
  else begin
         m:=1;
         for i:=2 to v do if anota[i].t < anota[m].t then m:=i;
         minimo_da_lista:=m
       end
end;

(*******************************************************************)
{
procedure vetor_unitario(a,b:h3.point; var res:h3.point);

var c:double; i:word;

begin
  res.c[0]:=1.0;
  for i:=1 to 3 do res.c[i]:=b.c[i] - a.c[i];
  c:=sqrt(1.0/(quad(res.c[1]) + quad(res.c[2]) + quad(res.c[3])));
  for i:=1 to 3 do res.c[i]:=c*res.c[i];
end;
}
(*******************************************************************)

procedure obtem_ponto_eq_parametrica(p1, p2:h3.point; t:real;
                                    var res:h3.point);
var i:word;

begin
  res.c[0]:=1.0;
  for i:=1 to 3 do res.c[i]:=p1.c[i] + t*(p2.c[i] - p1.c[i])
end;

(*******************************************************************)

procedure mapeia( pto:h3.point; raio:real; mat:R4X4.T; var res:h3.point);

var mat1:R4X4.T;
    i:byte;

begin
  mat1:=mat;
  for i:=2 to 3 do mat1[0,i]:=mat1[0,i]*raio;
  R4X4.map_row(pto.c, mat1, res.c);
end;

(*******************************************************************)

procedure resolve_eq_2_gr(a, b, c:double; var nsol:byte; var x1, x2:double);

var delta:double;

begin
  delta:=b*b - 4*a*c;
  if delta < 0
  then nsol:=0
  else if delta = 0
       then begin
              nsol:=1 ; x1 := (-b)/(2*a)
            end
       else begin
              nsol:=2 ;
              x1 := ((-b) + sqrt(delta))/(2*a) ;
              x2 := ((-b) - sqrt(delta))/(2*a)
            end
end;

(*******************************************************************)

function intercepta( c:h3.point; r:real; o, p:h3.point; var t:real): boolean ;

var a,b,d,x1,x2:double;  nsol:byte;

begin
  a := quad(p.c[1]-o.c[1]) + quad(p.c[2]-o.c[2]) + quad(p.c[3]-o.c[3]) ;
  b := 2*((p.c[1]-o.c[1])*(o.c[1]-c.c[1]) +
          (p.c[2]-o.c[2])*(o.c[2]-c.c[2]) +
          (p.c[3]-o.c[3])*(o.c[3]-c.c[3]));
  d := quad(o.c[1]-c.c[1]) + quad(o.c[2]-c.c[2]) + quad(o.c[3]-c.c[3]) - r*r;

  resolve_eq_2_gr(a, b, d, nsol, x1, x2);

  if (nsol=0) or ((nsol=1) and (x1<=0)) or ((nsol=2) and (x1<=0) and (x2<=0))
  then intercepta:=false
  else begin
         intercepta:=true;
         if nsol=1
         then  t:=x1
         else if x1 <= 0
              then t:=x2
              else if x2<=0
                   then t:=x1
                   else t:=min(x1,x2)
       end
end;

(*******************************************************************)

procedure varre_o_fractal(o, p:h3.point; imat:R4X4.T; nivel:integer);

var CE1, cz1, om, pm, p1, ptg:h3.point;
    imat1:R4X4.T;
    i:word;
    t:real;

begin
  R4X4.map_row(o.c, imat, om.c);
  R4X4.map_row(p.c, imat, pm.c);
  if intercepta(CE, RE, om, pm, t)
  then
       begin
         if intercepta(cz, rz, om, pm, t)
         then begin
                obtem_ponto_eq_parametrica(om, pm, t, p1);
                p1.c[0]:=-p1.c[0]*rz*rz;
                R4X4.map_col(imat, p1.c, ptg.c);
                anota_ponto(t, ptg, cor_nivel[nivel]);
              end;
         if nivel > 0
         then begin
                for i:=1 to n do
                  begin
                    R4X4.mul(imat, IM[i], imat1);
                    varre_o_fractal(o, p, imat1, nivel-1)
                  end
              end
       end
end;

(*******************************************************************)

function enxerga(org, dest:h3.point): boolean;

var q:word; p1:h3.point;

begin
  v:=0;
  varre_o_fractal(org, dest, id, nivelmax);
  q:=minimo_da_lista;
  if q=0 then enxerga:=true
         else begin
                obtem_ponto_eq_parametrica(org, dest, anota[q].t, p1);
                if muito_proximos(p1, dest) then enxerga:=true
                                            else enxerga:=false
              end
end;

(*******************************************************************)

procedure calcula_cor_do_ponto(q, ptg:h3.point; cor_fractal:r3.T; var res:r3.t);

var cor_prov,cori,f,n:r3.T; fator:double; i:word;
    f1,n1:h3.point;

begin
  cor_prov[0]:=0.0; cor_prov[1]:=0.0; cor_prov[2]:=0.0;
  n[0]:=ptg.c[1]; n[1]:=ptg.c[2]; n[2]:=ptg.c[3];
  r3.reduce(n);
  for i:=1 to L do
   if enxerga(PL[i].coord, q)
   then begin
          h3.dir(q, PL[i].coord, f);
          fator:=r3.dot(f, n);
          if fator > 0.0 then begin
                                r3.scale(PL[i].brilho, fator, cori);
                                r3.add(cor_prov, cori, cor_prov)
                              end
       end;
  res[0]:=cor_prov[0]*cor_fractal[0];
  res[1]:=cor_prov[1]*cor_fractal[1];
  res[2]:=cor_prov[2]*cor_fractal[2];
end;

(*******************************************************************)

procedure calcula_cor_do_pixel(x, y:real; var res:r3.T);

var p,res1,p1:h3.point; q:word; t:real;

begin
  h3.mk_point(1.0, x, y, 0.0, p);
  r4x4.map_row(p.c, mtp.inv, p.c);
  obtem_coord_cartesianas(p);
  v:=0;
  varre_o_fractal(o, p, id, nivelmax);
  q:=minimo_da_lista;
  if q=0 then begin
                res[0]:=0.0; res[1]:=0.0; res[2]:=0.0
              end
         else begin
                obtem_ponto_eq_parametrica(o, p, anota[q].t, p1);
                calcula_cor_do_ponto(p1, anota[q].ptg, anota[q].cor, res)
              end

{
  if intercepta(CE, PE, o, p, t)
  then begin
         v:=0; anota_ponto(t, CE);
         obtem_ponto_eq_parametrica(o, p, anota[1].t, p1);
         calcula_cor_do_ponto(p1, anota[1].c, res)
       end
  else begin
         res[0]:=0.0; res[1]:=0.0; res[2]:=0.0
       end
}
end;

(*******************************************************************)

procedure plota_o_pixel(x, y:real; var cor:r3.T);

var i:word;  c:real;

begin
  for i:=0 to 2 do if cor[i]<0.0 then cor[i]:=0.0;
  for i:=0 to 2 do if cor[i]>1.0 then cor[i]:=1.0;
  c:=0.0;
  for i:=0 to 2 do c:=c + cor[i];
  c:=c/3;
  grhacks.gr_set_color(c);
  grhacks.gr_draw_point(x,y);
end;

(*******************************************************************)
(*******************************************************************)

begin
  le_os_dados;
  h3.persp_map(o, f, u, mtp);

  cria_arquivo;

  grhacks.gr_begin_drawing(xmin, xmax, ymin, ymax);
  grhacks.gr_set_color(0.5);
  grhacks.gr_paint_plane;

  y:=sup; lin:=grhacks.gr_y(y); conty:=0;
  while y >= inf
  do begin

{       if conty=0 then begin
 }
       x:=esq; col:=grhacks.gr_x(x); contx:=0;
       while x <= dir
       do begin
  {          contx:=contx+1;
   }
            calcula_cor_do_pixel(x,y,cor);
            plota_o_pixel(x,y,cor);
            grava_cor_do_pixel(cor);

            x:=x+0.4;
            while (x<=dir) and (col=grhacks.gr_x(x)) do x:=x+0.4;
            col:=grhacks.gr_x(x)
          end;
{                                 end;
      conty:=conty + 1;}

       y:=y-0.4;
       while (y>=inf) and (lin=grhacks.gr_y(y)) do y:=y-0.4;
       lin:=grhacks.gr_y(y)
     end;

 closegraph;
 close(arq);
{ writeln('linhas:',conty,'  colunas:',contx);

  read(ch);}
end.