function [yp,dyp_dx] = eval_yp(x,v,eqmodel,myrgrid,gaussgrid)
% function [yp,dyp_dx] = eval_yp(x,v,model,whichgrid)
% evalate dy/drho for quantity that does not depend on g
% #codegen

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

switch eqmodel.method
    case 'state'
    % y is a state that is solved for
        if gauss
            dyp_dyhat = eqmodel.Lampgauss;
        else
            dyp_dyhat = eqmodel.Lamp;
        end
        yp = dyp_dyhat*x(eqmodel.xind,:);
        dyp_dx = zeros(ngrid,numel(x));
        dyp_dx(:,eqmodel.xind) = dyp_dyhat;
    case 'direct'
        % y is directly specified by v
        if gauss
            yp = eqmodel.Lampgauss*v(eqmodel.vind,:);
        else
            yp = eqmodel.Lamp*v(eqmodel.vind,:);
        end
        dyp_dx = zeros(ngrid,numel(x));
    otherwise
        error('eqmodel.method = %s not implemented',eqmodel.method)
end


return
