function disp_time = k_integrate(time,k0,h,disp_freq)
%   Function to integrate over k and return the free-surface elevation at the
%   time specified. Takes into account the values of the
%   k0    -   array of wavenumber values
%   disp_time   -  array of displacement values at all wavenumbers
%   time - time at which we want to calculate the displacement

% program does not work for time zero so move to close time
if time ==0
    time = 1e-6;
end

if h == 'infinity';
    k=1;
    k_up = (k0(k+1) + k0(k))/2;omega_up = sqrt(k_up);
    k_down = (k0(k) + 0)/2;omega_down = sqrt(k_down);
    disp_time  = 2*sqrt(k0(2))*real(disp_freq(k,1) ...
        *(exp(-i*omega_up*time) - exp(-i*omega_down*time))/(-i*time));
    for k = 2:(length(k0)-1)
        k_up = (k0(k+1) + k0(k))/2;omega_up = sqrt(k_up);
        k_down = (k0(k) + k0(k-1))/2;omega_down = sqrt(k_down);
        disp_time  = disp_time + 2*sqrt(k0(k))*real(disp_freq(k,:) ...
            *(exp(-i*omega_up*time) - exp(-i*omega_down*time))/(-i*time));
    end
    k=length(k0);
    k_up = k0(k);omega_up = sqrt(k_up*tanh(k_up*h));
    k_down = (k0(k) + k0(k-1))/2;omega_down = sqrt(k_down);
    disp_time  = disp_time + 2*sqrt(k0(k))*real(disp_freq(k,:)  ...
        *(exp(-1i*omega_up*time) - exp(-1i*omega_down*time))/(-1i*time));
else

    k=1;
    k_up = (k0(k+1) + k0(k))/2;omega_up = sqrt(k_up*tanh(k_up*h));
    k_down = (k0(k) + 0)/2;omega_down = sqrt(k_down*tanh(k_down*h));
    disp_time  = 2*sqrt(k0(2)*tanh(k0(2)*h))/(tanh(k0(2)*h)+k0(2)*h*sech(k0(2)*h)^2)*real(disp_freq(k,1) ...
        *(exp(-1i*omega_up*time) - exp(-1i*omega_down*time))/(-1i*time));
    for k = 2:(length(k0)-1)
        k_up = (k0(k+1) + k0(k))/2;omega_up = sqrt(k_up*tanh(k_up*h));
        k_down = (k0(k) + k0(k-1))/2;omega_down = sqrt(k_down*tanh(k_down*h));
        disp_time  = disp_time + 2*sqrt(k0(k)*tanh(k0(k)*h))/(tanh(k0(k)*h)+k0(k)*h*sech(k0(k)*h)^2)*real(disp_freq(k,:) ...
            *(exp(-1i*omega_up*time) - exp(-1i*omega_down*time))/(-i*time));
    end
    k=length(k0);
    k_up = k0(k);omega_up = sqrt(k_up*tanh(k_up*h));
    k_down = (k0(k) + k0(k-1))/2;omega_down = sqrt(k_down*tanh(k_down*h));
    disp_time  = disp_time + 2*sqrt(k0(k)*tanh(k0(k)*h))/(tanh(k0(k)*h)+k0(k)*h*sech(k0(k)*h)^2)*real(disp_freq(k,:)  ...
        *(exp(-i*omega_up*time) - exp(-i*omega_down*time))/(-i*time));
end