{$A+,B-,D+,E+,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S-,T-,V+,X+}
{$M 16384,0,655360}

{ Geracao de Texturas em uma Dimensao por Semelhanca }

{Pressione ESC a qualquer momento para terminar}

{ nova versao com graficos de pares}
{ nova versao com derivadas temporais normalizadas,
  permitindo melhor controle com o fator dt
  que determina a derivada temporal maxima
}
{ nova versao com vetores mostrando a evolucao }

program semelhan;
uses Graph,Crt,H3;

type real=double;  {usa double ao inves de real}

const
  tam_amostra=50;
  tam_entrada=100;
  raio_distr=0.1;
  dt_default=0.01/tam_entrada;

type
  amostra=array[0..tam_amostra-1]of real;
  entrada=array[0..tam_entrada-1]of real;
  generic=array[0..(64000 div sizeof(real))]of real;

{------------- Gera valores iniciais -------------}
procedure GeraAmostraQualquer(var a:amostra);
var i:integer;
begin
  for i:=0 to tam_amostra-1 do
    {a[i]:=-3*(exp(-sqr(i-tam_amostra div 2)/8)+
                    exp(-sqr(i-6-tam_amostra div 2)/8));}
    a[i]:=3.0*sin(2.0*Pi*i/tam_amostra);
    {a[i]:=3.0*((tam_amostra/2)-i)/tam_amostra;}
end;

procedure GeraEntradaAleatoria(var a:entrada);
const
  menor=-4.0;
  maior=+4.0;
var i:integer;
begin
  for i:=0 to tam_entrada-1 do
     {a[i]:=sin(2.0*Pi*i/tam_amostra);}
     {a[i]:=menor+(maior-menor)*random;}
     a[i]:=0.5*(sin(5.0*2.0*Pi*i/tam_entrada)+sin(4.0*Pi*i/tam_entrada));
end;

{-------------- Funcoes Matematicas -------------}
function pow3_2(x:real):real;
begin
  pow3_2:=x*sqrt(x);
end;


{------------- Calculo de Invariantes -----------}
procedure GeraInvariantes(var F; var U,V,W; size:integer);
var i,iprev,inext:integer;
    af:generic absolute F;
    au:generic absolute U;
    av:generic absolute V;
    aw:generic absolute W;
begin
  for i:=0 to size-1 do begin
     iprev:=(i-1+size) mod size;
     inext:=(i+1) mod size;
     au[i]:=af[i];
     av[i]:=af[iprev];
     aw[i]:=af[inext];
  end;
end;

{----------- Parte Grafica ---------------}
type
   janela=record
       x0,y0,vx,vy:integer;
       title:string;
       back,fore,border:integer;
   end;

procedure FazJanela(var J:janela; Atitle:string; Ax0,Ay0,
                      Avx,Avy,Aback,Afore,Aborder:integer);
begin
  with J do begin
    x0:=Ax0; y0:=Ay0; vx:=Avx; vy:=Avy;
    title:=Atitle;
    back:=Aback; fore:=Afore; border:=Aborder;
  end;
end;

procedure PoeJanela(var J:janela);
begin
  with J do begin
    SetFillStyle(SolidFill,back);
    Bar(x0-1,y0-10,x0+vx+1,y0+vy+1);
    SetColor(border);
    Rectangle(x0-1,y0-10,x0+vx+1,y0+vy+1);
    Line(x0-1,y0-1,x0+vx+1,y0-1);
    SetTextJustify(CenterText,CenterText);
    SetColor(fore);
    OutTextXY(x0+ vx div 2, y0-5, title);
  end;
end;

procedure TiraJanela(var J:janela; cor:integer);
begin
  with J do begin
    SetFillStyle(SolidFill,cor);
    Bar(x0-1,y0-10,x0+vx+1,y0+vy+1);
  end;
end;

procedure LimpaJanela(var J:janela);
begin
  with J do begin
    SetFillStyle(SolidFill,back);
    Bar(x0-1,y0-1,x0+vx+1,y0+vy+1);
    SetColor(border);
    Rectangle(x0-1,y0-10,x0+vx+1,y0+vy+1);
    Line(x0-1,y0-1,x0+vx+1,y0-1);
  end;
end;

procedure PutPixelJanela(var J:janela; x,y,cor:integer);
begin
  with J do begin
    if x>vx then exit;
    if y>vy then exit;
    if x<0 then exit;
    if y<0 then exit;
    PutPixel(x+x0,y+y0,cor);
  end;
end;


