function [nip,dnip_dx] = eval_nip(x,g,v,model,gaussgrid)
% function [nip,dnip_dx] = eval_nip(x,g,v,model,whichgrid)
% #codegen

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

prime = 1;

switch model.ni.method
    case 'nescal'
        % nip is scaled version of nep.
        % scaling of x coefficients is stored in v
        nep = eval_nep(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
        nip =  scal.*nep;
 
        if nargout>1
            [~,dnep_dx] = eval_nep(x,g,v,model,gauss);
            dnip_dx = bsxfun(@times,scal,dnep_dx);
        end
    case 'direct'
        % ni is directly specified by v
        if gauss
            nip = model.ni.Lampgauss*v(model.ni.vind);
        else
            nip = model.ni.Lamp*v(model.ni.vind);
        end
        dnip_dx = zeros(ngrid,numel(x));
    case 'qnze'
        % Calculate nip other unknown parameter from known Zeff(rho)/ne(rho)/ni(rho)/...          
        if nargout==1
            [nip] = eval_qnze('ni',x,g,v,model,gauss,prime);
        else
            [nip,dnip_dx] = eval_qnze('ni',x,g,v,model,gauss,prime);
        end
        
    otherwise
        error('model.ni.method = %s not implemented',model.ni.method)
end


return
