function a = wavemaker(f,alpha,N,h)
% wavemaker(f,z,h,alpha)
%
% program to solve the wavemaker problem. f is the displacement of the
% paddle assumed to be evenly spaced with f(1) the value at z = -h and
% f(end) the value at z = 0.
% alpha is the frequency squared, N+1 is the number of vertical
% eigenfunctions to use, and h is the water depth%
%
% a is the vector of coefficients for the eigenfunctions
% cos(k_n(z+h)/cos(k_n h)
%
% the boundary value problem is described in 
% http://www.wikiwaves.org/index.php/Wavemaker_Theory

% first we solve for the roots of the dispersion equation.
k = dispersion_free_surface(alpha,N,h);

a = zeros(N+1,1); % vector which will store the coefficents for the eigenfunctions

% we calculate the values of z corresponding to f.
z = -h:h/(length(f) -1):0;

% we need the weight for the trapazoidal rule used to integrate
weights = h/(length(f) -1)*ones(size(f));weights(1) = weights(1)/2;weights(end) = weights(end)/2;

for count = 1 : N+1
    eigenfunction = cos(k(count)*(z+h))/cos(k(count)*h);
    integral = sum(weights.*f.*eigenfunction);
    a(count) = -integral/(k(count)*(innerproduct_h_to_d(k(count),k(count),h,0))/cos(k(count)*h)^2);
end

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

% a set of commands used to test this program. Just solves for the
% coefficients and shows that they reproduce (approximately) the input f. 

% h=3;alpha = 2;N=20;
% z = -h:h/200:0;
% f = (z+h).^2;%f = cos(k(1)*(z+h))/cos(k(1)*h); 
% 
% a = wavemaker(f,alpha,N,h);
% 
% k = dispersion_free_surface(alpha,N,h);
% 
% f_approx = 0;
% for count = 1:N+1
%     f_approx = f_approx - k(count)*a(count)*cos(k(count)*(z+h))/cos(k(count)*h);
% end
% 
% plot(z,f,z,f_approx)

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

% this is a program to calculate the inner product of the eigenfunctions
% with themselves. 
   
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