function [vec_a,vec_b]= circle_dock_matching_one_n(alpha,h,a,n,N)
% this program solves the circlular dock problem for one
% value of n (the coefficient in the expansion of e^(i n\theta)
%
% alpha is the infinite depth wavenumber
% h is the water depth. 
% a is the radius of the dock. 
%
% vec_a is the vector of coefficients in the expansion of the water
% vec_b is the vector of coefficients in the expansion of the dock
% 
% this program was written by Michael Meylan 25/5/08 and the problem is discuss in
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Circular_Dock 
% 
% We also include a version of the code which uses modified bessel
% functions and is more stable at the end (the code is commented)

k=dispersion_free_surface(alpha,N,h); 
kappa = [0:N]*pi/h;

% 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

M = zeros(2*N+2);
for p = 1:N+1
    M(p,p) = -besselk(n,k(p)*a)*A(p,p);
    M(N+1+p,p) = -k(p)*dbesselk(n,k(p)*a)*A(p,p);
    m=1;
    M(p,N+1+m) = B(p,m);
    M(N+1+p,N+1+m) = abs(n)/a*B(p,m);
    for m = 2:N+1;
         M(p,N+1+m) = besseli(n,kappa(m)*a)*B(p,m);
         M(N+1+p,N+1+m) = kappa(m)*dbesseli(n,kappa(m)*a)*B(p,m);
    end
end

f = zeros(2*N+2,1);

f(1) = besseli(n,k(1)*a)*A(1,1);
f(N+2) = k(1)*dbesseli(n,k(1)*a)*A(1,1);

q = M\f;
vec_a = q(1:N+1);
vec_b = q(N+2:end);


% uncomment these lines to run the test (except this one)
% z = -h:h/200:0;
% 
% phi_out = cos(k(1)*(z+h))/cos(k(1)*h)*besseli(n,k(1)*a);
% phi_in = vec_b(1)*ones(size(phi_out));
% 
% for count = 1:N+1
%     phi_out = phi_out + vec_a(count)*cos(k(count)*(z+h))/cos(k(count)*h)*besselk(n,k(count)*a);
% end
% for count = 2:N+1
%     phi_in = phi_in + vec_b(count)*cos(kappa(count)*(z+h))*besseli(n,kappa(count)*a);
% end
% plot(z,abs(phi_out),z,abs(phi_in))
% 
% pause
% 
% phi_out_n = k(1)*cos(k(1)*(z+h))/cos(k(1)*h)*dbesseli(n,k(1)*a);
% phi_in_n = abs(n)/a*vec_b(1)*ones(size(phi_out));
% 
% for count = 1:N+1
%     phi_out_n = phi_out_n + vec_a(count)*k(count)*cos(k(count)*(z+h))/cos(k(count)*h)*dbesselk(n,k(count)*a);
% end
% for count = 2:N+1
%     phi_in_n = phi_in_n + vec_b(count)*kappa(count)*cos(kappa(count)*(z+h))*dbesseli(n,kappa(count)*a);
% end
% plot(z,abs(phi_out_n),z,abs(phi_in_n))

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

function out = dbesseli(n,z)
% this function calcualtes the derivative of besselk(n,z)

out = (besseli(n+1,z) + besseli(n-1,z))/2;

function out = dbesselk(n,z)
% this function calcualtes the derivative of besselk(n,z)

out = -(besselk(n+1,z) + besselk(n-1,z))/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% this is a version of the program using modified bessel function which is
% more stable

% function [vec_a,vec_b]= circle_dock_matching_one_n(alpha,h,a,n,N)
% % this program solves the circlular dock problem for one
% % value of n (the coefficient in the expansion of e^(i n\theta)
% %
% %
% % alpha is the infinite depth wavenumber
% % H is the water depth. 
% % a is the radius of the dock. 
% %
% % vec_a is the vector of coefficients in the expansion of the water
% % vec_b is the vector of coefficients in the expansion of the dock
% % 
% % this program was written by Michael Meylan 25/5/08 and the problem is discuss in
% % http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Circular_Dock
% 
% k=dispersion_free_surface(alpha,N,h); 
% kappa = [0:N]*pi/h;
% 
% % 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
% 
% % note that we are using the modified bessel functions to get greater
% % stability. 
% 
% M = zeros(2*N+2);
% for p = 1:N+1
%     M(p,p) = -besselk(n,k(p)*a,1)*A(p,p);
%     M(N+1+p,p) = -k(p)*dbesselk(n,k(p)*a)*A(p,p);
%     m=1;
%     M(p,N+1+m) = B(p,m);
%     M(N+1+p,N+1+m) = abs(n)/a*B(p,m);
%     for m = 2:N+1;
%          M(p,N+1+m) = besseli(n,kappa(m)*a,1)*B(p,m);
%          M(N+1+p,N+1+m) = kappa(m)*dbesseli(n,kappa(m)*a)*B(p,m);
%     end
% end
% 
% f = zeros(2*N+2,1);
% 
% f(1) = besseli(n,k(1)*a,1)*A(1,1);
% f(N+2) = k(1)*dbesseli(n,k(1)*a)*A(1,1);
% 
% q = M\f;
% vec_a = q(1:N+1);
% vec_b = q(N+2:end);
% 
% % uncomment these lines to run the test (except this one)
% z = -h:h/200:0;
% 
% phi_out = cos(k(1)*(z+h))/cos(k(1)*h)*besseli(n,k(1)*a,1);
% phi_in = vec_b(1)*ones(size(phi_out));
% 
% for count = 1:N+1
%     phi_out = phi_out + vec_a(count)*cos(k(count)*(z+h))/cos(k(count)*h)*besselk(n,k(count)*a,1);
% end
% for count = 2:N+1
%     phi_in = phi_in + vec_b(count)*cos(kappa(count)*(z+h))*besseli(n,kappa(count)*a,1);
% end
% plot(z,abs(phi_out),z,abs(phi_in))
% 
% pause
% 
% phi_out_n = k(1)*cos(k(1)*(z+h))/cos(k(1)*h)*dbesseli(n,k(1)*a);
% phi_in_n = abs(n)/a*vec_b(1)*ones(size(phi_out));
% 
% for count = 1:N+1
%     phi_out_n = phi_out_n + vec_a(count)*k(count)*cos(k(count)*(z+h))/cos(k(count)*h)*dbesselk(n,k(count)*a);
% end
% for count = 2:N+1
%     phi_in_n = phi_in_n + vec_b(count)*kappa(count)*cos(kappa(count)*(z+h))*dbesseli(n,kappa(count)*a);
% end
% plot(z,abs(phi_out_n),z,abs(phi_in_n))
% 
% 
% % now we have to convert back to the standard bessel functions
% vec_a = vec_a.*besselk(n,k*a,1).'./besselk(n,k*a).'*besseli(n,k(1)*a)/besseli(n,k(1)*a,1);
% vec_b(2:N+1) = vec_b(2:N+1).*(besseli(n,kappa(2:N+1)*a,1).')./(besseli(n,kappa(2:N+1)*a).');
% 
% 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);
% else
%     if k == 0
%         out = h-d;
%     else
%         out = (1/2)*((cos(k*h)*sin(k*h) + k*h)/k);
%     end
% end
% 
% function out = dbesseli(n,z)
% % this function calcualtes the derivative of besselk(n,z)
% 
% out = (besseli(n+1,z,1) + besseli(n-1,z,1))/2;
% 
% function out = dbesselk(n,z)
% % this function calcualtes the derivative of besselk(n,z)
% 
% out = -(besselk(n+1,z,1) + besselk(n-1,z,1))/2;