function w_t = h(t,w,k,lx,ltrunc,epsilon,mu)
% Define the equation  w_t = h(t,w,k,lx,ltrunc,epsilon,mu)
% To solve  u_tt = mu*u_xx - f_x,  we let v = u_t  and solve the
% first-order system
%   u_t = v,  v_t = mu*u_xx - f_x + viscosity

u = w(1:lx).';
v = w(lx+1:2*lx).';

u_t = v;

%v_t = local_heat(u,k);

%v_t = local_burgers(u,k);

%v_t = local_benjamin_ono(u,k,lx);

%v_t = local_burgers_hilbert(u,k,lx);

v_t = mu*dderiv(u,k) - deriv(local_hiz(u,k,lx),k);
if (epsilon>0)
    v_t = v_t + local_visc(v,k,lx,ltrunc,epsilon);
end

u_t = truncate(u_t,lx);
v_t = truncate(v_t,lx);

w_t = [u_t,v_t].';
end

%% Subfunction for heat equation
function u_t = local_heat(u,k)

  u_t = dderiv(u,k);                                    % Heat equation
end

%% Subfunction for burgers
function u_t = local_burgers(u,k)

  u_t = -deriv(0.5*u.*u, k);                            % Burgers
end

%% Subfunction for benjamin-ono
function u_t = local_benjamin_ono(u,k,lx)

  u_t = -deriv(0.5*u.*u,k) - hilbert(dderiv(u,k),lx);  % Benjamin-Ono
end

%% Subfunction for burgers-hilbert
function u_t = local_burgers_hilbert(u,k,lx)

  u_t = -deriv(0.5*u.*u, k)+hilbert(u,lx);           % Burgers-Hilbert  
end

%% Subfunction for hiz
function u_t = local_hiz(u,k,lx)

  h = hilbert(u,lx);
  u_t = -0.5*dderiv(hilbert(h.*h,lx),k)-h.*dderiv(u,k) ; % HIZ 
end


%% Subfunction for hiz
% Commutator form
function u_t = local_hizc(u,k,lx)

  h = hilbert(u,lx);
  hx = deriv(h,k);
  u_t = h.*mod_deriv(hx,k)- mod_deriv(h.*hx,k); % HIZ 
end

%% Subfunction for hiz -- evaluated as convolution of fft's
function u_t = local_hizb(u,k,lx)
% LOCAL_HIZB computes flux as convolution of fft's
% Very slow
n = lx;
u_hat = fft(u);
u_shift = fftshift(u_hat);
u_pad = [zeros(1,n/2) u_shift(2:n) zeros(1,n/2)];
u_t_hat = zeros(1,n);
lcf = zeros(1,(2*n-1));
for j=2:n
    jj = j-n/2-1;
    for k=1:(2*n-1)
        kk = k-n;
        lcf(k) = lambdah(-jj,kk,jj-kk);
    end
    u_shift = fliplr(circshift(u_pad,[0 -jj]));
    u_t_hat(j) = sum(lcf.*u_pad.*u_shift)/n;
end
u_t_hat = ifftshift(u_t_hat);
u_t = -real(ifft(u_t_hat));
end

function lcf = lambda(k1,k2,k3)
% LAMBDA computes three-wave interaction coefficient
mod = abs(k1) + abs(k2) + abs(k3);
if (mod==0)
    lcf=0;
else
    % coefficient for 0.5*[h^2]_{xx} - H[hu_{xx}]
    lcf = 2*abs(k1*k2*k3)/mod;
end
%lcf = 1;
end

function lcf = lambdah(k1,k2,k3)
% LAMBDA computes three-wave interaction coefficient
mod = abs(k1) + abs(k2) + abs(k3);
if (mod==0)
    lcf=0;
else
    % coefficient for 0.5*H[h^2]_{xx} + hu_{xx}
    lcf = 2*1i*sign(k1)*abs(k1*k2*k3)/mod;
end
%lcf = 1;
end

%% Subfunction for viscosity
function visc = local_visc(u,k,lx,ltrunc,epsilon)

visc = fft(u);
visc(1:ltrunc) = 0;
visc(lx-ltrunc+1:lx) = 0;
visc = epsilon*ifft(-(k.^2).*visc);
end