procedure PlotaPares(var J:janela; var aa,bb; size:integer;color:integer);
const escala_a=10.0;
      escala_b=40.0;
var  a:generic absolute aa;
     b:generic absolute bb;
     x0,y0,i:integer;

begin
  for i:=0 to size-1 do begin
    x0:=J.vx div 2 +round(escala_a*a[i]);
    y0:=J.vy div 2 +round(escala_b*b[i]);
    PutPixelJanela(J,x0,y0,color);
  end;
end;

Procedure PlotaTriplas(var J:janela; var pm:H3.PMap;
                       var xx,yy,zz; size:integer; color:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    p:H3.Point;
    i:integer;
begin
  for i:=0 to size-1 do begin
    p.c[0]:=1.0; p.c[1]:=x[i]*escala_x; p.c[2]:=y[i]*escala_y; p.c[3]:=z[i]*escala_z;
    map_point(p,pm,p);
    PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color);
  end;
end;

Procedure PlotaTriplasDistrib(var J:janela; var pm:H3.PMap;
                       var xx,yy,zz; size:integer; color:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    p,q:H3.Point;
    i,k:integer;
begin
  for i:=0 to size-1 do begin
    for k:=1 to 10 do begin
      p.c[0]:=1.0;
      p.c[1]:=(x[i]+raio_distr*sin(3*k))*escala_x;
      p.c[2]:=(y[i]+raio_distr*cos(5*k))*escala_y;
      p.c[3]:=(z[i]+raio_distr*sin(7*k))*escala_z;
      map_point(p,pm,p);
      PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color);
    end;
  end;
end;

Procedure PlotaTriplasVetores(var J:janela; var pm:H3.PMap;
                              var oxx,oyy,ozz,xx,yy,zz;
                              size:integer; color0,color,color1:integer);
const escala_x=40.0; escala_y=200.0; escala_z=40.0;
      escala_df=10.0;
var x:generic absolute xx;
    y:generic absolute yy;
    z:generic absolute zz;
    ox:generic absolute oxx;
    oy:generic absolute oyy;
    oz:generic absolute ozz;
    p,op:H3.Point;
    i:integer;
begin
  for i:=0 to size-1 do begin
    p.c[0]:=1.0; p.c[1]:=x[i]*escala_x; p.c[2]:=y[i]*escala_y; p.c[3]:=z[i]*escala_z;
    map_point(p,pm,p);
    op.c[0]:=1.0;
    op.c[1]:=(ox[i]+escala_df*(x[i]-ox[i]))*escala_x;
    op.c[2]:=(oy[i]+escala_df*(y[i]-oy[i]))*escala_y;
    op.c[3]:=(oz[i]+escala_df*(z[i]-oz[i]))*escala_z;
    map_point(op,pm,op);
    SetColor(color);
    Line(J.x0+ J.vx div 2 + round(p.c[1]),J.y0+ J.vy div 2 -round(p.c[2]),
         J.x0+ J.vx div 2 + round(op.c[1]),J.y0+ J.vy div 2 -round(op.c[2]));
    PutPixelJanela(J, J.vx div 2 + round (p.c[1]),
                   J.vy div 2 -round(p.c[2]), color1);
    PutPixelJanela(J, J.vx div 2 + round (op.c[1]),
                   J.vy div 2 -round(op.c[2]), color0);
  end;
end;


procedure TracaDados(var f; size:integer; x,y:integer); {sem janela}
const escala=10.0; {1:10}
var i:integer;
    a:generic absolute f;
begin
  for i:=0 to size-1 do putpixel(x+i,y,red);
  moveto(x,y+round(escala*a[0]));
  for i:=1 to size-1 do lineto(x+i,y+round(escala*a[i]));
end;

procedure PlotaDados(var J:janela;var f; size:integer);
const escala=10.0; {1:10}
var i:integer;
    a:generic absolute f;
begin
  for i:=0 to size-1 do putpixelJanela(J,i,J.vy div 2,red);
  for i:=0 to size-1 do putpixeljanela(J,i,J.vy div 2+round(escala*a[i]),cyan);
end;

