function [nep,dnep_dx] = eval_nep(x,g,v,model,gaussgrid)
% function [nep,dnep_dx] = eval_nep(x,g,v,model,whichgrid)
% #codegen
use_core_spline_base = hmode_spline_base_active(v,model);

if ~gaussgrid;
  gauss = false;
  ngrid = model.rgrid.nrho;
else
  gauss = true;
  ngrid = model.rgrid.nrhogauss;
end

switch model.ne.method
  case 'state'
    % ne is a state that is solved for
    if ~use_core_spline_base
      if gauss
        dnep_dnehat = model.ne.Lampgauss;
      else
        dnep_dnehat = model.ne.Lamp;
      end
      nep = dnep_dnehat*x(model.ne.xind,:);
      dnep_dx = zeros(ngrid,numel(x));
      dnep_dx(:,model.ne.xind) = dnep_dnehat;
    else
      if size(v,2)>1 
        if any(v(model.hmode.vind.activation, :)) && ~all(v(model.hmode.vind.activation, :))
          error('eval_nep for time dependent V not implemented for varying spline basis')
        end
      end
      if gauss
        dnep_dxx = model.ne.Lampgauss_hmode;
      else
        dnep_dxx = model.ne.Lamp_hmode;
      end
      xx = [x(model.ne.xind,:); v(model.hmode.vind.ne_rhoedge,:)./model.ne.scal];
      nep = dnep_dxx*xx;
      dnep_dnehat = dnep_dxx(:, 1:end-1);
      dnep_dx(:,model.ne.xind) = dnep_dnehat;
    end
  case 'direct'
    % ne is directly specified by v
    if gauss
      nep = model.ne.Lampgauss*v(model.ne.vind,:);
    else
      nep = model.ne.Lamp*v(model.ne.vind,:);
    end
    dnep_dx = zeros(ngrid,numel(x));
  otherwise
    error('model.ne.method = %s not implemented',model.ne.method)
end

%% check [TO REMOVE]
[nep2,dnep_dx2] = eval_ne(x,g,v,model,gaussgrid,1);
assert(norm(nep2-nep)<1e-6)
assert(norm(dnep_dx2-dnep_dx)<1e-6)

return
