function [nx,dnx_dx] = eval_nx(impindex,x,g,v,model,gauss,prime)
% function [nx,dnx_dx] = eval_nx(impindex,x,g,v,model,whichgrid,prime)
% evaluate value and jacobian for generic impurity (non-main ion
% species)
% prime = true -> returns d/drho of profile instead
% #codegen

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

switch impindex % may not use sprintf due to code generation issues
    case 1
        nxmodel = model.n1;
        nxstr = 'n1';
    case 2
        nxmodel = model.n2;
        nxstr = 'n2';
    case 3
        nxmodel = model.n3;
        nxstr = 'n3';
    otherwise
        error('improper choice of impindex')
end

switch nxmodel.method
    case 'nescal'
        % nx is scaled version of ne.
        % scaling of x coefficients is stored in v
        if gauss
            scal = nxmodel.Lamgauss*v(nxmodel.vind,:);
        else
            scal = nxmodel.Lam*v(nxmodel.vind,:);
        end
        
        if nargout<=1
            if prime==0
                ne = eval_ne(x,g,v,model,gauss);
            elseif prime==1
                ne = eval_nep(x,g,v,model,gauss);
            elseif prime==2
                ne = eval_ne(x,g,v,model,gauss,2);
            end
        else
            if prime==0
                [ne,dne_dx] = eval_ne(x,g,v,model,gauss);
            elseif prime==1
                [ne,dne_dx] = eval_nep(x,g,v,model,gauss);
            elseif prime==2
                [ne,dne_dx] = eval_nepp(x,g,v,model,gauss);
            end
            dnx_dx = bsxfun(@times,scal,dne_dx);
        end
        nx =  scal.*ne;

    case 'niscal'
        % nx is scaled version of ni.
        % scaling of x coefficients is stored in v
        if gauss
            scal = nxmodel.Lamgauss*v(nxmodel.vind,:);
        else
            scal = nxmodel.Lam*v(nxmodel.vind,:);
        end

        if nargout<=1
            if prime==0
                ni = eval_ni(x,g,v,model,gauss);
            elseif prime==1
              ni = eval_nip(x,g,v,model,gauss);
            elseif prime==2
              ni = eval_nipp(x,g,v,model,gauss);
            end
        else
            if prime==0
                [ni,dni_dx] = eval_ni(x,g,v,model,gauss);
            elseif prime==1
                [ni,dni_dx] = eval_nip(x,g,v,model,gauss);       
            elseif prime==2
                [ni,dni_dx] = eval_nipp(x,g,v,model,gauss);       
            end
            dnx_dx = bsxfun(@times,scal,dni_dx);
        end
        nx =  scal.*ni;

    case 'qnze'
        if nargout<=1
           nx = eval_qnze(nxstr,x,g,v,model,gauss,prime);
        else
           [nx,dnx_dx] = eval_qnze(nxstr,x,g,v,model,gauss,prime);
        end
    case 'direct'
        % density given directly via some v
        % ni is directly specified by v, use same matrices as ni.
        [Lam,Lamgauss] = get_Lam(nxmodel,prime);
        if gauss
            nx = Lamgauss*v(nxmodel.vind,:);
        else
            nx = Lam*v(nxmodel.vind,:);
        end
        dnx_dx = zeros(ngrid,numel(x));
    otherwise
        error('not implemented yet')
end

return

