function [tscaled,zeta_pulse,eta_s,u_s]=shelfresponsekdv(to,tf,dt,rho1,rho2,eta_i,S,L,H1,H)
%% A function to calculate the shelf response to an incident KdV soliton impulse.
%
% *Edwin Alfonso-Sosa
% *Ocean Physics Education Inc.
% *Puerto Rico
% *2013
%
%**************************************************************************************************************************************************************************************
% Arguments:
%   to:             Initial Time in minutes
%   tf:             Final Time in minutes
%   dt:             Time Interval in minutes
%   rho1:           Upper layer density kg/m^3
%   rho2:           Lower layer density kg/m^3
%   eta_i:          Amplitude of the Negative KdV soliton in meters. Please insert absolute value.
%   S:              Shelf Depth in meters
%   L:              Length of Shelf in meters
%   H1:             Deep-Ocean Upper Layer Depth in meters
%   H:              Deep-Ocean Total Depth in meters
%**************************************************************************************************************************************************************************************
% Output:
%   tscaled:        Scaled Time
%   zeta_pulse:     Scaled KdV Soliton Pulse Function
%   eta_s:          Scaled Surface displacement at the coast (x=-1)
%   u_s:            Scaled Cross-Shelf Current at the shelfbreak (x=0)
%**************************************************************************************************************************************************************************************
% Information and Reference:
%Calculates the Shelf Response to an incident KdV soliton pulse. Modified version of the analytical model presented in the article: A Model for the Generation of Coastal Seiches
%by Deep-Sea Internal Waves. J. Phys. Oceanogr.,20,1459-1467.doi: http://dx.doi.org/10.1175/1520-0485(1990)020<1459:AMFTGO>2.0.CO;2 The authors explain their model using these words:
%"The incident internal pulse at the interface generates a surface pulse at the shelf break with amplitude Q0/2 which travels across the shelf toward the coast. It's amplitude doubles
%to Q0 at the coast where it reflects from the coastal wall and travels back across the shelf toward the shelf break. Upon reaching the shelf break, part of the pulse is reflected
%shoreward while part continues into the deep ocean leaking energy to deep-ocean surface and internal waves." The incident pulse is a negative KdV soliton. The soliton amplitude can
%be set. The shelf depth and length, the deep-ocean upper layer depth, lower layer depth, total depth and respective densities can be changed to typical measured values.
%Restriction: the shelf depth is less than the deep-ocean upper layer depth: S < H1. To get the dimensional values of eta_s and zeta_pulse multiply each by eta_i; for u_s multiply
%by sqrt(g/H)*eta_i.
%The m-file was written for Matlab R12.1 by Edwin Alfonso-Sosa, Ocean Physics Education, Inc., 2013.
%**************************************************************************************************************************************************************************************
% CODE
%Scaling of Variables
t=[to:dt:tf]';%Time in minutes
H2=H-H1;%Lower Layer Depth in meters
h1_scaled=H1/H;%upper layer scaled
h2_scaled=H2/H;%lower layer scaled
d=S/H;%Scaled depth
g=9.81;%acceleration of gravity in m/s^2
epsilon=(rho2-rho1)/rho2%density difference between layers. No units.
CsquareMinus=0.5-0.5*sqrt(1-(4*epsilon*h1_scaled*h2_scaled));%Csquare-
CsquarePlus=0.5+0.5*sqrt(1-(4*epsilon*h1_scaled*h2_scaled));%Csquare+
mu_1=1/(1-(h1_scaled/CsquareMinus));%Equation 7
N=0;%N=0 means fundamental mode of oscillation, N=1,2,3,... means first, second, third modes of oscillation
omega_peak=sqrt(d)*pi*(.5+N);%Peak Response frequency ? (scaled)
Friction=0.001; %in m/s
r=(Friction/S)/(sqrt(g*H)/L);%Friction Coefficient Scaled
F=1-(r/omega_peak)*i;%Friction factor is an imaginary number.
l=0;%alongshelf wavenumber. If zero then wavefront is parallel to shelfbreak
alpha=sqrt(((omega_peak^2*F)/d)-l);
alpha_F=alpha/F;
BetaSquarePlus=((omega_peak^2)/CsquarePlus)-(l^2);
BetaSquareMinus=((omega_peak^2)/CsquareMinus)-(l^2);
lambda_numerador=(d/(h1_scaled*sqrt(BetaSquarePlus)))*(1-((sqrt(BetaSquarePlus)/sqrt(BetaSquareMinus))*(CsquareMinus/CsquarePlus)*(((epsilon*h1_scaled)-CsquarePlus)/((epsilon*h1_scaled)-CsquareMinus))));
lambda_denominador=(1-((CsquareMinus/CsquarePlus)*(((epsilon*h1_scaled)-CsquarePlus)/((epsilon*h1_scaled)-CsquareMinus))));
lambda=alpha_F*lambda_numerador/lambda_denominador;% Equation 9c. lambda con alpha/F.
%**************************************************************************************************************************************************************************************
%Steps to calculate Qn. Equation 13.
n=[0:100];
An=((4*mu_1)/(1+lambda));
Bn=((-1).^n);
Cn=((1-lambda)/(1+lambda)).^n;
Qn=An.*Bn.*Cn;
%figure;
%plot(n,Qn);
%**************************************************************************************************************************************************************************************
%Steps to calculate the KdV Soliton pulse function
delta_rho=rho2-rho1;
mean_rho=.5*(rho1+rho2);
c=sqrt(((g*delta_rho*H1*H2)/(rho2*H1+rho1*H2)))%linear phase speed in m/s
alfa=((3*c)/(2*H1*H2))*((rho2*H1^2-rho1*H2^2)/(rho2*H1+rho1*H2))%alfa units: 1/s
beta=((c*H1*H2)/6)*((rho1*H1+rho2*H2)/(rho2*H1+rho1*H2))%beta units: m^3/s
V=c+((alfa*(-eta_i))/3)%nonlinear soliton velocity in m/s. eta_i is the amplitude of the negative soliton
delta_square=(12*beta)/(alfa*(-eta_i));%unit: m^2
delta=sqrt(delta_square)%Characteristic Width of the soliton. delta unit: m
tscaled=(t.*60)./(L/sqrt(9.81*H));%time scaled
gamma=(delta/L);%scaled soltiton Characteristic Width (delta)
C_=V/sqrt(g*H);%scaled nonlinear soliton velocity. Check this later
argument=(0-C_.*tscaled)/gamma;
zeta_pulse=-sech(argument).^2;%soliton solution to KdV Equation. Unit amplitude.
%**************************************************************************************************************************************************************************************
%Steps to calculate the surface displacement eta_s. Equation 14a Chapman and Giese (1990)
p=length(tscaled);
for j=1:p,
    x=-1;%x=-1 at the coast
    term1=tscaled(j,1)+((x+1)/sqrt(d))-((2.*n+1)/sqrt(d)); %x=-1 at the coast
    term2=tscaled(j,1)-((x+1)/sqrt(d))-((2.*n+1)/sqrt(d));
    argument1=(0-C_.*term1)/gamma;
    zeta_pulse1=-sech(argument1).^2;
    argument2=(0-C_.*term2)/gamma;
    zeta_pulse2=-sech(argument2).^2;
    eta_s(j,1)=0.5*sum(Qn.*(zeta_pulse1+zeta_pulse2));
