function [a,b] = semiinfinite_dock(alpha,h,theta,N) % [a,b] = semiinfinite_dock(alpha,h,theta,N) % % this program solves the semiinfinite dock problem by an % eigenfunction matching. % alpha is the radian frequency squared % h is the water depth, theta is the angle of incidence % and N+1 is the number of eigenfunctions in each region % % a is 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 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_Semi-Infinite_Dock % 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 now assemble the equations to solve M = zeros(size(B)); for j = 1:N+1; for l = 1:N+1; M(j,l) = (k_hat(j) + kappa_hat(l))*B(j,l); end end RHS = zeros(N+1,1); RHS(1) = A(1,1); b = M\(2*k_hat(1)*RHS); a = A\(B*b - RHS); %-------------------------------------------------------------------------- % 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 = -h:h/200:0; % % phi_minus = cos(k(1)*(z+h))/cos(k(1)*h); % phi_plus = zeros(size(phi_minus)); % % phi_minus = phi_plus; % % for count = 1:N+1 % phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h)); % 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; % % for count = 1:N+1 % phi_minus_n = phi_minus_n + a(count)*k_hat(count)*cos(k(count)*(z+h))/cos(k(count)*h); % phi_plus_n = phi_plus_n - b(count)*kappa_hat(count)*cos(kappa(count)*(z+h)); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) %-------------------------------------------------------------------------- % 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