function [ne,dne_dx] = eval_ne(x,~,v,model,gauss,varargin)
% function [ne,dne_dx] = eval_nep(x,~,model,whichgrid,prime)
% #codegen

if numel(varargin)==0
  prime = 0; % default
else
  prime = varargin{1};
  assert(any(prime==[0,1,2]));
end

if ~gauss
  dne_dx = zeros(model.rgrid.nrho,numel(x)); % init
else
  dne_dx = zeros(model.rgrid.nrhogauss,numel(x)); % init
end
use_core_spline_base = hmode_spline_base_active(v,model) && strcmp(model.ne.method,'state');

% get correct Lambda depending on derivative order
[Lam,Lamgauss] = get_Lam(model.ne,prime,use_core_spline_base);

switch model.ne.method
  case 'state'
    % ne is a state that is solved for
    if ~use_core_spline_base
      if gauss
        dne_dnehat = Lamgauss;
      else
        dne_dnehat = Lam;
      end
      ne = dne_dnehat*x(model.ne.xind,:);
    else
      if size(v,2)>1 
        if any(v(model.hmode.vind.activation, :)) && ~all(v(model.hmode.vind.activation, :))
          error('eval_ne for time dependent V not implemented for varying spline basis')
        end
      end
      xx = [x(model.ne.xind,:); v(model.hmode.vind.ne_rhoedge, :)./model.ne.scal];
      if gauss
        dne_dxx = Lamgauss;
      else
        dne_dxx = Lam;
      end
      dne_dnehat = dne_dxx(:, 1:end-1);
      ne = dne_dxx*xx;
    end
    dne_dx(:,model.ne.xind) = dne_dnehat;
  case 'direct'
    % ne is directly specified by v
    if gauss
      ne = Lamgauss*v(model.ne.vind,:);
    else
      ne = Lam*v(model.ne.vind,:);
    end
  otherwise
    error('model.ne.method = %s not implemented',model.ne.method)
end


return
