function [Ew_lo] = p12ew(t,e,fp,rho);
% transform surface elevation to wave energy signal
%
% Ad Reniers 05 jun 2000

g = 9.81;

[nt] = length(t);
dt = t(2)-t(1);
T = t(nt);
df = 1/(nt*dt);
f_nyq = 1/2/dt;

% fourier transform of e

ed = detrend(e,0);
ef = fft(ed);

if (nt/2-floor(nt/2)) > eps
   f = [0:df:f_nyq -f_nyq:df:-df]';
else
   f = [0:df:f_nyq -f_nyq+df:df:-df]';
end

f_high = 1.0; %cut off frequency
i = find(abs(f)>f_high);
ni = length(i); 
f_low = 0.01; % high pass filter
j = find(abs(f)<f_low);
nj = length(j);

trans= ones(size(f));

% transfer function 
f(1) = 1e-4;   % to avoid singularities
elev = 3; % distance measured from the bed (positive upward)
h = 4.1; %local water depth at 20 m from the paddle
ht = ones(nt,1)*h;
hz = ones(nt,1)*elev;
k  = disper(abs(f)*2*pi,h);
trans = (cosh(k.*ht)./(cosh(k.*hz)));
% trans = sinh(k.*ht)./(cosh(k.*hz))./(2*pi.*abs(f)); % transfer function for velocities
trans(i,:) = zeros(size(trans(i,:)));
trans(j,:) = zeros(size(trans(j,:)));



ef = ef.*trans;
ed = ifft(ef);
ed = real(ed);

% hilbert transform to obtain envelope

eh = hilbert(ed);
Ew = 0.5*rho*g*abs(eh).^2;

% filter to remove high frequencies

fc    = 0.5*fp;
[eh_lo,eh_hi] = f_flt(abs(eh),fc,dt);
Ew_lo = 0.5*rho*g*abs(eh_lo).^2;

out(:,1) = t;
out(:,2) = Ew_lo;

