% test CELT.m and/or CELT_Q.m on a linear slope
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some model parameters
omega=2*pi/(12.4*3600); % Wave frequency rad/sec
f=0; % Inertial frequency
nu=0;%10^(-5); % Viscosity (can be zero)
rho0=1000; % Reference density
Nm=16; % Number of modes in solution 
Nx=32; % Number of steps

H_max=1000;
H_min=300;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create Stratification
dz=5;
Nz=H_max/dz;
N2=10*omega*ones([Nx+1 Nz]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create topography
s=1.3*omega/sqrt(N2(1)); % Slope
x2=(H_max-H_min)/s;    

% Topography
H0=linspace(H_max,H_min,Nx+1)';
X0=linspace(-x2,0,Nx+1)';X0=X0(1:end-1);

% Define output x- and z-grid
x=linspace(-x2*4,2*x2,400)';
z=linspace(0,H_max,100)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Option 1: Force with an incoming mode-1 wave
if 0

% Set incoming wave amplitudes [A(1) is barotropic wave, A(2) is mode-1, etc.]
A0=zeros([Nm 1]);A0(2)=1;
B0=zeros([Nm 1]);

% Call CELT model
[u p h A B c k uA uB pA pB]=CELT(H0,X0,A0,B0,N2,omega,f,nu,x,z);

else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Option 2: Force with a barotropic transport
    
Q=100; %Barotropic volume flux

% Set incoming wave amplitudes [A(1) is now mode-1, A(2) is mode-2, etc.]
A0=zeros([Nm 1]);
B0=zeros([Nm 1]);

% Call CELT model
[u p h A B c k uA uB pA pB]=CELT_Q(H0,X0,A0,B0,N2,omega,f,nu,x,z,100);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute fluxes (note factor of 1/2 is because A0 and B0 are complex
% harmonic amplitudes.)

% Incident energy flux
Finc=nansum(1/2*rho0*H0(1)*c(:,1).*abs(A0).^2);

% Right going energy flux
FA=1/2*rho0*repmat(H0',[Nm 1]).*c.*abs(A).^2;

% Left going energy flux
FB=1/2*rho0*repmat(H0',[Nm 1]).*c.*abs(B).^2;

% Conservation of energy 
% residual = energy flux in - energy flux out
res=nansum(FA(:,1)+FB(:,end))-nansum(FB(:,1)+FA(:,end));
disp(['The difference between the incoming and outgoing energy fluxes is ',num2str(res),' W/m']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot
ulims=[-1 1]*1/4*max(max(abs(u)));
xlims=[min(x) max(x)]*1e-3;

close;
close(figure(1));figure(1);clf;colormap(gray)

% Velocity field 
pcolor(x*1e-3,z,real(u)');axis ij;shading interp;caxis(ulims);hold on;
set(gca,'xlim',xlims,'xtick',-500:100:500,'ytick',0:200:1000)
xlabel('Distance [km]')
ylabel('Depth [m]')
title('Velocity at t=0')



