t=0;
dt = 1;
dx= 10e-3;

lX=150e-3;       %% Comprimento e profundidade da placa
lY=150e-3;
X = fix(lX/dx)
Y = fix(lY/dx)

if X*Y> 80000 fprintf(1,'Perdão, vetor de temperatura grande demais.\n');
return; end; 

[mX,mY] = meshgrid(0:dx:lX,0:dx,lY);
mX= 0:dx:lX;
mY= 0:dx:lY;

%% Constantes do modelo
espessura=10e-3;
cp=100;		% Calor Específico [J/kg*K]
ro=100;		% Densidade [kg/m^3]
k=50.0;		% Condutividade [J/s*m*K]
kk=k*espessura/(ro*cp); % Dispersividade térmica bidimensional [m^2/s]

l=1e-9;         % Ajuste da precisão da gaussiana
xSn=sqrt(4*kk*dt*log(1/l));
xS=floor((xSn/dx)/2)
xS=xS*2+1;

if xS^2 > X*Y fprintf(1,'A reslução espectral está baixa demais.\n');
return; end;

m = (xS-1)/2;
[Xg, Yg] = meshgrid(1:xS,1:xS);
Xg = Xg- mean(Xg(1,:));
Yg = Yg- mean(Yg(:,1));
Xg = Xg .* dx;
Yg = Yg .* dx;

S = exp(-(Xg.^2 + Yg.^2)/(4*kk*dt));

S = S/sum(sum(S));
figure(2); mesh(S);

u = zeros(X,Y);
pause
figure(1);
abac=0
while t < 100
	abac=abac+1;

%%% Regulação da temperatura nas bordas
	if t < 40 tref = 20;
	else tref = 0;
	end;
	u(X,:) = tref; u(1,:) = tref; u(:,1) = tref; u(:,Y) = tref;

	t = t+dt;

%%% Convolução pela gaussiana
	u = conv2(u,S);
	%%% Soma pra dentro o que saiu pra fora (reflexo)
	for a = 1:m
		u(m+a,:) = u(m+a,:)+   u( m+1-a   ,:);
		u(:,m+a) = u(:,m+a)+   u(:,   m+1-a );
		u(m+X+1-a,:) = u(m+X+1-a,:)+   u( m+X+a   ,:);
		u(:,m+Y+1-a) = u(:,m+Y+1-a)+   u(:,   m+Y+a );
	end
	%%% Esquece o resto que veio da convolução
	u = u( (m+1:X+m), (m+1:Y+m));

%%% Potagens
	min=floor(t/60);
	um=sum(sum(u))/(X*Y);	%%% Temperatura média
	fprintf(1,'[%02.0f:%05.2f]: t=%5.2f\r',min,t-60*min,um) ;

	mesh(u');
	title(sprintf('[%02.0f:%05.2f]: t=%5.2f',min,t-60*min,um)) ;
	caxis([0 20]);
	axis([0 X 0 Y 0 20]);

%%% Criação de imagens para a animação
%    eval(cat(2,'print -f1 -r74 -dpng ',sprintf('calor%03.0f',abac ) )) ;
 
	pause(1e-1);
end