function [potential,R,T]=finite_dock_green_function(alpha,h,L,N)
% [potential,R,T]=finite_dock_green_function(alpha,h,L,N)
% program to calculate the modes of a floating elastic plate
%--------------------------------------------------------------------------
% Variables:
% alpha  = frequency squared
% h      =  water depth
% L      =  beam/plate half-length
% N     =  number of calculation points along the plate is 2*N+1 (because
% we use Simpson's rule
%
% xi are the coefficients for the displacement, and displacement and
% potential and the values of the displacement and potential at the 2N+1 points
% evenly spaced between -L and L.
%
% xi = coefficients in the free modes 
% displacement = complex displacement
% potentiall = complex potential
% R = Reflection coefficient
% T = Transmission coefficient
% 
% Details can be found on http://www.wikiwaves.org/index.php/Green_Function_Method_for_a_Finite_Dock

% start of main

np = 2*N+1;                                                  

% define calculation points
x = -L:L/N:L;

dx = (2*L)/(np-1);
simp = 2*dx/3*ones(1,np);simp(2:2:np-1) = 4*dx/3;
simp(1) = dx/3;simp(end) = dx/3;
simp = diag(simp);

% calculate free-surface Green function
G = matrix_G_surface(alpha,h,L,N);

k = dispersion_free_surface(alpha,0,h);
% calculate incident wave potential
phi_i =exp(-k*x);

% solve for the diffraction wave potential
phi_d = inv(eye(length(G))-alpha*G)*phi_i.';

potential =  transpose(phi_d);

% we calculate the reflection and transmission coefficients.

R = -alpha/(tan(k*h) + k*h*sec(k*h)^2)*((exp(-k*x)).*(potential))*diag(simp);
T = 1 - alpha/(tan(k*h) + k*h*sec(k*h)^2)*((exp(k*x)).*(potential))*diag(simp);

% end of main