Untitled

import math
from scipy import optimize as opt 
import numpy as np
from scipy import integrate
from scipy.optimize import minimize
import matplotlib.pyplot as plt 
from scipy.integrate import quad
from scipy.special import zeta, gamma
mw2 = lambda phi: (m_w**2)*(phi**2/v**2)
mz2 = lambda phi: (m_z**2)*(phi**2/v**2)
mt2 = lambda phi: (m_t**2)*(phi**2/v**2)
cb = 223.6
cf = 14.8
pi = np.pi
m_h = 125
m_w = 80.385
m_z = 91.1876
m_t = 173.07
v = 246.22
M = 600
e = 0.01

# Low temperature approximation to the above integral
JLowT = lambda x: -np.sqrt(2)/4.*(x/pi)**1.5*np.exp(-x)*(1 + 15/(8.*x) + 105/(128.*x*x))
# High temperature
c = pi**2*np.exp(1.5-2*np.euler_gamma)
l = np.arange(20)
Cb = np.sqrt(pi)/2.*(-1)**l*zeta(2*l+3, 1)*gamma(l+1.5)/(gamma(l+4)*(2*pi)**(2*l+5)) 
Cf = (1-2**(2*l+3))*Cb   
def JbHT(x,n=5):
    return -pi**2/90. + x**2/24. - x**3/(12*pi)\
            - x**4*np.log(np.maximum(abs(x**2),1e-100)/(16*c))/(64*pi**2)\
            + np.sum([Cb[i]*x**(2*i+6) for i in range(n)])
def JfHT(x,n=0):
    return -7*pi**2/720. + x**2/48.\
           + x**4*np.log(np.maximum(x**2, 1e-100)/c)/(64*pi**2)\
           + np.sum([Cf[i]*x**(2*i+6) for i in range(n)])
def Jb(x):
    """ Interpolation for bosonic contribution to thermal potential
        (max error ~ 0.054%)
    """
    xlim = 3.72439780262
    if x.real < 0 or x.imag != 0:
        return np.nan
    if 0. <= x and x < xlim: return JbHT(x,n=5)
    elif x > xlim: return JLowT(x)
    elif np.isnan(x): return np.nan

def Jf(x):
    """ Interpolation for bosonic contribution to thermal potential
        (max error ~ 3.89%)
        With this interpolation, the numerical derivative is continuous
        and has a max error ~ 5.5% w.r.t. the full integral
    """
    if x.real < 0 or x.imag != 0:
        return np.nan
    if 0. <= x and x < 1.1: 
        return -(-7*pi**2/720. + x**2/48. + x**4*np.log(np.maximum(x**2, 1e-100)/c)/(64*pi**2))
    elif 1.1 <= x and x <= 3.4: return -(0.0855712017*x - 0.6086856492)/\
                                        (6.321477772 - 0.7249924836*x + x**2)
    elif x > 3.4: 
        return np.sqrt(2)/4.*(x/pi)**1.5*np.exp(-x)*(1 + 15/(8.*x))
    elif np.isnan(x): return np.nan
    
    
def jb(y_t):
  return ((2*pi**2)**-1)*quad((lambda x: x**2 * (np.log(1-np.exp(-np.sqrt(x**2+y_t**2))))), 0, np.infty)[0]
def jf(z_t):
  return ((2*pi**2)**-1)*quad((lambda x: x**2 * (np.log(1+np.exp(-np.sqrt(x**2+z_t**2))))), 0, np.infty)[0]
   

def V(phi,T, M):
    A = 3*v**4/(4*M**2) - m_h**2/2
    B = -3*v**2/(2*M**2) + m_h**2/(2*v**2)
    V0 = A*(phi**2)/2 + (B*phi**4)/4 + (1/(8*M**2))*phi**6  
    MASSA_W = ((mw2(phi)**2)*np.log(mw2(phi)/m_w**2)) -((3*mw2(phi)**2)/2) + (2*mw2(phi))*(m_w**2) - (m_w**4)/2
    MASSA_Z = ((mz2(phi)**2)*np.log(mz2(phi)/m_z**2)) -((3*mz2(phi)**2)/2) + (2*mz2(phi))*(m_z**2) - (m_z**4)/2
    MASSA_t = ((mt2(phi)**2)*np.log(mt2(phi)/m_t**2)) -((3*mt2(phi)**2)/2) + (2*mt2(phi))*(m_t**2) - (m_t**4)/2
    V1 = (1/(64*(pi**2)))*(6*MASSA_W+3*MASSA_Z-12*MASSA_t)
    if T!=0:
        del_V1 = (T**4*((6*Jb(np.sqrt(mw2(phi))/T) + (3*Jb(np.sqrt(mz2(phi))/T)) - (12*Jf(np.sqrt(mt2(phi))/T)))))         
        return V0+V1+del_V1
    else:
        return V0+V1
def v_linha(phi,T):
    return V(phi,T, M)-V(1e-15,T, M)


def find_tc(low, high):
    def vmin(T):
        phimin = minimize(lambda phi: v_linha(phi,T),v).x
        return v_linha(phimin, T)
    return opt.bisect(vmin, low, high)


#plot effective potential graph
#plt.figure(num=0,dpi=120.0)
philist = np.arange(1e-15,300,.1)
#plt.plot(philist, list(map(lambda x: v_linha(x,80), philist)))
#plt.xlabel('$\phi$')
#plt.ylabel('V($\phi$,T)')
#plt.title('Potencial Efetivo')
plt.grid(True)


def minv_T(Tc):
    """ Find the minimum value for phi (x.axis) """
    minimo_v = minimize(lambda phi: v_linha(phi,Tc),v).x
    csi = minimo_v / Tc
    return minimo_v, csi


# Sphaleron with boundary conditions: y0;v0
r = np.linspace(1e-5,2,100)
z0 = 0
T = 80
y0_max = minv_T(T)[0]
integer_y0_max =  math.ceil(y0_max)
print(y0_max, integer_y0_max) 


def deriv1(phi):
    return (V(phi+e,T,M) - V(phi-e,T,M))/(2*e)

def f(r,S):
    y,z = S
    dy_dr = z
    dz_dr = (-2*z)/r + deriv1(y)    
    return [dy_dr, dz_dr]

# Finding lowest valid y
def findY0():
    for i in range(1,integer_y0_max):
        y0 = i
        y, z = get_sol(y0)
        print("Trying y =",(y0), end = "... ")
        if (isCresc(y)):
            prev = y0
        else:
            print("That's it!\nFound integer", (prev), "\ngoing decimal\n")
            break
        print("----------")
        
    y0 = prev
    for dPlace in range (1,7):
        print("New decimal place: currently at", dPlace)
        for j in range (1,10):
            decy0 = y0 + j / 10 ** dPlace
            y, z = get_sol(decy0)
            print("trying dec y = " + (str)(decy0), end = "... ")
            if (isCresc(y)):
                prev = decy0
                print("----------")
            else:
                y0 = prev;
                print ("That's it!\ny must be", y0, "with further decimal places\n")
                break
    print ("Found y0 with", dPlace, "decimal places:", y0)
    return y


# Calc if graph is crescent for y
def isCresc(y):
    for indice in range(10):
        if(y[indice] < y[indice+1]):
            print("Still not it")
            return True
    else:
        return False

# IVP for y
def get_sol(y0):
    sol = integrate.solve_ivp(f, (1e-5,2),(y0,z0),t_eval=r)
    return sol.y

phi = findY0()

plt.plot(r,phi)
plt.grid(True)