% Plot dominant balance solution of LSWE for sinusoidal tsunami propagating
% onto a continental shelf H(x) = H0*tanh^2(x/W), -infty < x < 0.
% Solution has form
%  eta(x) = Re(h0*sqrt(c0/c(x))*cos(phi)),
%  phi(x) = om*(int_0^x dx/c(x) - t), with time t = 0.
% and should asymptote to h = h0*cos(om*x/c0) for x -> -infty.
   % Parameters in MKS units
   
   h0 = 1;
   g = 9.8;
   H0 = 4000;
   c0 = sqrt(g*H0);
   W = 100*1e3;
   om = 0.025;
   
   % Computations start here
 
   xilarge = 10;          % Value of xi = -x/W at which phase
                          % matched to desired x->-infty behavior.
   x = W*(-xilarge:0.0002:0);  % Go to -xilarge*W for accurate phase match
                          % with desired x->-infty behavior
			  % and use fine dx for near-coast plotting.
   nx = length(x);
   H = H0*tanh(-x/W);
   c = sqrt(g*H);
   hamp = h0*sqrt(c0./c);
   dx = x(2) - x(1);
   k  = om./c;
   km = [0 k(1:nx-1)];
   
   % Phase computation such that phi = (om/c0)*x at x = -xilarge*W.
   phi = 0.5*cumsum(k+km)*dx + x(1)*om/c0;
   eta= hamp.*cos(phi);

   % Aproximations for nonlinearity threshold xn (where hamp(xn) = H(xn))
   % and asymptotic accuracy threshold xom (where |dc/dx(xom)| = om)
   
   xn = -W*(h0/H0)^0.8
   xom = g*H0/(4*om^2*W)

   % Make plots over two different ranges of x, [-W 0} and [-0.05W 0]

   subplot(2,2,1) 
   plot(x(x<xn)/1000,eta(x<xn),'b-',...
	x(x<xn)/1000,hamp(x<xn),'b--',...
	x(x<xn)/1000,-hamp(x<xn),'b--',...
	x/1000,-H,'k-')
   xlim(W*[-1 0]/1000)
   ylim(h0*[-5 5])
   xlabel('x [km]')
   ylabel('\eta(x,0) [m]')
   legend('\eta(x,0)','\pm |h(x)|','Location','SouthWest')
   
   subplot(2,2,2)
   plot(x(x<xn)/1000,eta(x<xn),'b-',...
	x(x<xn)/1000,hamp(x<xn),'b--',...
	x/1000,-H,'k-',...
	x(x<xn)/1000,-hamp(x<xn),'b--',...
	[xn xn]/1000,h0*[-10 10],'k:')
   pt2km = 0.2;
   text(xn/1000-pt2km, 8*h0,'x_n')
   xlim(W*[-0.05 0]/1000)
   ylim(h0*[-10 10])
   xlabel('x [km]')
   ylabel('\eta(x,0) [m]')
   legend('\eta(x,0)','\pm |h(x)|','-H(x)','Location','SouthWest')