import dynode.main as dyn
import numpy
from pylab import *
import math
solver=dyn.CVOde(dyn.OdeBoussole,dyn.OdeAdams,dyn.OdeFunctional)
            
a1=0.2
a2=0.15
solver.set_cst([(2*math.pi*a1)**2,(2*math.pi*a2)**2,0])     
            
reltol=1e-5
abstol=1e-6
y0=[[0.1,0],[0.5,0],[1,0],[2,0]]
data=[]
for k in range(len(y0)):
    solver.init(0.0,y0[k],reltol,[abstol])
    data.append(solver.solve(1,500))
            
figure(0)
for k in range(len(y0)):
    plot(data[k][1],data[k][2],marker=".",linestyle=" ",markersize=2)
xlabel('theta')
ylabel('dtheta')
            
import numpy.fft
solver.init(0.0,[1,0],reltol,[abstol])
tmax=100
data=solver.solve(0.1,tmax)
tfd=numpy.fft.fft(data[1])
n=size(tfd)
f=numpy.zeros(n)
dB=numpy.zeros(n)
for k in range(n):
    f[k] = (k-1)*1.0/tmax
    dB[k] = 20*math.log10(abs(tfd[k]))
figure(1)
plot(f,dB)
xlabel('f')
ylabel('dB')
axis([0,3,-10,60])
            
solver.close()