function integrando_ec_calor
%Integrando la ecuación del calor, siguiendo la guia de difusividad
%térmica:  http://materias.df.uba.ar/labo4Aa2013c2/files/2013/08/Difusividad.pdf
%Modificado del ejemplo de la funcion pdepe, ejecutar: doc('pdepe').
%
%1- Comparamos la solución analítica que figura en la guía, con la solución
%numérica, obtenida integrando la ecuación diferencial en derivadas
%parciales.
%2- Agrego pérdidas
%3- Agrego fuente switcheada
%4- Agrego fuente switcheada y pérdidas
%
%Diego Shalom, Junio/2014


%defino tiempos y posiciones
x = 0:.5:10;
t = 0:.2:60;

%solucion analitica, que aparece en la guia
figure(1);clf
disp('Sol1: Solucion Analitica, que aparece en la guia')
sol1=solucion_analitica(x,t);
set(gcf,'name','Solucion analitica')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol1)

%solucion numerica, integrando la ecuacion diferencial
figure(2);clf
disp('Sol2: Solucion numerica,integrando la ecuacion del calor')
sol2 = pdepe(0,@fcn_ecdif,@fcn_condini,@fcn_condcontorno,x,t);%integro esta ecuación diferencial, con estas condiciones cond de contorno e iniciales, en esta grilla de x, t
set(gcf,'name','Solucion numerica')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol2)

% la diferencia entre la solucion analitica y la numerica
figure(3);clf
disp('Sol2-Sol1: Diferencia entre la solucion analitica y la numerica')
set(gcf,'name','Solucion numerica')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol2-sol1)

%solucion numerica, con perdidas!!!!
figure(4);clf
disp('Sol3: Solucion numerica, con perdidas')
sol3 = pdepe(0,@fcn_ecdif_conperdidas,@fcn_condini,@fcn_condcontorno,x,t);%integro esta ecuación diferencial, con estas condiciones cond de contorno e iniciales, en esta grilla de x, t
set(gcf,'name','Solucion numerica, con perdidas!!')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol3)

%solucion numerica, con fuente on/off!!!!
figure(5);clf
disp('Sol4: Solucion numerica, con fuente switcheada')
sol4 = pdepe(0,@fcn_ecdif,@fcn_condini,@fcn_condcontorno_onoff,x,t);%integro esta ecuación diferencial, con estas condiciones cond de contorno e iniciales, en esta grilla de x, t
set(gcf,'name','Solucion numerica, con fuente on/off!!')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol4)

%solucion numerica, con perdidas y fuente on/off!!!!
figure(6);clf
disp('Sol5: Solucion numerica, con perdidas y fuente switcheada')
sol5 = pdepe(0,@fcn_ecdif_conperdidas,@fcn_condini,@fcn_condcontorno_onoff,x,t);%integro esta ecuación diferencial, con estas condiciones cond de contorno e iniciales, en esta grilla de x, t
set(gcf,'name','Solucion numerica, con perdidas y fuente on/off!!')
set(gcf,'position',[100,300,1000,500])
grafico_una_solucion(x,t,sol5)

disp('listo!')
end
% --------------------------------------------------------------
function sol=solucion_analitica(x,t)
xmesh=repmat(x,length(t),1);%para la solucion analítica necesito grids
tmesh=repmat(t',1,length(x));%para la solucion analítica necesito grids
kappa=1/pi^2;%defino este kappa para que coincida con la ecuación que integro
fcn_solucion_analitica=@(x,t,kappa) 2*sqrt(kappa*t/pi).*exp(-x.^2./(4*kappa*t)) -  ...
                2*x./2.*erfc(x./(2*sqrt(kappa*t)));%la funcion de la solucion analitica
sol=fcn_solucion_analitica(xmesh,tmesh,kappa);
end
% --------------------------------------------------------------
function [c,f,s] = fcn_ecdif(x,t,u,DuDx) %ecuacion diferencial
%corresponde a la ecuacion del calor:
%   pi^2*du/dt=d2u/dx2
c = pi^2;
f = DuDx;
s = 0;%sin perdidas
end
% --------------------------------------------------------------
function [c,f,s] = fcn_ecdif_conperdidas(x,t,u,DuDx) %ecuacion diferencial con perdidas
%corresponde a la ecuacion del calor con perdidas:
%   pi^2*du/dt=d2u/dx2-u
c = pi^2;
f = DuDx;
s = -u;%con perdidas (asumo que la temperatura ambiente es cero)
end
% --------------------------------------------------------------
function u0 = fcn_condini(x) %condicion inicial
%condicion inicial, temperatura inicial toda cero (asumo temperatura ambiente=0)
u0=0;
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = fcn_condcontorno(xl,ul,xr,ur,t) %condiciones de contorno
%condiciones de contorno con flujo de calor constante entrando por la
%izquierda y tempertura=0 en el borde derecho
pl = 1;% en el borde izquierdo, du/dx=1 (asumo que F0/K=1), un flujo de calor entrante, constante 
ql = 1;
pr = ur;
qr = 0;%en el borde derecho, la funcion u es cero, la temperatura es cero
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = fcn_condcontorno_onoff(xl,ul,xr,ur,t) %condiciones de contorno con fuente on/off
%condiciones de contorno de la izquierda una fuente que se prende y se
%apaga con un dado periodo y dutycycle. y tempertura=0 en el borde derecho

periodo=10;%la duracion del periodo del on/off
dutycycle=0.5;%la fraccion del periodo que la fuente esta prendida
fase=mod(t,periodo)/periodo;%la fase dentro de cada periodo
pl=fase<dutycycle;%du/dx=1 o du/dx=0 , dependiendo de la fase y el dutycycle

ql = 1;
pr = ur;
qr = 0;%en el borde derecho, la funcion u es cero, la temperatura es cero
end
% --------------------------------------------------------------
function grafico_una_solucion(x,t,sol)
subplot(1,3,1)
surf(x,t,sol) 
title('Temperature vs time and position ')
xlabel('Position x')
ylabel('Time t')
zlabel('Temperature u(x,t)')

subplot(1,3,2)
plot(x,sol)
title('Solution at various times')
xlabel('Distance x')
ylabel('u(x,t)')

subplot(1,3,3)
plot(t,sol')
title('Solution at various positions')
xlabel('Time t')
ylabel('u(x,t)')
end