function [vec_a,vec_b,vec_c] = circle_submerged_dock_matching_one_n(alpha,h,d,a,n,N)
% [vec_a,vec_b,vec_c] = circle_submerged_dock_matching_one_n(alpha,h,d,a,n,N)
% 
% this program solves the submerged semiinfinite circular dock problem by an
% eigenfunction matching. 
% alpha is the radian frequency squared
% h is the water depth, d is the depth of the submerged plate,
% a is the dock radius
% and N+1 is the number of eigenfunctions in each region
%
% a is the vector of coefficients for the eigenfunctions
% cos(kh_n(z+h)/cos(kh_n h) where kh_n are the roots of the dispersion
% equation for depth h
%
% b is the vector of coefficients for the eigenfunctions
% cos(kd_n(z+h)/cos(kd_n h) where kd_n are the roots of the dispersion
% equation for depth d (b has size N_above)
%
% c is the vector of coefficients for the eigenfunctions
% cos(kappa_n(z+h) where kappa_n = n*pi/(h-d) (c has size N_below)
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Submerged_Semi-Infinite_Dock


% we need to allocate the modes on the right so that they are roughly in
% the same ratio as d and (h-d) while allowing at least one. 
% remember that there is actually N + 1 modes.

N_above = round(d/h*(N+1)) -1; % number of modes above is N_above + 1;
N_below = N-1-N_above; % numebr of modes below is N_below +1;
   if N_above <  0
       N_above = 0;
       N_below = N-1;
   elseif N_below < 0
       N_above = N-1;
       N_below = 0;
   end
       


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

% we calculate the roots of the dispersion equation for the dock covered
% region
kappa = [0:N_below]*pi/(h-d);

mu = [kd,kappa];

% 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_d(kh(j),kh(j),h,0)/cos(kh(j)*h)^2;
	for l = 1:N_above+1; 
		B(j,l) = innerproduct_d_to_0(kh(j),kd(l),h,d)/(cos(kh(j)*h)*cos(kd(l)*d));
    end
    % then the modes below
    for l = 1:N_below+1; 
		B(j,l+N_above+1) = innerproduct_h_to_d(kh(j),kappa(l),h,d)/cos(kh(j)*h);
	end
end

% We now assemble the equations to solve


M = zeros(2*N+2);
for p = 1:N+1
    M(p,p) = -besselk(n,kh(p)*a)*A(p,p);
    M(N+1+p,p) = -kh(p)*dbesselk(n,kh(p)*a)*A(p,p);
    m=N_above+2;
    M(p,N+1+m) = B(p,m);
    M(N+1+p,N+1+m) = abs(n)/a*B(p,m);
    for m = [1:N_above+1,N_above+3:N+1];
         M(p,N+1+m) = besseli(n,mu(m)*a)*B(p,m);
         M(N+1+p,N+1+m) = mu(m)*dbesseli(n,mu(m)*a)*B(p,m);
    end
end


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

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

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

vec_c = vec_b_temp(N_above+2:end);vec_b = vec_b_temp(1:N_above+1);
%--------------------------------------------------------------------------

% 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)) = 1 (from
% conservation of energy).

% uncomment these lines to run the test (except this one)
% z_full = -h:h/200:0;
% z_above = -d:h/200:0;
% z_below = -h:h/200:-d;
% 
% phi_out = cos(kh(1)*(z_full+h))/cos(kh(1)*h)*besseli(n,kh(1)*a);
% phi_in_above = zeros(size(z_above));
% phi_in_below = vec_c(1)*ones(size(z_below));
% 
% for count = 1:N+1
%     phi_out = phi_out + vec_a(count)*cos(kh(count)*(z_full+h))/cos(kh(count)*h)*besselk(n,kh(count)*a);
% end
% 
% for count = 1:N_above+1
%     phi_in_above = phi_in_above + vec_b(count)*cos(kd(count)*(z_above+d))/cos(kd(count)*d)*besseli(n,kd(count)*a);
% end
% 
% for count = 2:N_below+1
%     phi_in_below = phi_in_below + vec_c(count)*cos(kappa(count)*(z_below+h))*besseli(n,kappa(count)*a);
% end
% 
% plot(z_full,abs(phi_out),z_above,abs(phi_in_above),z_below,abs(phi_in_below))
% 
% pause
% 
% phi_out_n = kh(1)*cos(kh(1)*(z_full+h))/cos(kh(1)*h)*dbesseli(n,kh(1)*a);
% phi_in_above_n = zeros(size(z_above));
% phi_in_below_n = abs(n)/a*vec_c(1)*ones(size(z_below));
% % phi_out = phi_in;
% 
% for count = 1:N+1
%     phi_out_n = phi_out_n + kh(count)*vec_a(count)*cos(kh(count)*(z_full+h))/cos(kh(count)*h)*dbesselk(n,kh(count)*a);
% end
% 
% for count = 1:N_above+1
%     phi_in_above_n = phi_in_above_n - kd(count)*vec_b(count)*cos(kd(count)*(z_above+d))/cos(kd(count)*d)*dbesseli(n,kd(count)*a);
% end
% 
% for count = 1:N_below+1
%     phi_in_below_n = phi_in_below_n - kappa(count)*vec_c(count)*cos(kappa(count)*(z_below+h))*dbesseli(n,kappa(count)*a);
% end
% 
% plot(z_full,abs(phi_out_n),z_above,abs(phi_in_above_n),z_below,abs(phi_in_below_n))

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

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

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;

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_d_to_0(kh,kd,h,d)
% program to calculate 
% \int_{-d}^{0}\cos kh\left( z+h\right) \cos kd\left(z+d\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)*(((sin(kh*h-d*kd))*kh+(sin(kh*h-d*kd))*kd ...
        + (sin(kh*h+d*kd))*kh -(sin(kh*h+d*kd))*kd-2*(sin(-d*kh+kh*h))*kh)/(kh^2-kd^2));
else
    if kd == 0
        out = d;
    else
        out = (1/4)*((sin(kd*h+d*kd)+2*(cos(-d*kd+kd*h))*d*kd-sin(-d*kd+kd*h))/kd);
    end
end