%function [] = prob_13_4()
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; 

% Select a dummy point for example purposes:
%
qm = [ 0.5; 0.5 ]; 
alphas = linspace(0,1,200); 

g=1; % Gravity

% q^l lies on the two-shock centered on q^m:
ql = zeros(2,length(alphas));
ql(1,:) = qm(1) + alphas; 
ql(2,:) = qm(2) + alphas .* ( qm(2)/qm(1) + sqrt( g*qm(1)*(1+alphas/qm(1)).*(1+alphas/(2*qm(1))) ) );

% q^r lies on the one-shock centered on q^m:
qr = zeros(2,length(alphas));
qr(1,:) = qm(1) + alphas; 
qr(2,:) = qm(2) + alphas .* ( qm(2)/qm(1) - sqrt( g*qm(1)*(1+alphas/qm(1)).*(1+alphas/(2*qm(1))) ) );

% INDICES where we will take q^r and q^l to be defind from on thier respective Hugoniot curves
qri=60; 
qli=80;

% Plot the state q^m (in green):
figure; plot( qm(1), qm(2), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g' ); hold on;
% Plot the Hugoniot loci for the 
% -- two shock centered at q^m and the 
% -- one shock centered at q^m:
h1=plot( ql(1,:), ql(2,:), '-ob' ); 
h2=plot( qr(1,:), qr(2,:), '-or' ); 

% Plot the location of the state q^l:
plot( ql(1,qli), ql(2,qli), 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b' ); 
% Plot the location of the state q^r:
plot( qr(1,qri), qr(2,qri), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r' ); 

% Plot the two-shock through the point q^r:
qr2(1,:) = qr(1,qri) + alphas; 
qr2(2,:) = qr(2,qri) + alphas .* ( qr(2,qri)/qr(1,qri) + sqrt( g*qr(1,qri)*(1+alphas/qr(1,qri)).*(1+alphas/(2*qr(1,qri))) ) );
h3=plot( qr2(1,:), qr2(2,:), '-.xb' ); 

% Plot the one-shock through the point q^l:
ql1(1,:) = ql(1,qli) + alphas; 
ql1(2,:) = ql(2,qli) + alphas .* ( ql(2,qli)/ql(1,qli) - sqrt( g*ql(1,qli)*(1+alphas/ql(1,qli)).*(1+alphas/(2*ql(1,qli))) ) );
h4=plot( ql1(1,:), ql1(2,:), '-.xr' ); 

grid on; 
axis( [ 0.45 1.5 -0.5 3.5 ] ); 
xlabel( 'h' ); ylabel( 'h u' ); 

% Add a legend:
legend( [h1,h2,h3,h4], ...
	{'two-shock centered on q^m',...
	 'one-shock centered on q^m',...
	 'two-shock centered on q^r',...
	 'one-shock centered on q^l',}, 'Location', 'NorthWest' );

% Add some text:
text( 0.5, 0.25, 'q^m', 'FontSize', 15 )
text( 0.9, 1.75, 'q^l', 'FontSize', 15 )
text( 0.8, 0.1 , 'q^r', 'FontSize', 15 )
text( 1.3, 1.1 , 'q^{m''}', 'FontSize', 15 )

saveas( gcf, 'chap_13_prob_4.eps', 'psc2' );