#   RMA_nuevo.py
#   resolucion de una EDO de segundo r=orden
#   resorte masa amortiguador
#   mediante Runge Kutta de cuarto orden
# http://www.sc.ehu.es/sbweb/fisica3/numerico/diferencial/diferencial.html
# para mathplot  https://matplotlib.org/3.3.2/gallery/index.html

import numpy as np
import matplotlib.pyplot as plt

def f(x,y):
    return y

def g(x,y):
    return -(c*y+k_r*x)/m
    
#tiempo
# este no es un tiempo real, sino y¿una variable que indica el paso del tiempo
t0=0.0          # tiempo inicial de simulacion
tN=10.0         # tiempo final de la simulacion
h=0.01          # paso de tiempo
divisiones=int((tN-t0)/h)  
#datos
k_r=20.0        # cte resorte 
m=1.0           # masa del cuerpo
c=2.0           # cte de amortiguacion
 
#condiciones iniciales
x0=-5.0          #posicion inicial a t=0
v0=0.0          #velocidad incial a t=0

tp=np.zeros(divisiones)    # vector con valores del tiempo
xp=np.zeros(divisiones)    # vector con valores d e la posicion
yp=np.zeros(divisiones)    # vector con valores d e la velocidad
 
#primeros valores del vector solucion
x=x0           # x representa laposicion y la estoy igualando a la posicion en t0
y=v0           # y representa la veliciodad y la estoy igualando a la velocidad en t=0
n=0
t=t0           # pongo el tiempo en 0
while n<=divisiones-1:
    k1=h*f(x,y)
    l1=h*g(x,y)
    k2=h*f(x+k1/2,y+l1/2)
    l2=h*g(x+k1/2,y+l1/2)
    k3=h*f(x+k2/2,y+l2/2)
    l3=h*g(x+k2/2,y+l2/2)
    k4=h*f(x+k3,y+l3)
    l4=h*g(x+k3,y+l3)

    x=x+(k1+2*k2+2*k3*k4)/6
    y=y+(l1+2*l2+2*l3*l4)/6

    tp[n]=t
    xp[n]=x*1
    yp[n]=y*1
    n=n+1
    t=t+h

  
ax=plt.subplot(2,2,1)
ax.plot(tp, xp, color='lightblue', linewidth=1,label='lab')
ax.set_ylabel("Posicion")
ax.set_xlabel("Tiempo")


bx=plt.subplot(2,2,2)
bx.plot(tp, yp, color='red', linewidth=1)
bx.set_ylabel("Velocidad")
bx.set_xlabel("Tiempo")


cx=plt.subplot(2,2,3)
cx.plot(xp, yp, color='green', linewidth=1)
cx.set_ylabel("Velocidad")
cx.set_xlabel("Posicion")

plt.show()