In [2]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline

from astropy import units
from galpy.orbit import Orbit
from FerrersPotential import FerrersPotential as FP
from galpy.potential import PlummerPotential as PP

In [3]:
%%prun -s cumulative
pmw = FP(amp = 1, a = 8*units.kpc, b = 0.35, c = 0.2375, normalize = True, omegab = 10.*units.km/units.s/units.kpc)
ts = np.linspace(0,34,500)
omwa = Orbit(vxvv=[1,0.25,0.7,0,0.0,0]) #[R,vR,vT,z,vz,phi]
omwb = Orbit(vxvv=[1,0.25,0.7,0,0.0,0]) #[R,vR,vT,z,vz,phi]

#omwa.integrate(ts, pmw, method = 'leapfrog')
#omwb.integrate(ts, ps, method = 'leapfrog')


 

In [7]:
matplotlib.rcParams['figure.figsize'] = (7,7)
omwb.plot(d1 = 'x', d2 = 'y', overplot = False, color = 'red')
omwa.plot(d1 = 'x', d2 = 'y', overplot = True, color = '#00AF3F')



In [192]:
ts = np.linspace(0,500,5000)
omwf = Orbit(vxvv=[1,0.25,0.6,0,0.0,0]) #[R,vR,vT,z,vz,phi]

omwf.integrate(ts, pmw, method = 'leapfrog')

In [193]:
omwf.plot(d1 = 'x', d2 = 'y')



In [552]:
omwg.integrate(ts, ps, method = 'leapfrog')

====================== useful functions =========================


In [4]:
def rot(omega, t):
    temp = [[np.cos(t*omega), np.sin(t*omega)], [-np.sin(t*omega), np.cos(t*omega)]]
    return np.array(temp)

def inrotframe(x,y, ts, potential):
    xy = np.zeros([len(x),2])
    xy[:,0] = x
    xy[:,1] = y
    omega = potential.OmegaP()
    xrot, yrot = np.zeros(len(ts)), np.zeros(len(ts))
    for i in range(len(ts)):
        xrot[i],yrot[i] = np.dot(xy[i],rot(omega, ts[i]))
    return xrot, yrot

In [5]:
def allorbits(x,y):
    xo = [x[2*i] for i in range(int(len(x)/2))]
    xo1 = [x[2*i+1] for i in range(int(len(x)/2))]
    yo = [y[2*i] for i in range(int(len(x)/2))]
    yo1 = [y[2*i+1] for i in range(int(len(x)/2))]
    return [xo,yo],[xo1,yo1]

==================== Lyapunov exponent =======================


In [6]:
x = []
y = []
def evolveorbit(icon, ti, tau, pot):
    global x
    global y
    o = Orbit(vxvv=icon) # [R,vR,vT,z,vz,phi]
    tf = ti+tau
    ts = np.linspace(ti,tf,100)
    o.integrate(ts, pot, method = 'leapfrog')
    x.append(o.x(ts[0]))
    y.append(o.y(ts[0]))
    return [o.R(tf),o.vR(tf),o.vT(tf),o.z(tf),o.vz(tf),o.phi(tf)]  #for full 3D orbit (values can be zero though)

In [7]:
def dvector(o,d):
    return np.array(d)-np.array(o)

In [8]:
def alpha(delta):
    size = np.linalg.norm(delta)
    return size, delta/size

In [9]:
def LEs(time, size, initsize):
    return np.log(size/initsize)

