function [a,b,c,d] = finite_dock(alpha,h,L,theta,N)
% [a,b,c,d] = finite_dock(alpha,h,L,theta,N)
% 
% this program solves the finite dock problem by an
% eigenfunction matching. The floe occupies the region -L to L
% alpha is the radian frequency squared
% h is the water depth and N + 1 is the number of modes.
% theta is the angle of incidence
% the inner product is calculated using (j-1)*pi/H rather 
% than the water roots
%
% a and d 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
%
% b and c is the vector of coefficients for the eigenfunctions
% cos(kappa_n(z+h)) where kappa_n = n*pi/h
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Finite_Dock
%
% this program was written by Michael Meylan 20/4/08
% http://www.wikiwaves.org/index.php/Michael_Meylan

% 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
kappa = [0:N]*pi/h;

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

% if theta == 0
%     theta = 0.001;
%     disp('We have moved theta from zeros to a close value')
% end
% this problem could be fixed by making the eigenfunctions for j= 0/1
% not exp(i 0 x) and exp(-i 0 x) but 1 and x but I have not 
% done so. 

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

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+1; 
		B(j,l) = innerproduct_h_to_0(k(j),kappa(l),h)/cos(k(j)*h);
	end
end

% we do not calculate in the most efficient manner but we want the code to
% be easy to follow. We now assemble the equations to solve

M = zeros(4*N+4,4*N+4);

% we enter the equations in four groups corresponding to each matching

% first group
M(1:N+1,1:N+1) = -A;
M(1:N+1,N+2:2*N+2) = B;
M(1:N+1,2*N+3:3*N+3) = B*diag(exp(-2*L*kappa_hat));

% second group
M(N+2:2*N+2,1:N+1) = -A*diag(k_hat);
M(N+2:2*N+2,N+2:2*N+2) = -B*diag(kappa_hat);
M(N+2:2*N+2,2*N+3:3*N+3) = B*diag(kappa_hat.*exp(-2*L*kappa_hat));

% third group
M(2*N+3:3*N+3,3*N+4:4*N+4) = -A;
M(2*N+3:3*N+3,2*N+3:3*N+3) = B;
M(2*N+3:3*N+3,N+2:2*N+2) = B*diag(exp(-2*L*kappa_hat));

% fourth group
M(3*N+4:4*N+4,3*N+4:4*N+4) = A*diag(k_hat);
M(3*N+4:4*N+4,2*N+3:3*N+3) = B*diag(kappa_hat);
M(3*N+4:4*N+4,N+2:2*N+2) = -B*diag(kappa_hat.*exp(-2*L*kappa_hat));

if theta == 0
    % first group
    M(1:N+1,2*N+3) = zeros(size(M(N+2:2*N+2,N+2)));
    % second group
    M(N+2:2*N+2,N+2) = -B(:,1)/(2*L);
    M(N+2:2*N+2,2*N+3) = B(:,1)/(2*L);
    % third group
    M(2*N+3:3*N+3,N+2) = zeros(size(M(2*N+3:3*N+3,2*N+3)));;
    % forth group
    M(3*N+4:4*N+4,2*N+3) = B(:,1)/(2*L);
    M(3*N+4:4*N+4,N+2) = -B(:,1)/(2*L);
end
    
    
    
RHS = zeros(4*N+4,1);
RHS(1) = A(1,1);
RHS(N+2) = -k_hat(1)*A(1,1);

sol = M\RHS;
a = sol(1:N+1);
b = sol(N+2:2*N+2);
c = sol(2*N+3:3*N+3);
d = sol(3*N+4:4*N+4);


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

% a set of commands used to test this program.  We plot the potential and
% its derivative on each side of x = -L and show that they match.  You can also check
% that the program is correct by making sure abs(a(1))^2 + abs(d(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));
% % phi_minus = phi_plus;
% 
% count = 1;
% if theta == 0
%     phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h);
%     phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h));
% else
%     phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h);
%     phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h)) + c(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count));
% end
% 
% for count = 2:N+1
%     phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h);
%     phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h)) + c(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count));
% 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;
% 
% count = 1;
% if theta == 0
%     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)/(2*L)*cos(kappa(count)*(z+h)) + c(count)/(2*L)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count))
% else
%     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)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count))
% end
% 
% for count = 2: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)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count));
% 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_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) > 1e-6
    out = 1/2 * (sin(- kappa*h + k*h)) / (k - kappa) ...
        + 1/2 * (sin(kappa*h + k*h)) / (kappa + k);
else
    if k == 0
        out = h;
    else
        out = (1/2)*((cos(k*h)*sin(k*h) + k*h)/k);
    end
end