function [ze,dze_dx] = eval_ze(x,g,v,model,gaussgrid,varargin)
% function [ze,dze_dx] = eval_ze(x,g,model,whichgrid,prime)
% prime = true -> evaluate spatial derivative instead
% #codegen

prime = false; % default
if numel(varargin)>=1 && ~isempty(varargin{1})
    prime = varargin{1};
    assert(any(prime==[0,1,2]));
end

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

switch model.ze.method
    case 'direct'
        % ze is directly specified by v
        [Lam,Lamgauss] = get_Lam(model.ze,prime);
        if gauss
            ze = Lamgauss*v(model.ze.vind,:);
        else
            ze = Lam*v(model.ze.vind,:);
        end
        dze_dx = zeros(ngrid,numel(x));
    case 'qnze'
        % Calculate Ze + other unknown parameter from known Zeff(rho)/ne(rho)/ni(rho)/n1(rho)...          
        if nargout == 1
            [ze] = eval_qnze('ze',x,g,v,model,gauss,prime);
        else
            [ze,dze_dx] = eval_qnze('ze',x,g,v,model,gauss,prime);
        end
    otherwise
        error('model.ze.method = %s not implemented',model.ze.method)
end

return
