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')