%
%  w12shock.m -- djm -- 21 nov 2010 
%

%  grid parameters
L = 2*pi;  N = 36;  dx = L/N;

%  coordinates
xx = 0:dx:L;

%  initial condition
x0 = xx;  u0 = sin(x0);

figure(1);  clf
subplot(2,1,1)
plot(x0,u0,'b',x0,u0,'k.');  hold on
axis([0 L -1.1 1.1])
title(['\bf nonlinear wave'],'fontsize',12)
xlabel(['\bf x-axis'],'fontsize',10)
ylabel(['\bf u(x,t)'],'fontsize',10)

%  characteristics plot
te = 3/2;
xe = x0 + u0*te;
ue = u0;

subplot(2,1,2)
for jj = 1:size(xx,2)
	plot([x0(jj) xe(jj)],[0 te],'k');  hold on
	plot([x0 ],[0],'k.');
end
plot([0 L],[0 0],'b')

%  next times
for t1 = [1/3 2/3]; 
	x1 = x0 + u0*t1;
	u1 = u0;
	subplot(2,1,1)
	plot(x1,u1,'g',x1,u1,'kx');
	subplot(2,1,2)
	plot([0 L],[t1 t1],'g')
	plot(x1,t1,'kx')
end

for t1 = [1]; 
	x1 = x0 + u0*t1;
	u1 = u0;
	subplot(2,1,1)
	plot(x1,u1,'r',x1,u1,'kx');
	subplot(2,1,2)
	plot([0 L],[t1 t1],'r')
	plot(x1,t1,'kx')
end

for t1 = [4/3]; 
	x1 = x0 + u0*t1;
	u1 = u0;
	subplot(2,1,1)
	plot(x1,u1,'m',x1,u1,'kx');
	subplot(2,1,2)
	plot([0 L],[t1 t1],'m')
	plot(x1,t1,'kx')
end

te = 3/2;
xe = x0 + u0*te;
ue = u0;

subplot(2,1,2)
for jj = 1:size(xx,2)
	plot([x0(jj) xe(jj)],[0 te],'k');  hold on
	plot([x0 ],[0],'k.');
end
axis([0 L 0 3/2])
title(['\bf characteristic curves'],'fontsize',12)
xlabel(['\bf x-axis'],'fontsize',10)
ylabel(['\bf t-axis'],'fontsize',10)