function [a_s,b_s,c_s,d_s,e_s,f_s,a_a,b_a,c_a,d_a,e_a,f_a] = two_submerged_finite_docks_symmetry(alpha,h,h1,L1,L2,theta,N)
% [a,b,c,d,e,f] = submerged_finite_dock(alpha,h,L,theta,N)
% 
% this program solves the two finite submerged docks problem by an
% eigenfunction matching. The docks occupies the region -L2 to -L1 and L1
% to L2
% alpha is the radian frequency squared
% h is the water depth and h1 is the depth of the dock. 
% 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 f 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 with depth h
%
% b and c the vector of coefficients for the eigenfunctions
% cos(kd_n(z+h1)/cos(kd_n h) where k_n are the roots of the dispersion
% equation with depth h1
%
% d and e 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_Submerged_Finite_Dock
%
% this program was written by Michael Meylan 8/6/08
% http://www.wikiwaves.org/index.php/Michael_Meylan

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

% define L
L = (L2 - L1)/2;

N_above = round(h1/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,h1);

% we need to determine the wavenumber in the y direction
k_y = sin(theta)*kd(1);

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

% we now find the modified coefficients
kh_hat = sqrt(kh.^2 - k_y^2);
kd_hat = sqrt(kd.^2 - k_y^2);
kappa_hat = sqrt(kappa.^2 - k_y^2);


for j = 1:N+1;
    A(j,j) = innerproduct_h_to_d(kh(j),kh(j),h,0)/cos(kh(j)*h)^2;
    % first we do the modes above
	for l = 1:N_above+1; 
		B(j,l) = innerproduct_d_to_0(kh(j),kd(l),h,h1)/(cos(kh(j)*h)*cos(kd(l)*h1));
    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,h1)/cos(kh(j)*h);
	end
end

% we create a vector of the wave numbers for the middle region. 
k_dock = [kd,kappa];
k_dock_hat = [kd_hat,kappa_hat];

% first we solve the symmetric case
M_s = zeros(4*N+4,4*N+4);

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

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

% second group
M_s(N+2:2*N+2,1:N+1) = -A*diag(kh_hat);
M_s(N+2:2*N+2,N+2:2*N+2) = -B*diag(k_dock_hat);
M_s(N+2:2*N+2,2*N+3:3*N+3) = B*diag(k_dock_hat.*exp(-2*L*k_dock_hat));

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

% fourth group
M_s(3*N+4:4*N+4,3*N+4:4*N+4) = A*diag(kh_hat.*tanh(kh_hat*L1));
M_s(3*N+4:4*N+4,2*N+3:3*N+3) = B*diag(k_dock_hat);
M_s(3*N+4:4*N+4,N+2:2*N+2) = -B*diag(k_dock_hat.*exp(-2*L*k_dock_hat));

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

sol = M_s\RHS;
a_s = sol(1:N+1);
b1 = sol(N+2:2*N+2);
b_s = b1(1:N_above+1);
d_s = b1(N_above+2:end);
c1 = sol(2*N+3:3*N+3);
c_s = c1(1:N_above+1);
e_s = c1(N_above+2:end);
f_s = 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_s(1))^2  = 1 (from
% conservation of energy). We only include matching at the first edge, and
% we only consider the potential.  (please feel free to include the
% derivative and the other edge if you would like). 

% uncomment these lines to run the test (except this one)
% z_full = -h:h/200:0;
% z_above = -h1:h/200:0;
% z_below = -h:h/200:-h1;
% 
% phi_minus = cos(kh(1)*(z_full+h))/cos(kh(1)*h);
% phi_plus_above = zeros(size(z_above));
% phi_plus_below = zeros(size(z_below));
% 
% for count = 1:N+1
%     phi_minus = phi_minus + a_s(count)*cos(kh(count)*(z_full+h))/cos(kh(count)*h);
% end
% 
% for count = 1:N_above+1
%     phi_plus_above = phi_plus_above + b_s(count)*cos(kd(count)*(z_above+h1))/cos(kd(count)*h1) + ...
%         c_s(count)*cos(kd(count)*(z_above+h1))/cos(kd(count)*h1)*exp(-2*L*k_dock_hat(count));
% end
% 
% if theta == 0;
%     phi_plus_below = phi_plus_below + d_s(1)*cos(kappa(1)*(z_below+h));
%     for count = 2:N_below+1
%         phi_plus_below = phi_plus_below + d_s(count)*cos(kappa(count)*(z_below+h)) + ...
%             e_s(count)*cos(kappa(count)*(z_below+h))*exp(-2*L*k_dock_hat(count+N_above+1));
%     end
% else
%     for count = 1:N_below+1
%         phi_plus_below = phi_plus_below + d_s(count)*cos(kappa(count)*(z_below+h)) + ...
%             e_s(count)*cos(kappa(count)*(z_below+h))*exp(-2*L*k_dock_hat(count+N_above+1));
%     end
% end
% 
% plot(z_full,abs(phi_minus),z_above,abs(phi_plus_above),z_below,abs(phi_plus_below))



%  we solve the anti-symmetric case
M_a = zeros(4*N+4,4*N+4);

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

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

% second group
M_a(N+2:2*N+2,1:N+1) = -A*diag(kh_hat);
M_a(N+2:2*N+2,N+2:2*N+2) = -B*diag(k_dock_hat);
M_a(N+2:2*N+2,2*N+3:3*N+3) = B*diag(k_dock_hat.*exp(-2*L*k_dock_hat));

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

% fourth group
M_a(3*N+4:4*N+4,3*N+4:4*N+4) = A*diag(kh_hat./tanh(kh_hat*L1));
M_a(3*N+4:4*N+4,2*N+3:3*N+3) = B*diag(k_dock_hat);
M_a(3*N+4:4*N+4,N+2:2*N+2) = -B*diag(k_dock_hat.*exp(-2*L*k_dock_hat));

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

sol = M_a\RHS;
a_a = sol(1:N+1);
b1 = sol(N+2:2*N+2);
b_a = b1(1:N_above+1);
d_a = b1(N_above+2:end);
c1 = sol(2*N+3:3*N+3);
c_a = c1(1:N_above+1);
e_a = c1(N_above+2:end);
f_a = sol(3*N+4:4*N+4);


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

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

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