%
% Este programa resuelve numericamente el problema de
% difusion de calor en una barra inicialmente con
% temperatura uniforme, la cual se pone en contacto 
% en ambos extremos con reservorios caloricos
% a temperatura constante.
%
% La ecuacion original es 
% dT/dt=d/dx(Difusividad*dT/dx)
% donde las derivadas son parciales, T es la temperatura,
% t es el tiempo y x es la distancia desde el extremo izquierdo.
%
% Las condiciones iniciales y de frontera son
% T(x,t=0)=Tinicial, T(x=0,t)=Tizquierda, T(x=Longitud,t)=Tderecha
% donde "Longitud" es la longitud de la barra.
%
% El esquema numerico es de diferencias finitas centradas de segundo
% orden en el espacio y Euler adelantado de primer orden en el tiempo.
% O sea
% T(x,t+dt)=T(x,t)+Difusividad*dt/(dx*dx)*(T(x-dx,t)-2*T(x,t)+T(x+dx,t))
% donde dt es el paso de tiempo y dx es la distancia entre puntos de grilla.
%
%
%                                  Ken Takahashi (ktakahashi@geo.igp.gob.pe)
%
clear

Tinicial  =20; % Temperatura inicial de la barra
Tizquierda=10; % Temperatura impuesta en el extremo izquierdo
Tderecha  =30; % Temperatura impuesta en el extremo derecho

Longitud=1; % en metros
Tiempo=1000; % total de la integracion, en segundos
Difusividad=1e-4; % parametro de difusividad de temperatura para cobre

Nx=50; % Numero de puntos de grilla a lo largo de la barra
CFL=0.5; % Es el valor de Difusividad*dt/dx/dx, el cual debe ser positivo 
         & y menor o igual que 0.5 para que el metodo sea estable
nskip=20; % Cada cuantos pasos de tiempo actualizar el grafico

%%%%%%%%%%%%%%%%%

% Inicializar algunos valores
dx=Longitud/(Nx-1); x=(0:dx:Longitud)';
dt=CFL*dx*dx/Difusividad;
T=Tinicial*ones(Nx,1); T(1)=Tizquierda; T(end)=Tderecha;
A=spdiags([CFL*ones(Nx-2,1) (1-2*CFL)*ones(Nx-2,1) CFL*ones(Nx-2,1)],...
         [-1 0 1],Nx-2,Nx-2);
B=zeros(Nx-2,1); 
B(1)=CFL*Tizquierda; 
B(end)=CFL*Tderecha;



%%%%%%%%%%%%%%%%%

% Integrar en el tiempo

for t=0:dt:Tiempo

   % Graficar
   if (mod(t/dt,nskip)==0)
   bar(x,T)
   xlim([0 Longitud]),ylim([0 max([Tizquierda Tinicial Tderecha])+5])
   xlabel('Distancia desde extremo izquierdo (m)'),ylabel('Temperatura (grados)')
   title(['Tiempo = ' num2str(round(t)) ' segundos'])
   drawnow
   num= num2str(t/dt/nskip);
   if (t/dt/nskip<10); num=['0' num]; end
   print(['modelo.' num '.png'],'-dpng','-S320,240')
   end

   % Calcular valor siguiente
   T(2:end-1)=A*T(2:end-1)+B;

end