function iso_euler(rho_l,m_l,rho_r,m_r,c)
%
% Matlab-File zur Illustration der 1. Aufgabe des
% Projektblattes "Isotherme Euler-Gleichung" zur Vorlesung
% "Theorie und Numerik hyperbolischer Erhaltungsgleichungen"
% M.Ehrhardt, WS 2009/10
%
% Phasenraum der isothermen Eulergleichungen
% Plot der 1-,2- Schocks/Verduennungswellen
%
% linker Zustand (rho_l,m_l)
% rechter Zustand (rho_r,m_r)
% Schallgeschwindigkeit c
%
% Beispielaufruf: iso_euler(1,1,3,2,1)
%

holdflag = ishold; % 1 falls "hold on", 0 falls "hold off"
% Plotte linken Zustand u_l (gruen)
plot(rho_l,m_l,'go');hold on

% Plotte rechten Zustand u_r (magenta)
plot(rho_r,m_r,'mo');

% Berechne 1-Welle
[rho_m,rho_p,m_m,m_p] = kurven(rho_l,m_l,c,1);
S_1 = plot(rho_p,m_p,'b-');
V_1 = plot(rho_m,m_m,'r--');

% Berechne 2-Welle
[rho_m,rho_p,m_m,m_p] = kurven(rho_r,m_r,c,2);
S_2 = plot(rho_p,m_p,'b:');
V_2 = plot(rho_m,m_m,'r-.');

cx = (max(min(rho_l,rho_r)-2,0) + max(rho_l,rho_r)+2)/2;
cy = (max(m_l,m_r)+min(m_l,m_r))/2;
w  = max(max(rho_l,rho_r)+2-max(min(rho_l,rho_r)-2,0),max(m_l,m_r)-min(m_l,m_r)+4);
cx = max(cx,w/2);
axis([cx-w/2,cx+w/2,cy-w/2,cy+w/2]);
axis square;

legend([S_1 V_1 S_2 V_2],'1-Schock','1-Verd.welle','2-Schock','2-Verd.welle')
text(rho_l+w/50,m_l-w/50,'u_l')
text(rho_r+w/50,m_r-w/50,'u_r')
title('isotherme Eulergleichungen: Schocks/Verd.wellen im Phasenraum')
xlabel('Dichte \rho')
ylabel('Impuls(-dichte) m')

% "hold" wieder in alten Zusand setzen
if holdflag == 0
  hold off
end


% Unterroutinen
% --------------------------------------------------------------------------

function [rho_m,rho_p,m_m,m_p] = kurven(rho_l,m_l,c,p)
% berechnet 1-,2- Schocks/Verduennungswellen
% der isothermen Eulergleichungen
%
% Numerierung bezieht sich auf Buch von Leveque
%
if p==1
  cc = c;
elseif p==2
  cc = -c;
else
  error('p: invalid value!')    
end

xi    = 0:0.1:10;
rho_p = rho_l * (ones(size(xi))+xi); % (7.11)
% (7.12) bzw. (7.13)
m_p   =   m_l * (ones(size(xi))+xi) - cc*rho_l*xi.*sqrt(ones(size(xi))+xi);

if p==1
  rho_m = rho_l * exp(-xi/cc);                                % (8.13)
  m_m   = rho_l * (xi+m_l/rho_l*ones(size(xi))).*exp(-xi/cc); % (8.15)
elseif p==2
  % Analoge Formeln (1),(2) aus der 3.Aufgabe der 2.Uebung:
  rho_m = rho_l * exp(xi/cc);
  m_m   = rho_l * (-xi+m_l/rho_l*ones(size(xi))).*exp(xi/cc);
else
  error('p: invalid value!')  
end













