import matplotlib.pyplot as plt
import numpy as np
import math

#constants used
H = 2.27e-18
l = 1.5
G = 6.637*(10**(-11))
k = (8*3.14*G)**0.5
om_de = 0.75
omega_matter = 1 - om_de
w0 = -0.85
rho_c = 3*(H**2)/(k**2)
c = ((om_de*(1 + w0)*rho_c)/(H**2))**0.5
v0 = (om_de*(1 - w0)*rho_c)/(2*np.exp(-l*k))

#for iteration
p = 0.0
q = 1.0
n = 10000.0

def f(y1,y2,y3):
    return y1*(omega_matter/(y1**3) + (H**2)*(y3**2)/(2.0*rho_c) + (v0*np.exp(-l*y2*k))/(rho_c))**(1.0/2.0)

def g(y3):
    return y3

def h(y1,y2,y3):
    return -3*y3*(H**2)*(omega_matter/(y1**3) + (H**2)*(y3**2)/(2.0*rho_c) + (v0*np.exp(-l*y2*k))/(rho_c))**0.5 + (l*k*v0*np.exp(-l*y2*k))/(H**2)

#initial conditions
y1 = [1.0]
y2 = [1.0]
y3 = [c]

#the parameters for integration
dh = -(q-p)/n

for j in range(0,int(n)):
    k1 = f(y1[j],y2[j],y3[j])
    r1 = g(y3[j])
    s1 = h(y1[j],y2[j],y3[j])
    k2 = f(y1[j] + k1*dh/2.0,y2[j] + r1*dh/2.0,y3[j] + s1*dh/2.0)
    r2 = g(y3[j] + s1*dh/2.0)
    s2 = h(y1[j] + k1*dh/2.0,y2[j] + r1*dh/2.0,y3[j] + s1*dh/2.0)
    k3 = f(y1[j] + k2*dh/2.0,y2[j] + r2*dh/2.0,y3[j] + s2*dh/2.0)
    r3 = g(y3[j] + s2*dh/2.0)
    s3 = h(y1[j] + k2*dh/2.0,y2[j] + r2*dh/2.0,y3[j] + s2*dh/2.0)
    k4 = f(y1[j] + dh*k3,y2[j] + dh*r3,y3[j] + dh*s3)
    r4 = g(y3[j] + dh*s3)
    s4 = h(y1[j] + dh*k3,y2[j] + dh*r3,y3[j] + dh*s3)
    y1 == y1.append(y1[j] + (dh/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4))
    y2 == y2.append(y2[j] + (dh/6.0)*(r1 + 2.0*r2 + 2.0*r3 + r4))
    y3 == y3.append(y3[j] + (dh/6.0)*(s1 + 2.0*s2 + 2.0*s3 + s4))

def F(Y1,Y2,Y3):
    return Y1*(omega_matter/(Y1**3) + (H**2)*(Y3**2)/(2.0*rho_c) + (v0*np.exp(-l*Y2*k))/(rho_c))**(1.0/2.0)

def G(Y3):
    return Y3

def H(Y1,Y2,Y3):
    return -3*Y3*(H**2)*(omega_matter/(Y1**3) + (H**2)*(Y3**2)/(2.0*rho_c) + (v0*np.exp(-l*Y2*k))/(rho_c))**0.5 + (l*k*v0*np.exp(-l*Y2*k))/(H**2)

dH = (q-p)/n
Y1 = [y1[j]]
Y2 = [y2[j]]
Y3 = [y3[j]]
print(Y1,Y2,Y3)
abc = Y1*Y2
print (abc)

