#
#   resolucion de una EDO de segundo r=orden
#   resorte masa amortiguador
#   mediante Runge Kutta de cuarto orden
#
# https://mat.caminos.upm.es/wiki/Sistema_resorte-masa#M.C3.A9todo_de_Runge-Kutta_2
#
# el programa cuenta con dos ventanas una principal interactiva y una secundaria donde s
# se muestra la evolucion de la posicion respecto del tiempo

import tkinter as Tkinter
import time

#parametros ventana principal
w_max=400.0     # ancho de ventana
h_max=300.0     # alto de ventana


# parametros ventana grafico
w_max_wn=250     # ancho de ventana
h_max_wn=160     # alto de ventana
pos_x_wn=1000    # posicion horizontal
pos_y_wn=100     # posicion vertical

tp=[]
xp=[]
yp=[]

def f(x,y):
    return y

def g(x,y):
    #global m, k_r,c
    return -(c*y+k_r*x)/m
            
def dibuja(x,y,z):
    # esta rutina borra el canvas
    tama=10
    spring=20
    canvas.delete(Tkinter.ALL)
    wn_canvas.delete(Tkinter.ALL)
    eje=canvas.create_line(w_max/2-50.,h_max/2,w_max/2+50,h_max/2)
    acum_y=0.
    for i in range(len(x)-1): #iterativamente trazo lineas de un punto al siguiente
            
        #dibujo la masa
        bola=canvas.create_oval(w_max/2-tama,h_max/2+y[i]-tama,w_max/2+tama,h_max/2+y[i]+tama, fill="yellow", width=2)
        
        #dibujo resorte

        resorte1=canvas.create_line(w_max/2,h_max/2,w_max/2+spring,h_max/2+y[i]/4,fill="red")
        resorte2=canvas.create_line(w_max/2-spring,h_max/2+y[i]/4*3,w_max/2,h_max/2+y[i],fill="blue")
        resorte3=canvas.create_line(w_max/2+spring,h_max/2+y[i]/4,w_max/2-spring,h_max/2+y[i]/4*3,fill="green")

         
        canvas.update()
        
        # dibujo el grafico de posicion y velocidad
        if i!= len(x)-2:
            canvas.delete(resorte1)
            canvas.delete(resorte2)
            canvas.delete(resorte3)
            canvas.delete(bola)
        # posicion
        wn_canvas.create_line(x[i],h_max_wn/2+y[i],x[i+1],h_max_wn/2+y[i+1],fill="red")
        # velocidad
        wn_canvas.create_line(x[i],h_max_wn/2+z[i],x[i+1],h_max_wn/2+z[i+1],fill="blue")
    return


def boton():
    # rutina recalcula la solucion d ela ecuacion
    # y llama a rutina de dibujo
    calcula(t0,tN,h,x,v)
    dibuja(tp,xp,yp)
    return


def calcula(t0,tN,h,x0,y0):
    # calcula la respuesta del sistema llenando los vectores tp y xp.
    # para ello utiliza el metodo de Runge Kuta para resolver la ecuacion diferencial de segundo orden 
    for i in range(divisiones):
        tp.append(float(0.0))
        xp.append(float(0.0))
        yp.append(float(0.0))
    
 
    #primeros valores del vector solucion
    x=x0
    y=y0
    n=0
    t=t0
    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*3
        xp[n]=x*50     # amplifico la salida para poder graficar
        yp[n]=y*50
        n=n+1
        t=t+h
    return tp,xp,yp


def regla_masa(m_):
    # rutina que toma el valor del deslizante y lo trasnforma en un numero representano la masa
    global m
    m=float(m_)
    return 
    
def regla_resorte(k_):
    # rutina que toma el valor del deslizante y lo trasnforma en un numero representano la cte del resorte
    global k_r
    k_r=float(k_)
    return
    
def regla_amort(c_):
    # rutina que toma el valor del deslizante y lo trasnforma en un numero representano la cte de amortiguacion
    global c
    c=float(c_)
    return

