function [a,b] = vertical_fixed_plate(alpha,h,d1,d2,theta,N)
% [a,b] = vertical_fixed_plate(alpha,h,d1,d2,theta,N)
% 
% this program solves the vertical fixed plate problem by an
% eigenfunction matching. 
% alpha is the radian frequency squared
% h is the water depth, theta is the angle of incidence 
% the plate is form -d1 < z < d2
% and N+1 is the number of eigenfunctions in each region
%
% 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 on the left of the plate
%
% b 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 on the right of the plate
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Vertical_Fixed_Plate

% first we calculate the roots of the dispersion equation for the water
k = dispersion_free_surface(alpha,N,h); 

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

% we calculate the roots of the dispersion equation for the dock covered
% region

k_hat = [sqrt(k.^2 + k_y^2)];


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

for j = 1:N+1;
	for l = 1:N+1; 
		A(j,l) = innerproduct_h_to_d(k(j),k(l),h,d2)/(cos(k(j)*h)*cos(k(l)*h)) ...
            + innerproduct_d1_to_0(k(j),k(l),h,d1)/(cos(k(j)*h)*cos(k(l)*h)) ...
             + k_hat(l)* innerproduct_d1_to_0(k(j),k(l),h,d2)/(cos(k(j)*h)*cos(k(l)*h)) ...
             - k_hat(l)* innerproduct_d1_to_0(k(j),k(l),h,d1)/(cos(k(j)*h)*cos(k(l)*h));
         % this last line is becasue we only want to go from d2 to d1
	end
end

% We now assemble the equations to solve

f = zeros(N+1,1);

for j = 1:N+1
      l=1;
		f(j) = -innerproduct_h_to_d(k(j),k(l),h,d2)/(cos(k(j)*h)*cos(k(l)*h)) ...
            - innerproduct_d1_to_0(k(j),k(l),h,d1)/(cos(k(j)*h)*cos(k(l)*h)) ...
             + k_hat(l)* innerproduct_d1_to_0(k(j),k(l),h,d2)/(cos(k(j)*h)*cos(k(l)*h)) ...
             - k_hat(l)* innerproduct_d1_to_0(k(j),k(l),h,d1)/(cos(k(j)*h)*cos(k(l)*h));
          % this last line is becasue we only want to go from d2 to d1
end

LHS = A\f;

a = LHS/2;
a(1) = a(1) + 1/2;

b =  -LHS/2;
b(1) = b(1) + 1/2;


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

% a set of commands used to test this program.  We plot the potential and
% its derivative on each side and show that they match.  You can also check
% that the program is correct by making sure abs(a(1))^2 + abs(b(1))^2 = 1 (from
% conservation of energy).

% 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);
%     phi_plus = phi_plus + b(count)*cos(k(count)*(z+h))/cos(k(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));
% % phi_minus = phi_plus;
% 
% 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);
%     phi_plus_n = phi_plus_n + b(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h);
% end
% 
% plot(z,abs(phi_minus_n),z,abs(phi_plus_n))

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

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

function out = innerproduct_h_to_d(k,kappa,h,d)
% program to calculate 
% \int_{-h}^{-d}\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-d)) > 1e-6
    out = 1/2 * (sin(d*kappa - d*k - kappa*h + k*h)) / (k - kappa) ...
        + 1/2 * (sin(-d*kappa - d*k + kappa*h + k*h)) / (kappa + k);
else
    if k == 0
        out = h-d;
    else
        out = (1/2)*((cos(-k*d+k*h)*sin(-k*d+k*h)-k*d+k*h)/k);
    end
end

function out = innerproduct_d1_to_0(kh,kd,h,d)
% program to calculate 
% \int_{-d}^{0}\cos kh\left( z+h\right) \cos kd\left(z+h\right) dz

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