%% EJERCICIO 5.9: Dos cuerdas semiinfinitas de distinta densidad

%% EJERCICIO 5.9 grafico los coeficientes de reflexion y transmision
clear all

%defino variables
w=1;                     %frecuencia angular
Tension=1;                    %tension
mu1=1;                  %densidad de la cuerda 1
mu2=.001:.01:50;                  %densidad de la cuerda 2
v1=sqrt(Tension./mu1);                   %velocidad en medio 1
v2=sqrt(Tension./mu2);                 %velocidad en medio 2
k1=w./v1;                %numero de onda en medio 1
k2=w./v2;                %numero de onda en medio 2
I=1;                    %amplitud incidente
T=2*k1*I./(k1+k2);      %amplitud transmitida
R=I*(k1-k2)./(k1+k2);   %amplitud reflejada

%grafico los coeficientes de reflexion y transmision
figure(1);clf
plot(k2,T,k2,R)
legend('T','R')
line(xlim,[0 0],'linestyle',':','color','k')
line(xlim,[1 1],'linestyle',':','color','k')
line([1 1],ylim,'linestyle',':','color','k')
xlabel('k2')
ylabel('Coeficientes')
text(6,-.6,'Reflejado')
text(6,.4,'Transmitido')

%% EJERCICIO 5.9, grafico en funcion del tiempo
clear all

%defino variables
w=1;                    %frecuencia angular
Tension=1;              %tension
mu1=1;                  %densidad de la cuerda 1
mu2=10; %.001 0.1 1 10 100 1000  %densidad de la cuerda 2
v1=sqrt(Tension./mu1);        %velocidad en medio 1
v2=sqrt(Tension./mu2);        %velocidad en medio 2
k1=w/v1;                %numero de onda en medio 1
k2=w/v2;                %numero de onda en medio 2
I=1;                    %amplitud incidente
T=2*k1*I./(k1+k2);      %amplitud transmitida
R=I*(k1-k2)./(k1+k2);   %amplitud reflejada

%defino las funciones Incidente, Reflejada, Transmitida y Total
fiI=@(x,t) I*cos(k1.*x-w.*t) .* (x<0);
fiR=@(x,t) R*cos(k1.*x+w.*t) .* (x<0);
fiT=@(x,t) T*cos(k2.*x-w.*t) .* (x>=0);
fi=@(x,t) fiI(x,t)+fiR(x,t)+fiT(x,t);

X=-10:.01:10;
figure(1);clf
%grafico para cada tiempo
for t=0:.1:20
    subplot(4,1,1)
    plot(X,fiI(X,t))
    ylim([-2 2])
    title(sprintf('Incidente I=%0.2f',I))
    subplot(4,1,2)
    plot(X,fiR(X,t))
    ylim([-2 2])
    title(sprintf('Reflejada R=%0.2f',R))
    subplot(4,1,3)
    plot(X,fiT(X,t))
    ylim([-2 2])
    title(sprintf('Transmitida T=%0.2f',T))
    subplot(4,1,4)
    plot(X,fi(X,t))
    ylim([-2 2])
    title(sprintf('TOTAL t=%0.2f',t))
    text(min(xlim)+.1*range(xlim),1,sprintf('k1=%0.2f\nv1=%0.2f',k1,v1))
    text(max(xlim)-.2*range(xlim),1,sprintf('k2=%0.2f\nv2=%0.2f',k2,v2))
    line([0 0],xlim,'linestyle',':','color','k')
    drawnow
%     pause(0.01)
end

%% EJERCICIO 5.9 Grafico energia
clear all

%defino variables
w=1;                    %frecuencia angular
Tension=1;              %tension
mu1=1;                  %densidad de la cuerda 1
mu2=10; %.001 0.1 1 10 100 1000  %densidad de la cuerda 2
v1=sqrt(Tension./mu1);  %velocidad en medio 1
v2=sqrt(Tension./mu2);  %velocidad en medio 2
k1=w/v1;                %numero de onda en medio 1
k2=w/v2;                %numero de onda en medio 2
I=1;                    %amplitud incidente
T=2*k1*I./(k1+k2);      %amplitud transmitida
R=I*(k1-k2)./(k1+k2);   %amplitud reflejada

%defino las funciones Incidente, Reflejada, Transmitida y Total
fiI=@(x,t) I*cos(k1.*x-w.*t) .* (x<0);
fiR=@(x,t) R*cos(k1.*x+w.*t) .* (x<0);
fiT=@(x,t) T*cos(k2.*x-w.*t) .* (x>=0);
fi= @(x,t) fiI(x,t)+fiR(x,t)+fiT(x,t);

%defino la funcion energia potencial
Ep=@(x,t) (1/2*Tension*k1^2*(I*sin(k1.*x-w.*t)+R*sin(k1.*x+w.*t)).^2) .* (x<0) + ...
          (1/2*Tension*k2^2*(T*sin(k2.*x-w.*t)).^2) .* (x>=0); 

%defino la funcion energia cinetica
Ec=@(x,t) (1/2*mu1*w^2*(-I*sin(k1.*x-w.*t)+R*sin(k1.*x+w.*t)).^2) .* (x<0) + ...
          (1/2*mu2*w^2*(-T*sin(k2.*x-w.*t)).^2) .* (x>=0); 

%calculo la energia maxima, para la escala
Emax=1.1*max(abs(Tension*k1^2*(I+R)^2),abs(Tension*k2^2*T^2));


%
% set(gcf, 'Renderer','painters') 
% set(gcf, 'Renderer','opengl') 

figure(2);clf
% h(1)=axes('position',[.1 .6 .85 .35]);
% h(2)=axes('position',[.1 .1 .85 .35]);

% set(f, 'GraphicsSmoothing','off') 
X=-10:.01:10;
myxlim=X([1 end]);
%grafico para cada tiempo
tic
for t=0:.05:5
%     set(gcf,'currentaxes',h(1))
    subplot(2,1,1)
    plot(X,fi(X,t))
    ylim([-2 2])
    line([0 0],myxlim,'linestyle',':','color','k')    
    line(myxlim,[0 0],'linestyle',':','color','k')  
    xlabel('desplazamiento')
    ylabel('Energía')
    
%     set(gcf,'currentaxes',h(2))
    subplot(2,1,2)
    myEp=Ep(X,t);
    myEc=Ec(X,t);
    plot(X,myEp,X,myEc,X,myEp+myEc)
    ylim([0 Emax])
    xlabel('Tiempo')
    ylabel('Energía')
    line([0 0],myxlim,'linestyle',':','color','k')
    title(sprintf('TOTAL t=%0.2f',t))
    text(min(myxlim)+.1*range(myxlim),1,sprintf('k1=%0.2f\nv1=%0.2f',k1,v1))
    text(max(myxlim)-.2*range(myxlim),1,sprintf('k2=%0.2f\nv2=%0.2f',k2,v2))

%     legend('Potencial','Cinetica','Total')
    drawnow
%     pause(0.01)
end
toc