function [xi,displacement,potential,R,T]=rigid_plate_modes(alpha,gamma,h,L,N,option)
% [xi,displacement,potential,R,T]=rigid_plate_modes(alpha,h,L,N)
% program to calculate the modes of a floating elastic plate
%--------------------------------------------------------------------------
% Variables:
% alpha  = frequency squared
% gamma  = non-dimensional plate mass
% 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.
%
% option is an optional argument, if option = 'heave' then we solve
% for heave only and if option = 'pitch' we solve for pitch only.
%
% 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_Floating_Body_on_the_Surface#Equation_in_Terms_of_the_Modes


np = 2*N+1;


% Functions:
% beam_ev               ->  mode eigenvalues
% beam_em               ->  dry natural modes
% matrix_G_surface      ->  free-surface Green function
%--------------------------------------------------------------------------



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

if (nargin == 5) || ~(strcmp(option, 'heave') || strcmp(option, 'pitch'))
    wn = [1/sqrt(2*L) * ones(size(x)); x/L*sqrt(3/(2*L))];
elseif strcmp(option ,'heave')
    wn = [1/sqrt(2*L) * ones(size(x))];
else
    wn = [x/L*sqrt(3/(2*L))];
end




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

% solve for the plate radiation potentials
phi_r = inv(eye(size(G))-alpha*G)*G*(-i*sqrt(alpha))*wn.';

% calculate added mass and damping coefficients
% define simpsons coeff. matrix
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);

ab = transpose(wn*simp*phi_r);

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.';

% solve for the xi_n's
diff= wn*simp*phi_d; % calculate diffraction matrix
xi = (eye(size(wn,1)) - alpha*gamma*eye(size(wn,1)) + i*sqrt(alpha)*ab)\(-i*sqrt(alpha)*diff); % calculated alpha_n

displacement = transpose(xi)*wn;
potential =  transpose(xi)*transpose(phi_r) + 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 - i/sqrt(alpha)*displacement))*diag(simp);
T = 1 - alpha/(tan(k*h) + k*h*sec(k*h)^2)*((exp(k*x)).*(potential - i/sqrt(alpha)*displacement))*diag(simp);

% end of main