function Trafico

L = .1; % Tamanio del dominio
N = 500; %Numero de sub-intervalos en la discretizacion
T = .1;
x_ini = -L; x_fin = L;
dx = (x_fin-x_ini)/N;
u_max = 1;

%% Definir el eje x
x = zeros(N,1);
for j = 1:N
    x_j = x_ini+dx*( j-1/2 );
    x(j,1) = x_j;
end

%% Definir la condicion inicial

rho0 = zeros(N,1);
rho_l = 0.4; %0.6;
rho_r = 0.2; %0.8;
for j=1:N
    x_j=x(j,1);
    if x_j < 0
        rho0(j,1) = rho_l;
    else
        rho0(j,1) = rho_r;
    end
end
lambda_max = max(abs(u_max*(1-2*rho0))); %Es la nueva a
nu = 0.9; %Numero de Courant
dt = nu*dx/lambda_max;

%% Solucion exacta

rho_exacta = zeros(N,1);
if rho_l < rho_r
    s = 1-(rho_l+rho_r); %Velocidad de la onda de choque
    for j=1:N
        x_j=x(j,1);

        if x_j < s*T
            rho_exacta(j,1) = rho_l;
        else
            rho_exacta(j,1) = rho_r;
        end
    end
else %Aqui tendriamos una rarefaccion
    for j=1:N
        x_j=x(j,1);
        xi_l = u_max*(1-2*rho_l); %xi_l = f'(rho_l)
        xi_r = u_max*(1-2*rho_r); %xi_r = f'(rho_r)
        if x_j/T <= xi_l %Esta a la izquierda del abanico
            rho_exacta(j,1) = rho_l;
        end
        if xi_l <= x_j/T && x_j/T < xi_r %Esta dentro del abanico
            rho_exacta(j,1) = 1/2-1/(2*u_max)*x_j/T; %f'(a(xi)) = xi 
        end
        if  xi_r <= x_j/T %Esta a la derecha del abanico
            rho_exacta(j,1) = rho_r;
        end
    end
end

%% Metodo upwind

t = dt;
rho_n = rho0;
rho_np1 = rho_n;
while t < T
    for j=2:N-1
        lambda_bar = u_max*( 1- (rho_n(j-1,1)+rho_n(j,1)));
        if lambda_bar > 0
            rho_np1(j,1) = rho_np1(j,1) - lambda_bar*dt/dx*(rho_n(j,1)-rho_n(j-1,1));
        else
            rho_np1(j-1,1) = rho_np1(j-1,1) - lambda_bar*dt/dx*(rho_n(j,1)-rho_n(j-1,1));
        end
    end
    %Condiciones de frontera Nerhomann
    rho_np1(1,1) = rho_np1(2,1);
    rho_np1(N,1) = rho_np1(N-1,1);
    
    rho_n = rho_np1;
    lambda_max = max(abs(u_max*(1-2*rho_n)));
    dt = nu*dx/lambda_max;
    t = t+dt;
end
   

%% Graficar condicion inicial y solucion aproximada a tiempo T

figure(1)
hold on
plot(x,rho_exacta,'b -','LineWidth',2)
plot(x,rho_n,'r *','LineWidth',2)
xlabel('x','FontSize',20)
ylabel('Solucion','FontSize',20)
title('Condicion inicial y solucion','FontSize',20)
lg = legend('Condicion inicial','Solucion a tiempo T');
set(lg,'Location','NorthWest');
set(lg,'FontSize',15);
%axis([-L L 0.6 0.8])
%print('OndaChoque','-depsc',figure(1))
print('Rarefaccion','-depsc',figure(1))
