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


m = 1
k = 20
c = 1
l0 = 0.1
g = 9.81
fa = 0.2
fg = 0.1
vt = 0.1

u = lambda t : l0+t*vt

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

fg = 0

def Euler(x0 = 0, y0 = vt, t0 = 0,f = H, dt = 0.1**2, n = 100):
    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)

##LT,LX = Euler(dt = 0.001,n = 3000)

##plt.plot(LT,LX,label = 'sans glsmt')

def glsmt(x,y,t,G):
    T = k*(u(t)-x-l0)
    if G:
        f = 0.1
    else : #Cas adherence :
        f = 0.2
    return(T >= m*g*f)
    
fg = 0.1

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


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

vt = 0.3


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

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

l = 67
T = 24*3600
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 = 360*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)

LT,LX,LY = EDO2()

plt.plot(LX,LY)

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)
plt.plot(np.cos(T),np.sin(T))
plt.legend()
plt.show()