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