for i in range(0,int(n)):
    K1 = F(Y1[i],Y2[i],Y3[i])
    R1 = G(Y3[i])
    S1 = H(Y1[i],Y2[i],Y3[i])
    K2 = F(Y1[i] + K1*dH/2.0,Y2[i] + R1*dH/2.0,Y3[i] + S1*dH/2.0)
    R2 = G(Y3[i] + Y1*dH/2.0)
    S2 = H(Y1[i] + K1*dH/2.0,Y2[i] + R1*dH/2.0,Y3[i] + S1*dH/2.0)
    K3 = F(Y1[i] + K2*dH/2.0,Y2[i] + R2*dH/2.0,Y3[i] + S2*dH/2.0)
    R3 = G(Y3[i] + S2*dH/2.0)
    S3 = H(Y1[i] + K2*dH/2.0,Y2[i] + R2*dH/2.0,Y3[i] + S2*dH/2.0)
    K4 = F(Y1[i] + dH*K3,Y2[i] + dH*R3,Y3[i] + dH*S3)
    R4 = G(Y3[i] + dH*S3)
    S4 = H(Y1[i] + dH*K3,Y2[i] + dH*R3,Y3[i] + dH*S3)
    Y1 == Y1.append(Y1[i] + (dH/6.0)*(K1 + 2.0*K2 + 2.0*K3 + K4))
    Y2 == Y2.append(Y2[i] + (dH/6.0)*(R1 + 2.0*R2 + 2.0*R3 + R4))
    Y3 == Y3.append(Y3[i] + (dH/6.0)*(S1 + 2.0*S2 + 2.0*S3 + S4))

#other calculations
A = [(H)*((omega_matter/(y1[i]**3) + (H**2)*(y3[i]**2)/(2.0*rho_c) + (v0*np.exp(-l*y2[i]*k))/(rho_c))**(1.0/2.0)) for i in range(len(y3))]
A_1 = [(H)*((omega_matter/(Y1[i]**3) + (H**2)*(Y3[i]**2)/(2.0*rho_c) + (v0*np.exp(-l*Y2[i]*k))/(rho_c))**(1.0/2.0)) for i in range(len(Y3))]

B = [y2[i]/(10**4) for i in range(len(y2))]
B_1 = [Y2[i]/(10**4) for i in range(len(Y2))]

C = [y3[i]/(10**4) for i in range(len(y3))]
C_1 = [Y3[i]/(10**4) for i in range(len(Y3))]

w = [((H**2)*y3[i]**2 - 2.0*v0*np.exp(-l*k*y2[i]))/((H**2)*y3[i]**2 + 2.0*v0*np.exp(-l*k*y2[i])) for i in range(len(y3))]
w_1 = [((H**2)*Y3[i]**2 - 2.0*v0*np.exp(-l*k*Y2[i]))/((H**2)*Y3[i]**2 + 2.0*v0*np.exp(-l*k*Y2[i])) for i in range(len(Y3))]

omde = [(((H**2)*(y3[i]**2))/2.0 + (v0*np.exp(-l*k*y2[i]))/rho_c) for i in range(len(y3))]
omde_1 = [(((H**2)*(Y3[i]**2))/2.0 + (v0*np.exp(-l*k*Y2[i]))/rho_c) for i in range(len(Y3))]

omnr = [(1 - omde[i]) for i in range(len(omde))]
omnr_1 = [(1 - omde_1[i]) for i in range(len(omde_1))]

phi_dot = [y3[i]*H for i in range(len(y3))]
phi_dot_1 = [Y3[i]*H for i in range(len(Y3))]

a = np.linspace(start=1.0,stop=10**(-3),num=10001)
N = np.log(a)


#plt.plot(a,phi_dot,'y',label='phi_dot')
#plt.plot(N,phi_dot_1)
#plt.plot(a,omnr,'-r',label='matter density')
#plt.plot(a,omde,'-b',label='scalar density')
plt.plot(N,w,'-k',label='eos')
plt.plot(N,w_1,'-g',label='eos-f')
#plt.plot(N,B,'-g',label='y2')
#plt.plot(N,C,'-y',label='y3')
#plt.plot(a,A,'-y',label='hubble parameter')
#plt.plot(N,y1,'-g',label='y1')
#plt.plot(N,y2,'-m',label='y2')
#plt.plot(N,y3,'-y',label='y3')
plt.xlabel('N')
plt.title('exponential potential')
plt.legend(loc='best')
plt.show()