import math


def pseudo_root_search_mathias(eos, r, temp, a_mix, b_mix,
                               desired_phase,
                               kappa,
                               search_iterations=10):
    """
    Solves the Mathias constraint given by d_press_d_rho - 0.1*r*temp == 0
    for SRK and PR equations of state. This method is technically independent of the EOS
    provided the relevant derivatives exist.

    Input parameters are the eos, gas_constant as r, temp, a_mix and b_mix which are
    from the eos mixture calculation.

    The desired phase is either 'vap' or 'liq' indicated a vapor like phase or a liquid like phase
    to be returned. This setting also determines the interval for the root search

    search iterations is how long newton search continues before terminating, default is 10.
    Search should conclude in no more than 4-5 iterations.
    """
    # 0 is SRK
    # 1 is PR
    # Bounds for the rho for a given eos

    # This is method is independent of the actual EOS,
    # But need rho_lo, rho_hi, rho_mc, temp_mc and dpress_drho and d2press_drho2 for each each of equation
    # Only SRK and PR are implemented for now.

    # Only the rho changes for these equations, so no new mixing terms need to be calculated
    # Should converge very quickly
    SQRT_2 = math.sqrt(2)

    converged = False
    if eos == 0:
        u, w = 1, 0
        rho_lo, rho_hi = -1 / b_mix, 1 / b_mix
        # From original Mathias paper
        rho_mc = 0.25599 / b_mix
        temp_mc = 0.20268 * a_mix / (r * b_mix)
    elif eos == 1:
        u, w = 2, -1
        rho_lo, rho_hi = (1 - SQRT_2) / b_mix, 1 / b_mix
        # From Watson
        rho_mc = 0.25308 / b_mix
        temp_mc = 0.17014 * a_mix / (r * b_mix)
    else:
        return False, -1, -1, -1, -1, -1

    if desired_phase == 'vap':
        rho_interval_lo, rho_interval_hi = rho_lo, kappa * rho_mc
    elif desired_phase == 'liq':
        rho_interval_lo, rho_interval_hi = rho_mc, rho_hi
    else:
        return False, -1, -1, -1, -1, -1

    scaling = 1 / (r * temp)
    # scaling = 1

    if desired_phase == 'liq':
        # initial_estimate - given by Mathias
        rho_test = rho_hi - 0.4 * (rho_hi - rho_mc)
    else:
        rho_test = (rho_interval_lo + rho_interval_hi) * 0.5
        # rho_test = rho_hi - 0.4*(rho_hi - rho_mc)

    for j in range(search_iterations):
        # EOS in terms of rho (which 1/vol)
        # press = r*temp/(-b_mix + 1/rho) - a_mix/(w*b_mix**2 + u*b_mix/rho + rho**(-2))

        # Derivative of the EOS in terms of rho_test
        d_press_d_rho = r * temp / (rho_test ** 2 * (-b_mix + 1 / rho_test) ** 2) - (
            u * b_mix / rho_test ** 2 + 2 / rho_test ** 3) * a_mix / (
                            w * b_mix ** 2 + u * b_mix / rho_test + rho_test ** (-2)) ** 2
        f = (d_press_d_rho - 0.1 * r * temp)
        f *= scaling
        if f < 1e-6:
            converged = True
            break

        # 2nd Derivative of the EOS in terms of rho root
        d2_press_d_rho_2 = 2 * (
            -r * temp / (b_mix - 1 / rho_test) ** 2 - r * temp / (rho_test * (b_mix - 1 / rho_test) ** 3) + (
            u * b_mix + 3 / rho_test) * a_mix / (
                w * b_mix ** 2 + u * b_mix / rho_test + rho_test ** (-2)) ** 2 - (
                u * b_mix + 2 / rho_test) ** 2 * a_mix / (
                rho_test * (w * b_mix ** 2 + u * b_mix / rho_test + rho_test ** (-2)) ** 3)) / rho_test ** 3
        d2_press_d_rho_2 *= scaling

        df_drho = d2_press_d_rho_2
        rho_test_new = rho_test - f / df_drho
        if rho_test_new < rho_interval_lo:
            rho_test_new = (rho_test + rho_interval_lo) / 2
        elif rho_test_new > rho_interval_hi:
            rho_test_new = (rho_test + rho_interval_hi) / 2

        rho_test = rho_test_new

    # if not converged:
    #     print('press_rho did not converge')

    return converged, rho_mc, rho_lo, rho_hi, rho_test, temp_mc