function [a_s,b_s,c_s,d_s,a_a,b_a,c_a,d_a] = two_finite_docks_symmetry(alpha,h,L1,L2,theta,N) % [a_s,b_s,c_s,d_s,a_a,b_a,c_a,d_a] = two_finite_docks_symmetry(alpha,h,L1,L2,theta,N) % % this program solves the two symmetric finite dock 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 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 d the vector of coefficients for the eigenfunctions % cos(k_n(z+h))/cos(k_n h) where k_n are the roots of the dispersion % equation % % b and c is the vector of coefficients for the eigenfunctions % cos(kappa_n(z+h)) where kappa_n = n*pi/h % % the subscript s and a refer to symmetric and antisymmetric modes % % the boundary value problem is described in % http://www.wikiwaves.org/index.php/Two_Identical_Docks_using_Symmetry % % this program was written by Michael Meylan 22/4/08 % http://www.wikiwaves.org/index.php/Michael_Meylan % define L L = (L2 - L1)/2; % first we calculate the roots of the dispersion equation for the water k = dispersion_free_surface(alpha,N,h); % we need to determine the wavenumber in the y direction k_y = sin(theta)*k(1); % we calculate the roots of the dispersion equation for the dock covered % region kappa = [0:N]*pi/h; % we now find the modified coefficients k_hat = [sqrt(k.^2 - k_y^2)]; kappa_hat = sqrt(kappa.^2 - k_y^2); % 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 % 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 % first we find the symmetric solution 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*kappa_hat)); % second group M_s(N+2:2*N+2,1:N+1) = -A*diag(k_hat); M_s(N+2:2*N+2,N+2:2*N+2) = -B*diag(kappa_hat); M_s(N+2:2*N+2,2*N+3:3*N+3) = B*diag(kappa_hat.*exp(-2*L*kappa_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*kappa_hat)); % fourth group M_s(3*N+4:4*N+4,3*N+4:4*N+4) = A*diag(k_hat.*tanh(k_hat*L1)); M_s(3*N+4:4*N+4,2*N+3:3*N+3) = B*diag(kappa_hat); M_s(3*N+4:4*N+4,N+2:2*N+2) = -B*diag(kappa_hat.*exp(-2*L*kappa_hat)); if theta == 0 % first group M_s(1:N+1,2*N+3) = zeros(size(M_s(N+2:2*N+2,N+2))); % second group M_s(N+2:2*N+2,N+2) = -B(:,1)/(2*L); M_s(N+2:2*N+2,2*N+3) = B(:,1)/(2*L); % third group M_s(2*N+3:3*N+3,N+2) = 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+3) = B(:,1)/(2*L); M_s(3*N+4:4*N+4,N+2) = -B(:,1)/(2*L); end RHS = zeros(4*N+4,1); RHS(1) = A(1,1); RHS(N+2) = -k_hat(1)*A(1,1); sol = M_s\RHS; a_s = sol(1:N+1); b_s = sol(N+2:2*N+2); c_s = sol(2*N+3:3*N+3); d_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(1)) = 1 (from % conservation of energy). % uncomment these lines to run the test (except this one) % z = -h:h/200:0; % % % we begin my comparing the potential and normal derivative at x = -L2 % % phi_minus = cos(k(1)*(z+h))/cos(k(1)*h); % phi_plus = zeros(size(phi_minus)); % % phi_minus = phi_plus; % % count = 1; % if theta == 0 % phi_minus = phi_minus + a_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_s(count)*cos(kappa(count)*(z+h)); % else % phi_minus = phi_minus + a_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_s(count)*cos(kappa(count)*(z+h)) + c_s(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus = phi_minus + a_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_s(count)*cos(kappa(count)*(z+h)) + c_s(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus),z,abs(phi_plus)) % % pause % % phi_minus_n = -k_hat(1)*cos(k(1)*(z+h))/cos(k(1)*h); % phi_plus_n = zeros(size(phi_minus)); % % phi_minus = phi_plus; % % count = 1; % if theta == 0 % phi_minus_n = phi_minus_n + a_s(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_s(count)/(2*L)*cos(kappa(count)*(z+h)) + c_s(count)/(2*L)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % else % phi_minus_n = phi_minus_n + a_s(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus_n = phi_minus_n + a_s(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) % % % we now compare the potential and normal derivative at x = -L1 % pause % % phi_minus = zeros(size(cos(k(1)*(z+h))/cos(k(1)*h))); % phi_plus = zeros(size(phi_minus)); % % count = 1; % if theta == 0 % phi_minus = phi_minus + d_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_s(count)*cos(kappa(count)*(z+h)); % else % phi_minus = phi_minus + d_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_s(count)*cos(kappa(count)*(z+h)) + b_s(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus = phi_minus + d_s(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_s(count)*cos(kappa(count)*(z+h)) + b_s(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus),z,abs(phi_plus)) % % pause % % phi_minus_n = zeros(size(-k_hat(1)*cos(k(1)*(z+h))/cos(k(1)*h))); % phi_plus_n = zeros(size(phi_minus)); % % count = 1; % if theta == 0 % phi_minus_n = phi_minus_n - d_s(count)*k_hat(count)*tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_s(count)/(2*L)*cos(kappa(count)*(z+h)) - b_s(count)/(2*L)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % else % phi_minus_n = phi_minus_n - d_s(count)*k_hat(count)*tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) - b_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus_n = phi_minus_n - d_s(count)*k_hat(count)*tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) - b_s(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) % % pause %-------------------------------------------------------------------------- % now we solve for the antisymmetric case % 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 % first we find the symmetric solution 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*kappa_hat)); % second group M_a(N+2:2*N+2,1:N+1) = -A*diag(k_hat); M_a(N+2:2*N+2,N+2:2*N+2) = -B*diag(kappa_hat); M_a(N+2:2*N+2,2*N+3:3*N+3) = B*diag(kappa_hat.*exp(-2*L*kappa_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*kappa_hat)); % fourth group M_a(3*N+4:4*N+4,3*N+4:4*N+4) = A*diag(k_hat./tanh(k_hat*L1)); M_a(3*N+4:4*N+4,2*N+3:3*N+3) = B*diag(kappa_hat); M_a(3*N+4:4*N+4,N+2:2*N+2) = -B*diag(kappa_hat.*exp(-2*L*kappa_hat)); if theta == 0 % first group M_a(1:N+1,2*N+3) = zeros(size(M_a(N+2:2*N+2,N+2))); % second group M_a(N+2:2*N+2,N+2) = -B(:,1)/(2*L); M_a(N+2:2*N+2,2*N+3) = B(:,1)/(2*L); % third group M_a(2*N+3:3*N+3,N+2) = 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+3) = B(:,1)/(2*L); M_a(3*N+4:4*N+4,N+2) = -B(:,1)/(2*L); end RHS = zeros(4*N+4,1); RHS(1) = A(1,1); RHS(N+2) = -k_hat(1)*A(1,1); sol = M_a\RHS; a_a = sol(1:N+1); b_a = sol(N+2:2*N+2); c_a = sol(2*N+3:3*N+3); d_a = 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. % uncomment these lines to run the test (except this one) % z = -h:h/200:0; % % % we begin by comparing the potential and normal derivative at x = -L2 % % phi_minus = cos(k(1)*(z+h))/cos(k(1)*h); % phi_plus = zeros(size(phi_minus)); % % phi_minus = phi_plus; % % count = 1; % if theta == 0 % phi_minus = phi_minus + a_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_a(count)*cos(kappa(count)*(z+h)); % else % phi_minus = phi_minus + a_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_a(count)*cos(kappa(count)*(z+h)) + c_a(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus = phi_minus + a_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b_a(count)*cos(kappa(count)*(z+h)) + c_a(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus),z,abs(phi_plus)) % % pause % % phi_minus_n = -k_hat(1)*cos(k(1)*(z+h))/cos(k(1)*h); % phi_plus_n = zeros(size(phi_minus)); % % phi_minus = phi_plus; % % count = 1; % if theta == 0 % phi_minus_n = phi_minus_n + a_a(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_a(count)/(2*L)*cos(kappa(count)*(z+h)) + c_a(count)/(2*L)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % else % phi_minus_n = phi_minus_n + a_a(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus_n = phi_minus_n + a_a(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) + c_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) % % % we now compare the potential and normal derivative at x = -L1 % pause % % phi_minus = zeros(size(cos(k(1)*(z+h))/cos(k(1)*h))); % phi_plus = zeros(size(phi_minus)); % % count = 1; % if theta == 0 % phi_minus = phi_minus + d_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_a(count)*cos(kappa(count)*(z+h)); % else % phi_minus = phi_minus + d_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_a(count)*cos(kappa(count)*(z+h)) + b_a(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus = phi_minus + d_a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + c_a(count)*cos(kappa(count)*(z+h)) + b_a(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus),z,abs(phi_plus)) % % pause % % phi_minus_n = zeros(size(-k_hat(1)*cos(k(1)*(z+h))/cos(k(1)*h))); % phi_plus_n = zeros(size(phi_minus)); % % count = 1; % if theta == 0 % phi_minus_n = phi_minus_n + d_a(count)*k_hat(count)/tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_a(count)/(2*L)*cos(kappa(count)*(z+h)) - b_a(count)/(2*L)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % else % phi_minus_n = phi_minus_n + d_a(count)*k_hat(count)/tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) - b_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % for count = 2:N+1 % phi_minus_n = phi_minus_n + d_a(count)*k_hat(count)/tanh(k_hat(count)*L1)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n + c_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h)) - b_a(count)*kappa_hat(count)*cos(kappa(count)*(z+h))*exp(-2*L*kappa_hat(count)); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) % % pause %-------------------------------------------------------------------------- % this is a program to calculate the inner product of the eigenfunctions % with each other. 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