t=0;
f=0.01;
%f=440;
dt= 1/(f*8);
%dt=1;
dta = 1/ (8*2*f);
dx= 2e-3;

lX=15e-2;	%% Comprimento e profundidade da placa
lY=7.5e-2;
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=890;		% Calor Específico [J/kg*K]
ro=2700;	% Densidade [kg/m^3]
k=225.0;	% Condutividade [J/s*m*K]
kk=k*espessura/(ro*cp); % Dispersividade térmica bidimensional [m^2/s]
%h=60;
h=1000;

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;

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 1
%while t<1200
	abac=abac+1;

%%% Cálculo da tensão, corrente e potências nos transistores
	ta = [t:dta:t+dt];
	V = 20*sin(2*pi*ta*f);
	I = V/8;
	P1 = sum((V>0) .* (+25-V).*I)/length(ta);
	P2 = sum((V<0) .* (-25-V).*I)/length(ta);

	t = t+dt;

%%% Calor entrando
	Dx=find((mX>=0.10).*(mX<=0.11));
	Dy=find((mY>=0.02).*(mY<=0.03));
	area = length(Dx)*length(Dy);
	u(Dx,Dy)=u(Dx,Dy)+  dt*(P1/area)/(ro*cp*espessura*dx^2);
	Dx=find((mX>=0.05).*(mX<=0.06));
	Dy=find((mY>=0.02).*(mY<=0.03));
	area = length(Dx)*length(Dy);
	u(Dx,Dy)=u(Dx,Dy)+  dt*(P2/area)/(ro*cp*espessura*dx^2);

%%% 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));

%%% Convecção
	uv = u;
	u = u .* exp(-h*(dx^2)*dt);

%%% Potagens
	qconv=cp*ro*espessura*dx^2*sum(sum(u-uv))/(dt);

	[umin,gr]=min(u-1);[umin,gr]=min(umin);

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

	mesh(u');
	title(sprintf('[%02.0f:%05.2f]: t=%5.2f qconv=%05.2f qtrans=%05.2f',min,t-60*min,um, qconv, P1+P2   )) ;
	caxis([umin*0.8 100]);
	axis([0 X 0 Y umin*0.85 100]);
%%% Para gerar imagens
%		eval(cat(2,'print -f1 -r74 -dpng ',sprintf('ampAB%03.0f',abac ) )) ;
	pause(1e-1);
end