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,1,100)
z0 = 0
T = 80
y0_max = minv_T(T)[0]
integer_y0_max =  math.ceil(y0_max)
print(y0_max, integer_y0_max) 
y0 = 229.05478

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

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

phi, dphi = get_sol(y0)
plt.plot(r,phi)
plt.grid(True)