def regla_posi(p_):
    # rutina que toma el valor del deslizante y lo trasnforma en un numero representano la cte de amortiguacion
    global x
    x=float(p_)
    return


######################
# cuerpo principal
#########################


#tiempo
t0=0.0          # tiempo inicial de simulacion
tN=100.0         # tiempo final de la simulacion
h=0.01          # paso de tiempo
divisiones=int((tN-t0)/h)  
#datos
k_r=4.0        # cte resorte 
m=2.0           # masa del cuerpo
c=1.0           # cte de amortiguacion
 
#condiciones iniciales
x=0.5          #posicion inicial a t=0
v=0.0          #velocidad incial a t=0


#   Invoco el creador de Tk, que me devuelve una ventana a la que nombro root
root = Tkinter.Tk()

#   a ese objeto le modifico la propiedad titulo
root.title("RMA")
#   Invoco el creador del objeto canvas dentro del objeto root y lo llamo canvas
canvas = Tkinter.Canvas(root,width=w_max, height=h_max,bg='white')
#   Hubico el  objeto canvas dentro del onbejto ventana
#   ocupando las dos primeras filas y columnas
canvas.grid(column=0, row=0,columnspan=2,rowspan=3) 

# creo control boton para mandar a calcular
b = Tkinter.Button(root, text ="Lanzar Simulacion", command = boton, activebackground="yellow",activeforeground="Black",width=25)
b.grid(row=3, column=0)

lab1=Tkinter.Label(root,text="PRUEBA")
lab1.grid(row=3,column=1,sticky=Tkinter.E)

#   creo dentro del objeto root una etiqueta y lapongo en la fila 0 columna 4
lab1=Tkinter.Label(root,text="MASA")
lab1.grid(row=0,column=9,sticky=Tkinter.W)
#   creo dentro del objeto root una barra deslizante
#   esta barra cuando es movida, llama a la rutina definida en command
#   las propiedade de la barra de los valores por los que se desliza y su resolucion se definen aqui.
#   se hubica la barra en forma horizontal en la fila 0 columna 5
root_x=Tkinter.Scale(root,orient="vertical",command=regla_masa,from_=1.0,to=100.0,resolution=1.0)
root_x.grid(row=0,column=10)
root_x.set(m)


#   idem para la barra de resorte
lab2=Tkinter.Label(root,text="RESORTE")
lab2.grid(row=1,column=4,sticky=Tkinter.E)
root_y=Tkinter.Scale(root,orient="horizontal",command=regla_resorte,from_=1.0,to=30.0,resolution=1.0)
root_y.grid(row=1,column=5)
root_y.set(k_r)

#   idem para la barra de amortiguacion
lab3=Tkinter.Label(root,text="AMORT")
lab3.grid(row=2,column=4,sticky=Tkinter.E)
root_g=Tkinter.Scale(root,orient="horizontal",command=regla_amort,from_=1.0,to=30.0,resolution=1)
root_g.grid(row=2,column=5)
root_g.set(c)

#   idem para la barra de posicion inicial
lab4=Tkinter.Label(root,text="POSI")
lab4.grid(row=3,column=4,sticky=Tkinter.E)
root_p=Tkinter.Scale(root,orient="horizontal",command=regla_posi,from_=0.0,to=1.0,resolution=0.1)
root_p.grid(row=3,column=5)
root_p.set(x)


# creo otra ventana para el grafico
wn=Tkinter.Toplevel()
wn.title("Grafico")
text_wn=str(w_max_wn)+"x"+str(h_max_wn)+"+"+str(pos_x_wn)+"+"+str(pos_y_wn)
wn.geometry(text_wn)
wn_canvas = Tkinter.Canvas(wn,width=w_max_wn, height=h_max_wn,bg='white')
wn_canvas.grid(row=0,column=0,columnspan=1,rowspan=4) 


root.mainloop()