function [a,b] = semiinfinite_plate(alpha,beta,gamma,h,N,theta,nu)
% [a,b] = semiinfinite_plate(alpha,beta,gamma,h,N,theta)
% 
% this program solves the semiinfinite plate problem by an
% eigenfunction matching method. 
% alpha, beta and gamma are non-dimensionnal parameters depending on the
% plate (details can be found on the wikiwave site).
% h is the water depth, N is the number of real roots taken when solving
% the dispersion equations, and theta is the wave incident angle.
%
% a is the vector of coefficients for the eigenfunctions
% cos(k_n(z+h)/cos(k_n h) where k_n are the roots of the dispersion
% equation for a free surface
%
% b is the vector of coefficients for the eigenfunctions
% cos(kappa_n(z+h)/cos(kappa_n h) where kappa_n are the roots of the 
% dispersion equation for a floating plate
%
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Semi-Infinite_Floating_Elastic_Plate


if nargin == 6
    nu = 0.3;
end

% first we calculate the roots of the dispersion equation for a free
% surface. N is the number of real roots.
k = dispersion_free_surface(alpha,N,h); 

% we need to determine the wavenumber in the y direction
k_y = sin(theta)*k(1);

% Then we calculate the roots of the dispersion equation for a floating
% plate
kappa = dispersion_elastic_surface(alpha,beta,gamma,N+2,h);

% we now find the modified coefficients
k_hat = sqrt(k.^2 - k_y^2);
kappa_hat = sqrt(kappa.^2 - k_y^2);

% we calculate the matrix of inner products. 
A = zeros(N+1,N+1);
B = zeros(N+1,N+3);

for j = 1:N+1;
    A(j,j) = innerproduct_h_to_0(k(j),k(j),h)/cos(k(j)*h)^2;
	for l = 1:N+3; 
		B(j,l) = innerproduct_h_to_0(k(j),kappa(l),h)/...
                 (cos(kappa(l)*h)*cos(k(j)*h));
	end
end


% We define the different blocks of the matrix
M1 = A;
M2 = -B;
M3 = zeros(2,N+1);
M4 = zeros(2,N+3);
for j = 1:N+3
    M4(1,j) = kappa(j)*(kappa_hat(j)^3+kappa_hat(j)*k_y^2*(2-nu))*tan(kappa(j)*h);
    M4(2,j) = kappa(j)*(kappa_hat(j)^2+k_y^2*nu)*tan(kappa(j)*h);
end
M5 = zeros(N+1,N+1);
M6 = zeros(N+1,N+3);
for j = 1:N+1;
	for l = 1:N+3; 
		M6(j,l) = (k_hat(j) + kappa_hat(l))*B(j,l);
	end
end

% We now assemble the full matrix defining the linear system to solve
M = [M1,M2;M3,M4;M5,M6];

% We can define the right hand term of the equation as well
RHT = zeros(2*N+4,1);
RHT(1) = -A(1,1);
RHT(N+4) = 2*k_hat(1)*A(1,1);

% Finally we can solve the linear system M*S=RHT
S = M\RHT;

a = S(1:N+1);
b = S(N+2:2*N+4);

%--------------------------------------------------------------------------

% a set of commands used to test this program.  We plot the potential and
% its derivative on each side and show that they match.  
% uncomment these lines to run the test (except this one)
% z = -h:h/200:0;
% 
% phi_minus = cos(k(1)*(z+h))/cos(k(1)*h); 
% phi_plus = zeros(size(phi_minus));
% 
% for count = 1:N+1
%     phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h);
% end
% for count = 1:N+3
%     phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h))/...
%         cos(kappa(count)*h);
% end
% 
% plot(z,abs(phi_minus),z,abs(phi_plus))
% 
% pause
% 
% phi_minus_n = -k_hat(1)*cos(k(1)*(z+h))/cos(k(1)*h); 
% phi_plus_n = zeros(size(phi_minus));
% 
% for count = 1:N+1
%     phi_minus_n = phi_minus_n + a(count)*k_hat(count)*...
%         cos(k(count)*(z+h))/cos(k(count)*h);
% end
% for count = 1:N+3
%     phi_plus_n = phi_plus_n + b(count)*kappa_hat(count)*...
%         cos(kappa(count)*(z+h))/cos(kappa(count)*h);
% end
% 
% plot(z,abs(phi_minus_n),z,abs(phi_plus_n))
% 
% 
% energy conservation 
% 
% D=(kappa_hat(3)*k(1)*cos(k(1)*h)^2/(k_hat(1)*kappa(3)*cos(kappa(3)*h)^2))...
% *(((beta/alpha)*4*kappa(3)^4*(-kappa(3))*sin(kappa(3)*h)^2+...
% 0.5*sin(2*(kappa(3))*h)+kappa(3)*h)/(0.5*sin(2*k(1)*h)+k(1)*h));
% 
% abs(a(1))^2 + D*abs(b(3))^2 % this should be one

%--------------------------------------------------------------------------

% this is a program to calculate the inner product of the eigenfunctions
% with each other. 

function out = innerproduct_h_to_0(k,kappa,h)
% program to calculate 
% \int_{-h}^{0}\cos \kappa\left( z+h\right) 
% \cos k\left(z+h\right) dz

% first we test for the case when \kappa and k are equal (or nearly equal)
if abs((k-kappa)/h) > 1e-6
%     out = 1/2 * (sin(- kappa*h + k*h)) / (k - kappa) ...
%         + 1/2 * (sin(kappa*h + k*h)) / (kappa + k);
      out = (k*sin(k*h)*cos(kappa*h)-kappa*cos(k*h)*sin(kappa*h))/...
          (k^2-kappa^2);
else
    if k == 0
        out = h-d;
    else
        out = (1/2)*((cos(k*h)*sin(k*h) + k*h)/k);
    end
end