function GetString(var J:janela; message,default:string):string;
var i,l,key,stop:integer;
begin
  J.title:=message;
  message:=default;
  l:=length(default);
  i:=l+1;
  PoeJanela(J);
  SetTextJustify(LeftText,TopText);
  stop:=0;
  repeat
    LimpaJanela(J);
    SetColor(J.Fore);
    OutTextXy(J.x0,J.y0,default);
    SetFillStyle(SolidFill,J.Fore);
    Bar(J.x0+8*(i-1),J.y0,J.x0+8*i,J.y0+8);
    SetColor(J.Back);
    if i<=l then OutTextXy(J.x0+8*(i-1),J.y0,default[i]);
    key:=ord(readkey);
    if key=0 then key:=256+ord(readkey);
    case key of
      13: begin GetString:= default; stop:=1; end;
      27: begin GetString:= message; stop:=1; end;
      8:  begin default:=copy(default,1,l-1); dec(l); dec(i); end;
      else begin
        if (key>31) and (key<127) then begin
          if stop=0 then begin default:=chr(key); stop:=2; i:=2; l:=1;end
          else begin default:=default+chr(key); inc(i); inc(l); end;
        end;
      end;
    end;
  until stop=1;
end;

function GetReal(var J:janela; message,default:string):real;
var r:real; s:string;
    i,cod:integer;
