function [a,b] = semiinfinite_plate(alpha,beta,gamma,h,N,theta,nu) % [a,b] = semiinfinite_plate(alpha,beta,gamma,h,N,theta) % % this program solves the semiinfinite plate problem by an % eigenfunction matching method. % alpha, beta and gamma are non-dimensionnal parameters depending on the % plate (details can be found on the wikiwave site). % h is the water depth, N is the number of real roots taken when solving % the dispersion equations, and theta is the wave incident angle. % % 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 for a free surface % % b is the vector of coefficients for the eigenfunctions % cos(kappa_n(z+h)/cos(kappa_n h) where kappa_n are the roots of the % dispersion equation for a floating plate % % % the boundary value problem is described in % http://www.wikiwaves.org/index.php/Eigenfunction_Matching_for_a_Semi-Infinite_Floating_Elastic_Plate if nargin == 6 nu = 0.3; end % first we calculate the roots of the dispersion equation for a free % surface. N is the number of real roots. k = dispersion_free_surface(alpha,N,h); % we need to determine the wavenumber in the y direction k_y = sin(theta)*k(1); % Then we calculate the roots of the dispersion equation for a floating % plate kappa = dispersion_elastic_surface(alpha,beta,gamma,N+2,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+3); 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+3; B(j,l) = innerproduct_h_to_0(k(j),kappa(l),h)/... (cos(kappa(l)*h)*cos(k(j)*h)); end end % We define the different blocks of the matrix M1 = A; M2 = -B; M3 = zeros(2,N+1); M4 = zeros(2,N+3); for j = 1:N+3 M4(1,j) = kappa(j)*(kappa_hat(j)^3+kappa_hat(j)*k_y^2*(2-nu))*tan(kappa(j)*h); M4(2,j) = kappa(j)*(kappa_hat(j)^2+k_y^2*nu)*tan(kappa(j)*h); end M5 = zeros(N+1,N+1); M6 = zeros(N+1,N+3); for j = 1:N+1; for l = 1:N+3; M6(j,l) = (k_hat(j) + kappa_hat(l))*B(j,l); end end % We now assemble the full matrix defining the linear system to solve M = [M1,M2;M3,M4;M5,M6]; % We can define the right hand term of the equation as well RHT = zeros(2*N+4,1); RHT(1) = -A(1,1); RHT(N+4) = 2*k_hat(1)*A(1,1); % Finally we can solve the linear system M*S=RHT S = M\RHT; a = S(1:N+1); b = S(N+2:2*N+4); %-------------------------------------------------------------------------- % a set of commands used to test this program. We plot the potential and % its derivative on each side and show that they match. % 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)); % % for count = 1:N+1 % phi_minus = phi_minus + a(count)*cos(k(count)*(z+h))/cos(k(count)*h); % end % for count = 1:N+3 % phi_plus = phi_plus + b(count)*cos(kappa(count)*(z+h))/... % cos(kappa(count)*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)); % % 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); % end % for count = 1:N+3 % phi_plus_n = phi_plus_n + b(count)*kappa_hat(count)*... % cos(kappa(count)*(z+h))/cos(kappa(count)*h); % end % % plot(z,abs(phi_minus_n),z,abs(phi_plus_n)) % % % energy conservation % % D=(kappa_hat(3)*k(1)*cos(k(1)*h)^2/(k_hat(1)*kappa(3)*cos(kappa(3)*h)^2))... % *(((beta/alpha)*4*kappa(3)^4*(-kappa(3))*sin(kappa(3)*h)^2+... % 0.5*sin(2*(kappa(3))*h)+kappa(3)*h)/(0.5*sin(2*k(1)*h)+k(1)*h)); % % abs(a(1))^2 + D*abs(b(3))^2 % this should be one %-------------------------------------------------------------------------- % 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)/h) > 1e-6 % out = 1/2 * (sin(- kappa*h + k*h)) / (k - kappa) ... % + 1/2 * (sin(kappa*h + k*h)) / (kappa + k); out = (k*sin(k*h)*cos(kappa*h)-kappa*cos(k*h)*sin(kappa*h))/... (k^2-kappa^2); else if k == 0 out = h-d; else out = (1/2)*((cos(k*h)*sin(k*h) + k*h)/k); end end