#Solucion Par
from matplotlib.pyplot import *
from numpy import *
from scipy.integrate import odeint
from __future__ import division
from math import pi
%matplotlib inline 

# Constantes
B=25*10**(-3)   
L=0.25    
m=0.1    
xi=pi/4
g=9.8
vi=0  


# Calculo de parametros:
wn=sqrt(g/L)       # frecuencia natural
z=B/(2*wn*m)         # constante de amortiguamiento
wd=sqrt(1-z**2)*wn # frecuencia de oscilacion forzada

tf= 9*2*pi/wd # tiempo final de simulacion: 9 periodos
t = linspace(0, tf, 500)

def x(t): # solucion al P.V.I x(0)=xi, x'(0)=vi
    c2=xi;  c1=(z*wn*xi)/wd
    return exp(-z*wn*t)*(c1*sin(wd*t)+c2*cos(wd*t))
def v(t):
    c2=xi  
    c1=(z*wn*xi)/wd
    return exp(-z*wn*t)*((-z*wn)*(c1*sin(wd*t)+c2*cos(wd*t)) + c1*wd*cos(wd*t)-c2*wd*sin(wd*t))


# graficas
fig, (ax0, ax1) = subplots(ncols=2,figsize=(12, 5))
ax0.plot(t,v(t),'k--',label=r'$v(t)$',lw=2)
ax0.set_xlim(0,tf);
ax0.legend(loc='best')
ax1.plot(x(t),v(t),'k',lw=2);
ax1.set_title(u'Diagrama de fase')

fig2=figure()
plot(t,x(t),'k',label=r'$x(t)$',lw=2)
xlim(0,tf);
legend(loc='best')