function [ntmres,dntmres_dxdot,dntmres_dx,dntmres_du,dntmres_dd,dntmres_dw] ...
    = ntm_equation(stap,trap,geop,u,model,params)

if strcmp(model.ntm.method,'state')
    if params.ntm.implicit % recompute including jacobians
        error('ntm implicit not yet implemented')
    else % explicit
        dwdt = trap.dwntmdt; % computed from MRE(xk)
        wdot = stap.wdotntms; % present dw/dt from xk1 and xk
        fntm = wdot - dwdt; % residual
        
        dntmres_dx = zeros(model.ntm.n_ntm,model.dims.nx);
        dntmres_dxdot = zeros(model.ntm.n_ntm,model.dims.nx);
        dntmres_dxdot(1:model.ntm.n_ntm,model.ntm.xind) = eye(model.ntm.n_ntm);
        dntmres_du = zeros(model.ntm.n_ntm,model.dims.nu);
        dntmres_dd = zeros(model.ntm.n_ntm,model.dims.nd);
        dntmres_dw = zeros(model.ntm.n_ntm,model.dims.nw);
    end
    ntmres    = model.ntm.eqscal*fntm;
else
    ntmres = [];
    dntmres_dx    = zeros(0,model.dims.nx);
    dntmres_dxdot = zeros(0,model.dims.nx);
    dntmres_du    = zeros(0,model.dims.nu);
    dntmres_dd = zeros(0,model.dims.nd);
    dntmres_dw = zeros(0,model.dims.nw);
end