function [a,b] = semi_infinite_change_in_depth(alpha,h,d,N_left,N_right)
% [a,b] = semi_infinite_change_in_depth(alpha,h,d,N_left,N_right)
% 
% this program solves the submerged semiinfinite dock problem by an
% eigenfunction matching. 
% alpha is the radian frequency squared
% h is the water depth on the left (incident wave direction)
% d is the depth of right. Note we assume d > h.
% theta is the angle of incidence 
% and N_left+1 or N_right 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 on the left
%
% b is the vector of coefficients for the eigenfunctions
% cos(kd_n(z+d))/cos(kd_n d) where kd_n are the roots of the dispersion
% equation for depth d (b has size N_above) on the right
%
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Semi-Infinite_Change_in_Depth


       


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


% we calculate the matrix of inner products. 
A_h = zeros(N_left+1,N_left+1);
A_d = zeros(N_right+1,N_right+1);
B = zeros(N_right+1,N_left+1);


for j = 1:N_right+1;
    Ad(j,j) = innerproduct_d_to_0(kd(j),kd(j),d,d)/cos(kd(j)*d)^2;
    % first we do the modes above
	for l = 1:N_left+1; 
        Ah(l,l) = innerproduct_d_to_0(kh(l),kh(l),h,h)/cos(kh(l)*h)^2;
		B(j,l) = innerproduct_d_to_0(kd(j),kh(l),d,h)/(cos(kh(l)*h)*cos(kd(j)*d));
    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 = [-Ah,transpose(B);B*diag(kh),diag(kd)*Ad];

RHS = zeros(N_left+N_right+2,1);
RHS(1) = Ah(1,1);
RHS(N_left+2:end) = kh(1)*B(:,1);

LHS = M\RHS;
a = LHS(1:N_left+1);
b = LHS(N_left+2:end);

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

% 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 
% (tan(kh(1)*h) + kh(1)*h*sec(kh(1)*h)^2)* abs(a(1)).^2 ...
% + (tan(kd(1)*d) + kd(1)*d*sec(kd(1)*d)^2)* abs(b(1)).^2 
% = 
% (tan(kh(1)*h) + kh(1)*h*sec(kh(1)*h)^2)
% (from from conservation of energy).

% uncomment these lines to run the test (except this one)
z_left = -h:h/200:0;
z_right = -d:h/200:0;

phi_left = cos(kh(1)*(z_left+h))/cos(kh(1)*h);
phi_right = zeros(size(z_right));

% phi_minus = phi_plus;

for count = 1:N_left+1
    phi_left = phi_left + a(count)*cos(kh(count)*(z_left+h))/cos(kh(count)*h);
end

for count = 1:N_right+1
    phi_right = phi_right + b(count)*cos(kd(count)*(z_right+d))/cos(kd(count)*d);
end


plot(z_left,abs(phi_left),z_right,abs(phi_right))

pause

phi_left_n = -kh(1)*cos(kh(1)*(z_left+h))/cos(kh(1)*h);
phi_right_n = zeros(size(z_right));


for count = 1:N_left+1
    phi_left_n = phi_left_n + kh(count)*a(count)*cos(kh(count)*(z_left+h))/cos(kh(count)*h);
end

for count = 1:N_right+1
    phi_right_n = phi_right_n - kd(count)*b(count)*cos(kd(count)*(z_right+d))/cos(kd(count)*d);
end


plot(z_left,abs(phi_left_n),z_right,abs(phi_right_n))

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

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

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