import Funciones
from math import sin, cos, sqrt, pi
from tabulate import tabulate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

def resolverProblema():
    #resolverProblemaPenduloLibre(.5, 0, 0, 33, 300, 16, 1, "Pendulo libre")
    #resolverProblemaPenduloAmortiguado(0, 12, 0, 20, 100, 10, 1, 0.2, 1, "Pendulo Amortiguado")
    #resolverProblemaPenduloForzado(2, 1, 0, 40, 400, 10, 1, 0.2, 1, 1.2, 1, "Pendulo Forzado")

    #ejercicio1()
    #ejercicio2(0.1,1000)
    ejercicio3(1, 5000)
    return

def resolverProblema2Orden(wp, xo, wo, a, b, N):
    xo = (xo + pi) % (2 * pi) - pi
    (t, x, w) = Funciones.RK4Vectorial(wp, float(xo), float(wo), N, float(a), float(b))
    tabla = []
    for i in range(0,len(t)):
        tabla.append((t[i],x[i],w[i]))
    print (tabulate(tabla))
    return (t,x,w)

def resolverProblemaPenduloLibre(xo, wo, a, b, N, g, l, nombre):
    def wp(t,x,w):
        return - float(g)/float(l) * sin(x)
    xo = (xo + pi) % (2 * pi) - pi
    (t, x, w) = Funciones.RK4Vectorial(wp, float(xo), float(wo), N, float(a), float(b))
    tabla = []
    for i in range(0,len(t)):
        tabla.append((t[i],x[i],w[i]))
    print (tabulate(tabla))
    Funciones._plot(nombre, "Tiempo", 'Tita(rad)', t, x, 2.0, "Pos", w, "Vel")
    diagramaDeFase2D("Posicion en funcion del Tiempo",x, w)
    diagramaDeFase3D(t,x,w)
    return (t, x)

def resolverProblemaPenduloAmortiguado(xo, w0, a, b, N, g, l, cteb, m, nombre):
    def wp(t,x,w):
        wo = sqrt(float(g) / float(l))
        gamma = 2 * float(cteb) / float(m)
        return - (wo ** 2) * sin(x) - 2 * gamma * w
    xo = (xo + pi) % (2 * pi) - pi
    (t, x, w) = Funciones.RK4Vectorial(wp, float(xo), float(w0), N, float(a), float(b))
    tabla = []
    for i in range(0,len(t)):
        tabla.append((t[i],x[i],w[i]))
    print (tabulate(tabla))
    Funciones._plot(nombre,"Tiempo", 'Tita(rad)', t, x, 2.0, "Pos", w, "Vel")
    diagramaDeFase2D("Posicion en Funcion del Tiempo)",x, w)
    diagramaDeFase3D(t,x,w)

def resolverProblemaPenduloForzado(xo, w0, a, b, N, g, l, cteb, m, wf, f, nombre):
    def wp(t,x,w):
        wo = sqrt(float(g) / float(l))
        gamma = 2 * float(cteb) / float(m)
        return - (wo ** 2) * sin(x) - 2 * gamma * w + f * cos(wf * t)
    xo = (xo + pi) % (2 * pi) - pi
    (t, x, w) = Funciones.RK4Vectorial(wp, float(xo), float(w0), N, float(a), float(b))
    tabla = []
    for i in range(0,len(t)):
        tabla.append((t[i],x[i],w[i]))
    print (tabulate(tabla))
    Funciones._plot(nombre, "Tiempo", 'Tita(rad)', t, x, 2.0, "Pos", w, "Vel")
    diagramaDeFase2D("Posicion en funcion del Tiempo", x, w)
    diagramaDeFase3D(t,x,w)

def diagramaDeFase2D(titulo, x,w):
    Funciones._plot(titulo, "Tita(rad)", 'Vel(rad/s)', x, w, 2.0)

def diagramaDeFase3D(t,x,w):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(x, w, t, c='r', marker='o')
    plt.ylabel("Vel(rad/s)")
    plt.xlabel("Tita(rad)")
    ax.set_zlabel("Tiempo(s)")
    plt.show()

def ejercicio1():
    def wp1(t, x, w):
        return -16 * sin(x)

    def wp11(t, x, w):
        return -16 * x

    (t , x, w) = resolverProblema2Orden(wp1, 3, 0, 0, 33, 300)
    ej1chico = resolverProblema2Orden(wp1, 0.1, 0, 0, 33, 300)
    ej1sinsen = resolverProblema2Orden(wp11, 3, 0, 0, 33, 300)
    ej1sinsenchico = resolverProblema2Orden(wp11, 0.1, 0, 0, 33, 300)

    Funciones._plot("Tita en Funcion del Tiempo", "Tiempo(s)", 'Tita(rad)', t, x, 2.0, "Pos")
    Funciones._plot2("Comparacion", "Tiempo(s)", 'Tita(rad)', t, x, 2.0, "Ej1", ej1chico[1], "Amplitud baja", otro2=ej1sinsen[1], otroName2="Sin Sen()", otro3=ej1sinsenchico[1], otroName3="A baja y sin Sen()")

    eMax = Funciones.buscarErrorDif(wp1,float(0),float(33),float(3),float(0),2400)
    print (eMax)

    diagramaDeFase2D("Velocidad en Funcion de la Posicion", x, w)
    diagramaDeFase3D(t,x,w)

def ejercicio2(aSobreG, N):
    def wp2(t, x, w):
        return -sin(x) * (1 - aSobreG * sin(t))

    (t,x,w) = resolverProblema2Orden(wp2, 0.1, 0, 0, 200, N)

    emax = Funciones.buscarErrorDif(wp2,float(0),float(200),float(0.1),float(0),N)
    print (emax)

    Funciones._plot("Tita(rad) en Funcion del Tiempo(A/g = " + str(aSobreG) + ")", "Tiempo(s)", 'Tita(rad)', t, x, 2.0, "Pos")
    diagramaDeFase2D("Velocidad en Funcion de la Posicion para A/g = " + str(aSobreG),x,w)
    diagramaDeFase3D(t,x,w)

    return

def ejercicio3(wf, N):
    def wp3(t, x, w):
        return - sin(x) * (1 - wf ** 2 * sin(wf * t))

    (t,x,w) = resolverProblema2Orden(wp3, 0.1, 0, 0, 200, N)

    emax = Funciones.buscarErrorDif(wp3,float(0),float(200),float(0.1),float(0),N)
    print (emax)

    Funciones._plot("Tita(rad) en Funcion del Tiempo (w = " + str(wf) + ")", "Tiempo(s)", "Tita(rad)", t, x, 2.0, "Pos")
    diagramaDeFase2D("Velocidad en Funcion de la Posicion para w = " + str(wf),x,w)
    diagramaDeFase3D(t,x,w)

    return


resolverProblema()