In [10]:
def lyapunov(o,tau, potential, Tm):
    global x
    global y
    x,y = [],[]
    time,LE = [],[] # array large enough
    continuing = True
    i = 1
    w = [1.,0.,0.,0.,0.,0.] # initial setting of deviation vector
    initsize = 1e-5 # may be useful to have it as a parameter of the function
    while continuing:
        icdo = list(np.array(o)+initsize*np.array(w)) # initial conditions for evolved orbit
        newo = evolveorbit(o, tau*i, tau, potential) # integration of tested orbit
        newdo = evolveorbit(icdo, tau*i, tau, potential) # integration of deviated orbit
        wj = dvector(o=newo,d=newdo) # deviation vector after integration
        size, veps0 = alpha(wj) # size of the upper, new deviation vector of size eps0
        LE.append(LEs(tau*i, size, initsize))
        time.append(tau*i)
        if i*tau > Tm:
            break
        i += 1
        o = newo 
        w = veps0
    A = np.array([sum(LE[:i]) / time[i] for i in range(len(time))])
    return A, time

In [12]:
%%prun -s cumulative
#pmw = PP()
pmw = FP(amp = 1, a = 8*units.kpc, b = 0.35, c = 0.2375, normalize = True, omegab = 10.*units.km/units.s/units.kpc)
icon = np.array([1.,0.25,0.7,0.,0.,0.])
tau = 0.01
Tm = 100 #680s
les, time = lyapunov(icon, tau, pmw, Tm)


 

In [13]:
matplotlib.rcParams['figure.figsize'] = (10.0, 5.5)
plt.loglog(time, les)
plt.xlabel(r'$log_{10}\,(t)$', fontsize = 15)
plt.ylabel(r'$log_{10}\,(\lambda)$', fontsize = 15)
#plt.xlim(0,3)


Out[13]:
<matplotlib.text.Text at 0x7f1def9397b8>

In [14]:
[xo,yo], [xo1,yo1] = allorbits(x, y)

In [15]:
matplotlib.rcParams['figure.figsize'] = (7,7)
plt.plot(xo,yo)
plt.plot(xo1,yo1)
plt.xlim(-1.2,1.2)
plt.ylim(-1.2,1.2)


Out[15]:
(-1.2, 1.2)

In [16]:
diff = []
for i in range(len(xo)):
    diff.append((xo[i]-xo1[i])**2+(yo[i]-yo1[i])**2)

In [17]:
#plt.semilogy(diff)


Out[17]:
[<matplotlib.lines.Line2D at 0x7f1def420a90>]

In [18]:
pmw = FP(amp = 1, a = 8*units.kpc, b = 0.35, c = 0.2375, normalize = True, omegab = 10.*units.km/units.s/units.kpc)
orb = Orbit(vxvv = [1,0.25,0.6,0,0.0,0])
timespan = np.linspace(0,100,1000)
orb.integrate(timespan, pmw, method = 'leapfrog')

In [19]:
matplotlib.rcParams['figure.figsize'] = (7,7)
orb.plot(d1 = 'x', d2 = 'y')



In [20]:
%%prun -s cumulative
pmw = FP(amp = 1, a = 8*units.kpc, b = 0.35, c = 0.2375, normalize = True, omegab = 10.*units.km/units.s/units.kpc)
icon = np.array([1.,0.25,0.6,0.,0.,0.])
tau = 0.01
Tm = 200
les, time = lyapunov(icon, tau, pmw, Tm)


 

In [41]:
matplotlib.rcParams['figure.figsize'] = (10.0, 5.5)
plt.loglog(time, les)
plt.xlabel(r'$log_{10}\,(t)$', fontsize = 15)
plt.ylabel(r'$log_{10}\,(\lambda)$', fontsize = 15)


Out[41]:
<matplotlib.text.Text at 0x7f1ded7d0cf8>

In [54]:
pmw = FP(amp = 1, a = 8*units.kpc, b = 0.35, c = 0.2375, normalize = True, omegab = 10.*units.km/units.s/units.kpc)
o23 = Orbit(vxvv = [0.3,0.,0.514355155709012,0.,0.3722961706156539,0.])
timespan = np.linspace(0,40,1000)
o23.integrate(timespan, pmw, method = 'leapfrog')

In [55]:
matplotlib.rcParams['figure.figsize'] = (7,7)
o23.plot(d1 = 'x', d2 = 'y')



In [56]:
o23.plot(d1 = 'x', d2 = 'z')