function rhoinv = find_rhoinv(q,rhogauss,shear,cp)

indql1 = find(q<1);
indq1 = indql1(end);

rhoinvi = rhogauss(indq1) + (rhogauss(indq1+1) - rhogauss(indq1))*...
    (1. - q(indq1))/(q(indq1+1) - q(indq1));
% set minimum rho_in=0.1
rhoinvi = max(cp.rhoinv, rhoinvi);

if ~cp.sawtooth
    rhoinv = rhoinvi;
else
    % Check shear at q=1
    shear_q1 = shear(indq1);
    % rho_inv fixed at the value before the crash
    if shear_q1>cp.s1crit
        rhoinv = rhoinvi;
    else
        rhoinv = cp.rhoinv;
    end
end