function [ni,dni_dx] = eval_ni(x,g,v,model,gauss)
% function [ni,dni_dx] = eval_ni(x,g,v,model,whichgrid)
% #codegen

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

prime = false; 
switch model.ni.method
    case 'nescal'
        % ni is scaled version of ne.
        % scaling of x coefficients is stored in v
        ne = eval_ne(x,g,v,model,gauss);
        if gauss
            scal = model.ni.Lamgauss*v(model.ni.vind,:);
        else
            scal = model.ni.Lam*v(model.ni.vind,:);
        end
        ni =  scal.*ne;

        if nargout>1
            [~,dne_dx] = eval_ne(x,g,v,model,gauss);
            dni_dx = bsxfun(@times,scal,dne_dx);
        end
    case 'direct'
        % ni is directly specified by v
        if gauss
            ni = model.ni.Lamgauss*v(model.ni.vind,:);
        else
            ni = model.ni.Lam*v(model.ni.vind,:);
        end
        dni_dx = zeros(ngrid,numel(x));
    case 'qnze'
        % Calculate ni+ other unknown parameter from known Zeff(rho)/ne(rho)/ni(rho)/...          
        if nargout<=1
           ni = eval_qnze('ni',x,g,v,model,gauss,prime);
        else
           [ni,dni_dx] = eval_qnze('ni',x,g,v,model,gauss,prime);
        end
    otherwise
        error('not implemented yet')
end


return
