%
%  w10shock.m -- (djm - 22 mar 10)
%

clear

global alpha gam rr

cc = 1.0;  gam = 7/5;  alpha = cc/3;  rr = cc^(2/(gam-1));
Box = 6;

do_more = 3;

figure(1);  clf
hold on

%  piston motion
tp = 0.0:0.05:Box;
plot(alpha/2*tp.^2,tp,'b','linewidth',3)
axis([0 5 0 6])

% undisturbed characteristics, plot from axis
xa = 0.2:0.2:2*Box;
for xxa = xa
    plot([xxa/2 xxa],cc*[xxa/2 0],'g','linewidth',1)
    plot(xxa + [0 Box],cc*[0 Box],'k','linewidth',1)
end

switch do_more
    case{1,2,3}
        %  compression zone R+ characteristics from piston 
        tp_list = 0.2:0.2:Box;
        
        for tp = tp_list;
            %  piston endpoint
            xp = alpha/2*tp^2;
            up = alpha*tp;
            cp = (gam-1)/2 * up + cc;
            
            %  far point
            tp2 = Box;
            xp2 = xp + (up + cp)*(tp2-tp);
            
            plot([xp xp2],[tp tp2],'k','linewidth',1)
        end
end

switch do_more
    case{2,3}
        %  envelope boundary of crossing R+ characteristics
        tp_l = 0.0:0.01:Box/2;
        
        xp_l = (alpha/2)*tp_l.^2;
        
        %  magic envelope formula
        te = (2/(gam+1))*(gam*tp_l + cc/alpha);

        up_l = alpha*tp_l;
        cp_l = (gam-1)/2 * up_l + cc;
        
        xe = xp_l + (up_l + cp_l).*(te-tp_l);
                    
        plot(xe,te,'m--','linewidth',2)
end

%  first characteristic from piston
plot([0 Box],cc*[0 Box],'r','linewidth',2)

%  shock point (t = 2.5)
plot(xe(1),te(1),'k.','markersize',10)

title('\bf R^+ characteristics')
xlabel('\bf x-axis')
ylabel('\bf t-axis')

switch do_more
    case{3}
        figure(2);  clf;  axis([0 5 -0.5 5.5])

        %  solution profiles
        for tt = 0.5:0.5:4
            t0 = 0:tt/80:tt;
            
            %  piston location at tt
            xx = (alpha/2)*t0.^2;
            
            %  parametric solution
            u0 = alpha*t0;
            c0 = (gam-1)/2 *u0 + cc;
            x0 = xx + (u0 + c0).*(tt-t0);
            
            %  pre-pend constant state
            x0 = [5 x0];  u0 = [0 u0];
            
            figure(2);  hold on
            plot([0 5],tt*[1 1],'k:');  hold on
            plot(x0,tt+u0,'k','linewidth',2)
            plot(x0(end)*[1 1],tt+[0 u0(end)],'k:')
            plot(x0(end),tt,'b.','markersize',10)
            plot(x0(2),tt,'r.','markersize',10)

            %  envelope point
            if tt > 2/(gam+1)/alpha
                te = (tt*(gam+1)/2 - cc/alpha)/gam;
                ue = alpha*te;
                xe = (alpha/2)*te^2 + ((gam+1)/2*ue + cc).*(tt-te);            
    
                plot(xe*[1 1],tt+[0 ue],'k:')
                plot(xe,tt,'m.','markersize',10)
            end
        end
        
        figure(2);
        plot([0 Box],[0 0],'k','linewidth',2)
        title('\bf a problem with the solution ...')
        ylabel('\bf shifted u')
        xlabel('\bf x-axis')
end