import numpy as np
import matplotlib.pyplot as plt

alpha = 2.*np.sqrt(3.)
beta = (17./5.)**1.5 - alpha

def f_approx(lam):
    global alpha, beta
    return (1.-(alpha*np.cosh(lam/2.)+beta)**(-2./3.))*np.tanh(lam/2.)

def disratio(m1, m2, m3):
    global alpha, beta
    mu = abs(m1-m3)/m2/2.
    nu = (m1+m3)/m2/2.
    tanh_lower = mu/(1.+nu)
    tanh_upper = mu/(12./17.+nu)
    while(tanh_upper - tanh_lower > 1.e-11):
        cosh_lower = 1./np.sqrt((1.-tanh_lower)*(1.+tanh_lower))
        cosh_upper = 1./np.sqrt((1.-tanh_upper)*(1.+tanh_upper))
        tanh_lower = mu/(1.+nu-(alpha*cosh_upper+beta)**(-2./3.))
        tanh_upper = mu/(1.+nu-(alpha*cosh_lower+beta)**(-2./3.))
    lam = np.arctanh((tanh_lower+tanh_upper)/2.)*2.
    x3 = np.exp(lam)
    x2 = (x3+1.)*f_approx(lam)
    if(m1 < m3):
        r= (1.+x2)/(x3-x2)
    else:
        r= (x3-x2)/(1.+x2)
    #use exact equation to correct
    c = np.empty(6)
    c[5] = m1+m2
    c[4] = (3*m1+2*m2)
    c[3] = 3*m1+m2
    c[2] = -(m2+3*m3)
    c[1] = -(2*m2+3*m3)
    c[0] = -(m2+m3)
    p = np.polynomial.Polynomial(c)
    pp = p.deriv()
    der = pp(r)
    if(abs(der)>0.):
        r -= p(r)/der
        der = pp(r)
        if(abs(der)>0.):
            r -= p(r)/der        
    return r
        
m1 = float(input("enter m_1 = "))
m2 = float(input("enter m_2 = "))
m3 = float(input("enter m_3 = "))
print("dis(2, 3)/dist(1, 2) = ", disratio(m1, m2, m3))