%% 1D LSWE dam break simulation using Fourier decomposition
%  Assumes f = H = 1 (x in units of R = sqrt(gH)/f, t in units of 1/f).

%% Initialize dam break simulation

global ETARH ETALH ETAGH VRH VLH VGH X K XMIN XMAX YMAX CLIM BETA U

% User-specified parameters

Lx = 8*2*pi; % Length of periodic domain
nx = 128;    % Number of gridpoints

eta0 = 0.1;  % Half-amplitude of initial step in eta
x1 = Lx/4;   % Location of initial step
CLIM = 0.1;

% Other parameters

BETA = 0;  % f-plane
U = 0;  % no mean flow
dx = Lx/nx;
X = (0:nx)*dx; % Extended gridpoints (includes x=Lx as well as 0) for plot 
xp = X(1:nx);  % Fourier Gridpoints of periodic domain of length Lx
XMIN = X(1);
XMAX = X(nx/2 + 1);
YMAX = 1.5;
K = (2*pi/Lx)*[0:(nx/2-1) (-nx/2):(-1)];

eta = eta0*sign(x1-xp).*sign(x1+Lx/2-xp); % Even extn of eta(x,0) to [0 Lx]

% Fourier components 

ETAH = fft(eta);
zero = zeros(1,nx);
UH = zero;
VH = zero;

cp = sqrt(1 + K.^-2);
ag = (ETAH - 1i.*K.*VH)./(K.^2 + 1);
ar = 0.5*(ETAH + UH./cp) - 0.5*ag;
al = 0.5*(ETAH - UH./cp) - 0.5*ag;

ETARH = ar;
ETALH = al;
ETAGH = ag;

VRH = -1i*ar./K;
VLH = -1i*al./K;
VGH = 1i*K.*ag;

% Initialize m=1 (constant) Fourier component to overwrite divide-by-zeros

ETARH(1) = 0;
ETALH(1) = 0;
ETAGH(1) = 0;

VRH(1) = 0;
VLH(1) = 0;
VGH(1) = 0;

%% Plot initial condition to explain color scheme

clf
waveplot(0)

%% Make movie

Dmovie = makemovie('waveplot',0:0.02:10);
