function Re_theta = Fasihfar_Johnson(lambda,Tu)
    % Fasihfar and Johnson 1992

    p = [1/82301184 1/4286520 -59/2381400 -74/297675 1369/99225 0];

    l_d = zeros(size(Tu));
    theta_delta = zeros(size(Tu));
    for i = 1:length(l_d)
        if lambda(i) == 0
            l_d(i) = 0;
        else
            p(end) = -lambda(i);
            rr = roots(p);
            rr = rr(real(rr)~=0);
            [~,ii] = min(abs(real(rr)));
            l_d(i) = real(rr(ii));
        end
        theta_delta(i) = (37/315)-l_d(i)/945-l_d(i)^2/9072;
    end

    Re_theta = 111*((l_d+12)./(Tu/100)).^(0.5).*theta_delta;

end