begin
  repeat
    s:=GetString(J,message,default);
    val(s,r,cod);
    if cod<>0 then write(#7);
  until cod=0;
  GetReal:=r;
end;

{///////// Parte Principal do Programa /////////}

var af,au,av,aw:amostra;
    f,u,v,w:entrada;
    t,key:integer;
    obs,foc,upp:H3.Point;
    pm:H3.PMap;

    GD,GM:integer;

procedure AskPMap;
var J:janela;
begin
  FazJanela(J,'',GetMaxX div 2,9*GetMaxY div 10, 200, 8,DarkGray,Green,Red);
  obs.c[0]:=1.0; foc:=obs; upp:=obs;
  obs.c[1]:=GetReal(J,'Entre observador eixo X','3.0');
  obs.c[2]:=GetReal(J,'Entre observador eixo Y','-3.0');
  obs.c[3]:=GetReal(J,'Entre observador eixo Z','3.0');
  foc.c[1]:=GetReal(J,'Entre foco eixo X','0.0');
  foc.c[2]:=GetReal(J,'Entre foco eixo Y','0.0');
  foc.c[3]:=GetReal(J,'Entre foco eixo Z','0.0');
  upp.c[1]:=GetReal(J,'Entre zenith eixo X','0.0');
  upp.c[2]:=GetReal(J,'Entre zenith eixo Y','0.0');
  upp.c[3]:=GetReal(J,'Entre zenith eixo Z','1.0');
  H3.persp_map(obs,foc,upp,pm);
  TiraJanela(J,Black);
end;


procedure Evolve;
const
  alfa:array[-1..1,1..3] of real=
       (( 0,  1   ,  0),
        ( 1,  0   ,  0),
        ( 0,  0   ,  1));
      dt=0.05;
      poda=6.0;
      tol=1.0;
var r,du,dv,dw,max:real;
    dedu,dedv,dedw,dfdt:entrada;
    i,iprev,inext,j:integer;
  function epp(i,j:integer):real;
  begin
    epp:=exp(-0.5*tol*sqr(f[i]-f[j]));
  end;
  function epq(i,j:integer):real;
  begin
    epq:=exp(-0.5*tol*sqr(f[i]-af[j]));
  end;

begin
  GeraInvariantes(f,u,v,w,tam_entrada);
  for i:=0 to tam_entrada-1 do begin
    iprev:= (i-1+tam_entrada) mod tam_entrada;
    inext:= (i+1) mod tam_entrada;
    dfdt[i]:=epp(i,iprev)*(
                     (u[i]-u[iprev]) *(alfa[0,1]-alfa[-1,1])
                     +(v[i]-v[iprev])*(alfa[0,2]-alfa[-1,2])
                     +(w[i]-w[iprev])*(alfa[0,3]-alfa[-1,3]))
             +epp(inext,iprev)*(
                     (u[inext]-u[iprev]) *(alfa[1,1]-alfa[-1,1])
                     +(v[inext]-v[iprev])*(alfa[1,2]-alfa[-1,2])
                     +(w[inext]-w[iprev])*(alfa[1,3]-alfa[-1,3]))
             +epp(inext,i)*(
                     (u[inext]-u[i])   *(alfa[1,1]-alfa[0,1])
                     +(v[inext]-v[i])  *(alfa[1,2]-alfa[0,2])
                     +(w[inext]-w[i])  *(alfa[1,3]-alfa[0,3]));

    for j:=0 to tam_amostra-1 do begin
      dfdt[i]:=dfdt[i]+ (tam_entrada/tam_amostra)*(
             epq(iprev,j)*(
               (u[iprev]-aU[j]) *alfa[-1,1]
               +(v[iprev]-aV[j])*alfa[-1,2]
               +(w[iprev]-aW[j])*alfa[-1,3])
             +epq(i,j)*(
               (u[i]-aU[j]) *alfa[0,1]
               +(v[i]-aV[j])*alfa[0,2]
               +(w[i]-aW[j])*alfa[0,3])+
             epq(inext,j)*(
               (u[inext]-aU[j]) *alfa[1,1]
               +(v[inext]-aV[j])*alfa[1,2]
               +(w[inext]-aW[j])*alfa[1,3])
             );

    end;
  end;

  max:=0.0;
  for i:=0 to tam_entrada-1 do begin
    if abs(dfdt[i])>max then max:=abs(dfdt[i]);
  end;
  for i:=0 to tam_entrada-1 do begin
    f[i]:=f[i]-dfdt[i]*dt/max;
    if f[i]>poda then f[i]:=poda;
    if f[i]<-poda then f[i]:=-poda;
  end;
end;


procedure Evolvesteps(t:integer);
var i:integer;
begin
  for i:=1 to t do Evolve;
end;

var JanelaEntrada,     JanelaAmostra, JanelaDerivada,
    JanelaLaplaciano,  Janela3d  :janela;
    oldf,oldv,oldw:entrada;

begin

  GD:=Detect;
  InitGraph(GD,GM,'');

  AskPMap;

  FazJanela(JanelaDerivada,'Derivadas', 1,GetMaxY div 2,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaLaplaciano,'Laplaciano',GetMaxX div 4,GetMaxY div 2,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaAmostra,'Funcao Inicial',1,10,120,200,
            DarkGray,Green,Red);
  FazJanela(JanelaEntrada,'Fc em Evolucao',GetMaxX div 4,10,120,200,
            DarkGray,Green,Red);
  FazJanela(Janela3d,'Triplas f x f` x f``',GetMaxX div 2,10,
             4*GetMaxX div 10, 7*GetMaxY div 10,
             DarkGray,Green,Red);

  GeraAmostraQualquer(af);
  GeraInvariantes(af,au,av,aw,tam_amostra);
  GeraEntradaAleatoria(f);
  GeraInvariantes(f,u,v,w,tam_entrada);

  PoeJanela(JanelaAmostra);
  PoeJanela(JanelaEntrada);
  PlotaDados(JanelaAmostra,af,tam_amostra);
  PlotaDados(JanelaEntrada,f, tam_entrada);
  PoeJanela(JanelaDerivada);
  Plotapares(JanelaDerivada,af,av,tam_amostra,LightMagenta);
  Plotapares(JanelaDerivada,f,v,tam_entrada,Cyan);
  PoeJanela(JanelaLaplaciano);
  Plotapares(JanelaLaplaciano,af,aw,tam_amostra,LightMagenta);
  Plotapares(JanelaLaplaciano,f,w,tam_entrada,Cyan);
  PoeJanela(Janela3d);
  PlotaTriplas(Janela3d,pm,af,av,aw,tam_amostra,LightMagenta);
  PlotaTriplas(Janela3d,pm, f, v, w,tam_entrada,cyan);

  readkey;
  while keypressed do readkey;
  t:=1;
  repeat
    oldf:=f; oldv:=v; oldw:=w;
    key:=0;
    EvolveSteps(t);
    inc(t);
{    LimpaJanela(JanelaAmostra);                 }
{    PlotaDados(JanelaAmostra,af,tam_amostra);   }
    LimpaJanela(JanelaEntrada);
    PlotaDados(JanelaEntrada,f, tam_entrada);
    LimpaJanela(JanelaDerivada);
    Plotapares(JanelaDerivada,af,av,tam_amostra,LightMagenta);
    Plotapares(JanelaDerivada,f,v,tam_entrada,Cyan);
    LimpaJanela(JanelaLaplaciano);
    Plotapares(JanelaLaplaciano,af,aw,tam_amostra,LightMagenta);
    Plotapares(JanelaLaplaciano,f,w,tam_entrada,Cyan);
    LimpaJanela(Janela3d);
    PlotaTriplas(Janela3d,pm,af,av,aw,tam_amostra,LightMagenta);
    PlotaTriplasVetores(Janela3d,pm,oldf,oldv,oldw,f,v,w,tam_entrada,
      Blue,LightBlue,Cyan);
    if keypressed then begin
      key:=ord(readkey);
      if key=0 then key:=256+ord(readkey);
    end;
  until key=27;


  CloseGraph;
end.