from scipy.integrate import romb,trapz, quad

integral = quad

def integrand1(s,N,x):
    if x == 1:
        return 2*N*s * (1 - exp(-2*N*s))
    return (1 - exp(-2*N*s*x) - exp(-(1-x)) - exp(-2*N*s))/(x*(1-x))

def integrand2(s,N,x):
    if x == 0:
        return 0
    return (exp(2*N*s*x) - 1) * (1 - exp(-2*N*s*x))/(x*(1-x))
def T_kimura(s,N,x):
    J1 = 2/(s*(1 - exp(-2*N*s))) * integral(lambda t: integrand1(s,N,t), x, 1)[0]
    u = (1-exp(-2*N*s*x))/(1-exp(-2*N*s))
    J2 = 2/(s*(1-exp(-2*N*s))) * integral(lambda t: integrand2(s,N,t), 0, x)[0]
    return J1 + ((1-u)/u) * J2
T_kimura(0.01, 1e5, 1e-5)