import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint


m = 10
k = 200
c = 10
l0 = 0.1
g = 9.81
fa = 0.5
fg = 0.2
vt = 0.6
XY0 = [0,0]

u = lambda t : l0+t*vt

def H(x,y,t):
    return (y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*fg)

def H2(XY,t):
    x,y = XY
    return [y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*fg]


def Euler(x0 = 0, y0 = 0, t0 = 0,f = H, dt = 0.1**2, n = 2000):
    x,y,t = x0,y0,t0
    LX,LY,LT =[x],[y],[t]
    for i in range(n):
        dx,dy = f(x,y,t)
        t = t+dt
        x = x + dt*dx
        y = y + dt*dy
        LX.append(x)
        LY.append(y)
        LT.append(t)
    return(LT,LX,LY)

# LT,LX,LY = Euler(dt = 0.01,n = 2000)
# ODEI = odeint(H2,XY0,LT)
# oLX,oLY = ODEI[:,0],ODEI[:,1]


# ##plt.figure('Resolution sans frottement')
# ##
# ##plt.plot(LT,LX,label = 'Euler sf')
# ##plt.plot(LT,oLX,label = 'odint sf')
# ##plt.legend()
# ##
# ##plt.figure('Phase')
# ##
# ##plt.plot(LX,LY,label = 'Euler phase sf')
# ##plt.plot(oLX,oLY,label = 'odeint phase sf')
# ##plt.legend()
# ##plt.show()
# 
# 
# 
# ##plt.show()
# ##plt.show()
#fg = 0.3
def HG(x,y,t):
    T = (k/m)*(u(t)-x-l0)
    if y<10**(-5):
        if T > (c/m)*y+g*(fa):
            print(1)
            return (y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*(fa))
        else :
            print(2)
            return (y,0)
    else :
        if y>0:
            print(3)
            return (y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*(fg))
        # else :
        #     print(4)
        #     return (0,0)

    

def Eulerglsmt(x0 = 0, y0 = 0, t0 = 0,f = HG, dt = 0.01, n = 500):
    x,y,t,G = x0,y0,t0,1
    LX,LY,LT =[x],[y],[t]
    for i in range(n):
        dx,dy = f(x,y,t)
        dx = max(0,dx)
        t = t+dt
        x = x + dt*dx
        y = y + dt*dy
        LX.append(x)
        LY.append(y)
        LT.append(t)
    return(LT,LX,LY)


LT,LX,LY = Eulerglsmt()
plt.figure()
plt.plot(LT,LX, label = 'avec glsmt')
plt.legend()

plt.figure('Phase')
plt.plot(LX,LY)
plt.show()


vt = 0.3


##LT,LX,LY = Eulerglsmt(dt = 0.001,n = 3000)
##
##plt.plot(LT,LX, label = 'avec glsmt')
##plt.legend()
##plt.show()

##################

l = 67
T = 24*360
lamb = 48.86
w = 2*np.pi*np.sin(2*np.pi*lamb/360)/T


def G(x,y,dx,dy,t):
    return(dx,dy,2*w*dy-g*x/l,-2*w*dx-g*y/l)

def EDO2(x0=1,y0=0,dx0 = 0,dy0 = 0,t0 = 0,n = 24, dt = 0.01):
    x,y,dx,dy,t = x0,y0,dx0,dy0,t0
    LX,LY,LT =[x],[y],[t]
    for i in range(n):
        px,py,pdx,pdy = G(x,y,dx,dy,t)
        t = t+dt
        x = x + dt*px
        y = y + dt*py
        dx = dx + dt*pdx
        dy = dy + dt*pdy
        LX.append(x)
        LY.append(y)
        LT.append(t)
    return(LT,LX,LY)
m = 10
LT,LX,LY = EDO2(n = 24*m,dt = 1/m)

#plt.plot(LX,LY,label = 'Euler')

def G_odeint(U,t):
    x,y,dx,dy = U[0],U[1],U[2],U[3]
    return(np.array([dx,dy,2*w*dy-g*x/l,-2*w*dx-g*y/l]))

U0 = np.array([1,0,0,0])

LU = odeint(G_odeint,U0,LT)

LX,LY = LU[:,0],LU[:,1]

plt.plot(LX,LY,'--',label = 'Odeint')
T = np.linspace(0,2*np.pi,200)
plt.plot(np.cos(T),np.sin(T))
plt.legend()
plt.show()