##Stick and Slip que j'ai (Sadek) fait il y a quelques années et remodifié cette année (2022)

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.02
vt = 0.3
XY0 = [0,0]
Tmax = 6
dt = 0.001

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


ListeT = []

def HG(x,y,t):
    T = (k/m)*(u(t)-x-l0)
    if y<10**(-3):
        if T > (c/m)*y+g*(fa):
            return (y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*(fa))
        else :
            return (y,0)
    else :
        if y>0:
            return (y,(k/m)*(u(t)-x-l0)-(c/m)*y-g*(fg)*y/abs(y))

    

def Eulerglsmt(x0 = 0, y0 = 0, t0 = 0,f = HG, dt = dt, n = int(Tmax/dt)):
    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()
LT = np.array(LT)

tractant = LX-u(LT)


plt.figure()
plt.plot(LT,LX)
plt.xlabel("Temps")
plt.title("Position de la masse tirée")
#plt.plot(LT[1:],ListeT)
plt.show()

plt.plot(LT, -tractant), plt.xlabel("Temps"),plt.title("Distance Masse tirée - Support tirant"), plt.show()


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