% ---------------------------------------------------------------
%     Forward-backward approximation to the rotating SW eqns
%     in a periodic domain using a staggered mesh
% ---------------------------------------------------------------

clear all;

% user specified parameters

nx = 1000;          % number of grid points (not including the repeated point 
                  % at the begining of the periodic domain).  
min_x = 0;        % coordinate of left edge of computational domain
max_x = 5000e3;        % coordinate of right edge of computational domain

cdx=0.9;            % Courant number

t_stop = 110000;       % end time for 2-waves

% windowing parameters for the plot

min_u = -1.2;   % minimum value for plot of dependent variable
max_u =  1.2;   % maximum value for plot of dependent variable

min_v = -1.5;   % minimum value for plot of dependent variable
max_v =  1.5;   % maximum value for plot of dependent variable

min_h = -1.5;   % minimum value for plot of dependent variable
max_h =  1.5;   % maximum value for plot of dependent variable

% min_x_wndo = 0;   % coordinate of left edge of plotted subdomain
% max_x_wndo = 5000e3;   % coordinate of right edge of plotted subdomain

min_x_wndo = 1200e3;   % coordinate of left edge of plotted subdomain
max_x_wndo = 3800e3;   % coordinate of right edge of plotted subdomain


% ---------------------------------------------------------------
%                     Initialization
% ---------------------------------------------------------------

% Physical constants
fcor=1e-4;  % Coriolis force
g=9.8;      % gravitational acceleration
cp=10;      % SW phase speed
H=cp^2/g;   % Mean fluid depth


% generate the discrete x axis

dx=(max_x-min_x)/nx;
xpt = linspace(min_x+dx,max_x,nx);
xpt_h = xpt-dx/2;   % staggered x points for the height field

dt = cdx*dx/cp;     % recover time step from the Courant number
n_steps = t_stop/dt;  % number of time steps

% generate the initial field 

for ix = 1:nx
        initial_h(ix,1) = square(2*pi*xpt_h(ix)/max_x);
end
 
% smooth the height field

for ism = 1:4
    initial_h=smooth(initial_h);
end
    
 initial_u = zeros(nx,1);
 initial_v = zeros(nx,1);

% plot initial condition

figure(1);
plot(xpt,initial_h,'k:');
title('Initial free surface');
axis([min_x_wndo max_x_wndo min_h max_h]);


% -------------------------------------------------------------
% solve forward-backward system
% -------------------------------------------------------------

u_old=initial_u;
v_old=initial_v;
h_old=initial_h;
time =0;

for it=1:n_steps
    
    time=time+dt;
% -------------------------------------------------------------
% solve forward-backward system
% -------------------------------------------------------------
 
    [u_new,v_new,h_new] =forward_backward(u_old,v_old,h_old,g,fcor,...
        H,dt,dx);

% -------------------------------------------------------------
% display animated plots 
% -------------------------------------------------------------
    figure(2);
    plot(xpt_h,h_new,'r');
    nstr1 = num2str(time,3); title(['Free surface height h at time = ',nstr1]);
    axis([min_x_wndo max_x_wndo min_h max_h]);
    %pause(t_anim); 
    
    figure(3);
    plot(xpt,u_new);
    nstr1 = num2str(time,3); title(['u at time = ',nstr1]);
    axis([min_x_wndo max_x_wndo min_u max_u]);
    %pause(t_anim); 
    
    figure(4);
    plot(xpt,v_new);
    nstr1 = num2str(time,3); title(['v at time = ',nstr1]);
    axis([min_x_wndo max_x_wndo min_v max_v]);
    %pause(t_anim); 
    
% -------------------------------------------------------------
% update arrays for next step
% -------------------------------------------------------------
    
    u_old=u_new; v_old=v_new; h_old=h_new;
end


