
function DEN = beta_den(ru,rl,phi,omega)
    wphi = phi * omega; % Set omega = 1 for non-shear
    a1 = 3 * (rl * (1 - ru) + ru * (1 - rl));
    a2 = 9 * rl * ru - wphi ^ 2 * (1 - rl) * (1 - ru);
    a3 = 18 * rl * ru + a1 * wphi ^ 2;
    DEN = (a1 * omega ^ 2 - a2) * wphi * sin(wphi) + 18 * rl * ru * omega ^ 2 - (18 * rl * ru * omega ^ 2 + a1 * wphi ^ 2) * cos(wphi);
    %beta_den = 9 * ru * rl * ((omega ^ 2) * (2 - 2 * cos(wphi)) - wphi * sin(wphi)) + _
    %           a1 * ((omega ^ 2) * wphi * sin(wphi) - wphi ^ 2 * cos(wphi)) + _
    %           (1 - ru) * (1 - rl) * wphi ^ 3 * sin(wphi)
    
    % Modification - include the numerator as well so that 1/beta goes to zero
    % Edit: the modification didn't work; it becomes too non-linear to use Newton's on
    % DEN = DEN / (a1 * wphi * cos(wphi) + a2 * sin(wphi))
end