end
%**************************************************************************************************************************************************************************************
%Steps to calculate the currents at the shelf break u_s. Equation 14b Chapman and Giese (1990)
for j=1:p,
    x=0;%x=0 at the shelf break
    term1=tscaled(j,1)+((x+1)/sqrt(d))-((2.*n+1)/sqrt(d)); %x=-1 at the coast
    term2=tscaled(j,1)-((x+1)/sqrt(d))-((2.*n+1)/sqrt(d));
    argument1=(0-C_.*term1)/gamma;
    zeta_pulse1=-sech(argument1).^2;
    argument2=(0-C_.*term2)/gamma;
    zeta_pulse2=-sech(argument2).^2;
    u_s(j,1)=(-1/(2*sqrt(d)))*sum(Qn.*(zeta_pulse1-zeta_pulse2));
end
%**************************************************************************************************************************************************************************************
%Plot scaled variables (no units)
figure;
subplot(3,1,1);
plot(tscaled,eta_s);
title('Time series of sea surface elevation at the coast (x=-1)');
ylabel('Sea Surface (scaled) x=-1');
subplot(3,1,2);
plot(tscaled,u_s);
title('Time series of cross-shelf velocity at the shelf break (x=0)');
ylabel('Velocity (scaled) x=0');
subplot(3,1,3)
plot(tscaled,zeta_pulse);
title('Interface displacement at the shelf break (x=0) KdV soliton');
ylabel('Displacement (scaled)');
xlabel('Time (scaled)');
%**************************************************************************************************************************************************************************************
%This m-file script is a courtesy of Ocean Physics Education, Inc., 2013. http://oceanphysics.weebly.com/