function [vec_a,vec_b]= circle_plate_matching_one_n(alpha,beta,gamma,nu,h,a,n,N)
% [vec_a,vec_b]= circle_plate_matching_one_n(alpha,beta,gamma,h,a,n,N)
%
% this program solves the circlular plate problem for one
% value of n (the coefficient in the expansion of e^(i n\theta))
%
% alpha is the infinite depth wavenumber
% beta and gamma are coefficients which depend on the plate properties
% nu is the Poisson's ratio
% h is the water depth. 
% a is the radius of the plate. 
% n is an integer.
% N is the number of real roots of the dispersion equation for a free 
% surface.
%
% 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 plate
% 
% 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=dispersion_elastic_surface(alpha,beta,gamma,N+2,h);

% 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 calculate the coefficients which intervene in the boundary conditions
% of the plate equations.

c_1=zeros(1,N+3);
c_2=zeros(1,N+3);

for l=1:N+3
    c_1(l)=(kappa(l)^2*besseli(n,kappa(l)*a)-((1-nu)/a)*...
        (kappa(l)*dbesseli(n,kappa(l)*a)-n^2*besseli(n,kappa(l)*a)/a))/...
        (beta*kappa(l)^4+1-alpha*gamma);
    
    c_2(l)=(kappa(l)^3*dbesseli(n,kappa(l)*a)+((1-nu)*n^2/a^2)*...
        (kappa(l)*dbesseli(n,kappa(l)*a)+besseli(n,kappa(l)*a)/a))/...
        (beta*kappa(l)^4+1-alpha*gamma);
end

% We set up now the matrix of the linear system to solve

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

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

f(1) = besseli(n,k(1)*a)*A(1,1);
f(N+4) = 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)*cos(kappa(1)*(z+h))/cos(kappa(1)*h)*besseli(n,kappa(1)*a);
% 
% 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+3
%     phi_in = phi_in + vec_b(count)*cos(kappa(count)*(z+h))/cos(kappa(count)*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 = vec_b(1)*kappa(1)*cos(kappa(1)*(z+h))/cos(kappa(1)*h)*dbesseli(n,kappa(1)*a);
% 
% 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+3
%     phi_in_n = phi_in_n + vec_b(count)*kappa(count)*cos(kappa(count)*(z+h))/cos(kappa(count)*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)/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) + 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_plate_matching_one_n(alpha,beta,gamma,nu,h,a,n,N)
% % [vec_a,vec_b]= circle_plate_matching_one_n(alpha,beta,gamma,h,a,n,N)
% %
% % this program solves the circlular plate problem for one
% % value of n (the coefficient in the expansion of e^(i n\theta))
% %
% % alpha is the infinite depth wavenumber
% % beta and gamma are coefficients which depend on the plate properties
% % nu is the Poisson's ratio
% % h is the water depth. 
% % a is the radius of the plate. 
% % n is an integer.
% % N is the number of real roots of the dispersion equation for a free 
% % surface.
% %
% % 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 plate
% % 
% % 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=dispersion_elastic_surface(alpha,beta,gamma,N+2,h);
% 
% % 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 calculate the coefficients which intervene in the boundary conditions
% % of the plate equations.
% 
% c_1=zeros(1,N+3);
% c_2=zeros(1,N+3);
% 
% for l=1:N+3
%     c_1(l)=(kappa(l)^2*besseli(n,kappa(l)*a,1)-((1-nu)/a)*...
%         (kappa(l)*dbesseli(n,kappa(l)*a)-n^2*besseli(n,kappa(l)*a,1)/a))/...
%         (beta*kappa(l)^4+1-alpha*gamma);
%     
%     c_2(l)=(kappa(l)^3*dbesseli(n,kappa(l)*a)+((1-nu)*n^2/a^2)*...
%         (kappa(l)*dbesseli(n,kappa(l)*a)+besseli(n,kappa(l)*a,1)/a))/...
%         (beta*kappa(l)^4+1-alpha*gamma);
% end
% 
% % We set up now the matrix of the linear system to solve
% 
% M = zeros(2*N+4);
% for p = 1:N+1
%     M(p,p) = -besselk(n,k(p)*a,1)*A(p,p);
%     M(N+2,p) = 0;
%     M(N+3,p) = 0;
%     M(N+3+p,p) = -k(p)*dbesselk(n,k(p)*a)*A(p,p);
%     for m = 1:N+3;
%          M(p,N+1+m) = besseli(n,kappa(m)*a,1)*B(p,m);
%          M(N+2,N+1+m) = c_1(m);
%          M(N+3,N+1+m) = c_2(m);
%          M(N+3+p,N+1+m) = kappa(m)*dbesseli(n,kappa(m)*a)*B(p,m);
%     end
% end
% 
% f = zeros(2*N+4,1);
% 
% f(1) = besseli(n,k(1)*a,1)*A(1,1);
% f(N+4) = 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)*cos(kappa(1)*(z+h))/cos(kappa(1)*h)*besseli(n,kappa(1)*a);
% % 
% % 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+3
% %     phi_in = phi_in + vec_b(count)*cos(kappa(count)*(z+h))/cos(kappa(count)*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 = vec_b(1)*kappa(1)*cos(kappa(1)*(z+h))/cos(kappa(1)*h)*dbesseli(n,kappa(1)*a);
% % 
% % 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+3
% %     phi_in_n = phi_in_n + vec_b(count)*kappa(count)*cos(kappa(count)*(z+h))/cos(kappa(count)*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;