In [2]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
#matplotlib.rcParams['figure.figsize'] = (7.0, 5.5)
from astropy import units
from galpy.orbit import Orbit
from FerrersPotential import FerrersPotential as FP
In [3]:
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 [4]:
def allorbits(x,y):
xo = [x[3*i] for i in range(int(len(x)/3))]
xo1 = [x[3*i+1] for i in range(int(len(x)/3))]
xo2 = [x[3*i+2] for i in range(int(len(x)/3))]
yo = [y[3*i] for i in range(int(len(x)/3))]
yo1 = [y[3*i+1] for i in range(int(len(x)/3))]
yo2 = [y[3*i+2] for i in range(int(len(x)/3))]
return [xo,yo],[xo1,yo1],[xo2,yo2]
In [5]:
# function for normalization of vectors because numpy version is more than 2 times slower!! - but there are other
# complications which make this approach almost just as efficient
def norm(vec):
temp = 0
for item in vec:
temp += item*item
temp = math.sqrt(temp)
for i in range(len(vec)):
vec[i] = vec[i]/temp
return vec
def substr(vec1,vec2):
temp = []
for i in range(len(vec1)):
temp.append(vec1[i]-vec2[i])
return temp
def add(vec1,vec2):
temp = []
for i in range(len(vec1)):
temp.append(vec1[i]+vec2[i])
return temp
def vsize(vec):
temp = 0
for item in vec:
temp += item*item
return math.sqrt(temp)
Set the stopping flag, SF ← 0, the counter, i ← 1, and the orbit characterization variable, OC ← ‘regular’. While ( SF = 0) Do Evolve the orbit and the deviation vectors from time t = (i − 1) τ to t = i τ , i. e. Compute x(i τ ) and w 1 (i τ ), w 2 (i τ ). Normalize the two vectors, i.e. Set w 1 (i τ ) ← w 1 (i τ )/kw 1 (i τ )k and w 2 (i τ ) ← w 2 (i τ )/kw 2 (i τ )k. Compute and Store the current value of the SALI: SALI(i τ ) = min{kw 1 (i τ ) + w 2 (i τ )k, kw 1 (i τ ) − w 2 (i τ )k}. Set the counter i ← i + 1. If [SALI((i − 1) τ ) < S m ] Then Set SF ← 1 and OC ← ‘chaotic’. End If If [(i τ > T M )] Then Set SF ← 1. End If End While Report the time evolution of the SALI and the nature of the orbit.
In [60]:
""" #1
name: devolist
input:
list of initial conditions for the orbit [R,vR,vT,z,vz,phi]; icon
in other than first iteration also 'indev', a list which contains normalized deviation vectors from previous
iteration to be applied to the current state of the Orbit
WARNING: Due to the method of creating the deviation vectors none of the coordinates can be zero!
Otherwise at least two of them are linearly dependent!
output:
list of initial conditions for deviated orbit; devos
"""
def devolist(icon, indev=None):
devos = []
if indev == None:
reldev = 0.001
dorb = icon[:]
for i in range(2):
dorb = icon[:]
dorb[i] += dorb[i]*reldev
devos.append(dorb)
else:
for i in range(2):
dorb = icon[:]
dorb = list(np.array(devo)+np.array(indev[i]))
devos.append(devo)
return devos
In [43]:
""" #2
name: evolveorbit
input:
list of initial conditions for the orbit [R,vR,vT,z,vz,phi]; icon
final time; tf
potential; pot
output:
list of orbit parameters in time tf;
"""
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)]
In [44]:
""" #3
name: dvector
input:
evolved orbit o and deviated orbit d
output:
deviation vector w
"""
def dvector(o,d):
return np.array(d)-np.array(o)
In [96]:
"""
name: sali
input:
initial conditions for the orbit
renormalization time ~ timestep
potential in which the Orbit is inspected
output:
current value of sali, new IC for the Orbit, normalized dev. vectors in matrix for next step
"""
def sali(o, time, tau, potential, w=None):
if type(w) == type(None):
v = devolist(list(o), indev = w)#deviated orbits from the Orbit
else:
v = w
#print(w)
newo = evolveorbit(o, time, tau, potential) #evolves the Orbit
wj, wjnorm = [], []
for dorbit in v:
#print(dorbit+ o*(type(w) != type(None)))
temp = evolveorbit(dorbit + o*(type(w) != type(None)), time, tau, potential) #evolves the deviated orbits to time tau
item = dvector(o=newo,d=temp)
wj.append(item)#creates list of deviation vectors in time tau
wjnorm.append(item/np.linalg.norm(item))
wj = np.array(wj)
wjnorm = np.array(wjnorm)
salivalue = min([np.linalg.norm(wjnorm[0]+wjnorm[1]), np.linalg.norm(wjnorm[0]-wjnorm[1])])
return salivalue, newo, wj #returns value of sali, new IC for the Orbit, x-normalized-x dev. vectors in matrix
In [114]:
"""
name: salievol
input:
initial conditions for the orbit
renormalization time ~ timestep
potential in which the Orbit is inspected
maximum time Tm and small treshold value of the SALI Gm
output:
42 (and list of SALI values during the integration)
"""
def salievol(o,tau, potential, Tm, Gm):
global x
global y
x,y = [],[]
salis = np.zeros([1+int(Tm/tau),2]) # array large enough
continuing = True # stopping flag
i = 1 # counter
w = None #initial setting
while continuing:
salis[i-1][1], newo, a = sali(o, tau*i, tau, potential, w)
salis[i-1][0] = tau*i
if salis[i-1][1]<Gm:
print('if this code works, the orbit is chaotic')
continuing = False
return salis[:i-1][:]
if i*tau > Tm:
print('time reached treshold, orbit is regular if this code works')
continuing = False
return salis[:][:]
i += 1
o = newo
w = a
In [ ]:
#%%prun -s cumulative
Tm=200 #
tau = 1*10**-2
Gm = 10**-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)
icon = np.array([0.5,0.3,0.3,0.,0.,np.pi/50])
results2 = salievol(icon, tau, pmw, Tm, Gm)
In [ ]:
matplotlib.rcParams['figure.figsize'] = (15.0, 5.5)
plt.plot((results1[:,0]), np.log10(results1[:,1]), c= 'blue')
plt.plot((results2[:,0]), np.log10(results2[:,1]), c= 'crimson')
plt.xlabel(r'$t$', fontsize = 15)
plt.ylabel(r'$log_{10}\,SALI_3$', fontsize = 15)
#plt.xlim(-1.5,4.5)
In [ ]:
matplotlib.rcParams['figure.figsize'] = (10, 10)
[xo,yo],[xo1,yo1],[xo2,yo2] = allorbits(x,y)
matplotlib.rcParams['figure.figsize'] = (10.0, 10.0)
plt.plot(xo[:], yo[:], c = 'crimson')
plt.plot(xo1[:], yo1[:], c = 'blue')
plt.plot(xo2[:], yo2[:], c = 'orange')
plt.xlim(-0.6,0.6)
plt.ylim(-0.6,0.6)
plt.show()
In [20]:
matplotlib.rcParams['figure.figsize'] = (15.0, 5.5)
plt.plot(np.log10(results[:,0]), np.log10(results[:,1]), c= 'blue')
plt.xlabel(r'$log_{10}\,t$', fontsize = 15)
plt.ylabel(r'$log_{10}\,SALI_3$', fontsize = 15)
plt.xlim(-1.5,4.5)
Out[20]:
In [21]:
matplotlib.rcParams['figure.figsize'] = (15.0, 5.5)
plt.plot((results[:,0]), np.log10(results[:,1]), c= 'blue')
plt.xlabel(r'$t$', fontsize = 15)
plt.ylabel(r'$log_{10}\,SALI_3$', fontsize = 15)
plt.xlim(0,25000)
plt.ylim(-16,1)
Out[21]:
In [201]:
matplotlib.rcParams['figure.figsize'] = (15.0, 5.5)
plt.plot((results[:,0]), np.log10(results[:,1]), c= 'blue')
plt.xlabel(r'$t$', fontsize = 15)
plt.ylabel(r'$log_{10}\,SALI_3$', fontsize = 15)
plt.xlim(4000,4100)
Out[201]:
In [36]:
matplotlib.rcParams['figure.figsize'] = (10.0, 10.0)
orb = Orbit(vxvv = 0.5,0.3,0.3,0.,0.,np.pi/50])
timespan = np.linspace(0,100,800)
orb.integrate(timespan, pmw, method = 'leapfrog')
orb.plot(d1 = 'x', d2 = 'y')
#plt.xlim(-0.5,0.5)
#plt.ylim(-0.3,0.3)
In [166]:
#%%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,68,1000)
omwa = Orbit(vxvv=[0.4,0.238,0.7,0,0.0,0])
omwa1 = Orbit(vxvv=[1,0.238,0.698,0,0.0,0])
omwa2 = Orbit(vxvv=[1,0.238,0.697,0,0.0,0])
omwa.integrate(ts, pmw, method = 'leapfrog')
#omwa1.integrate(ts, pmw, method = 'leapfrog')
omwa2.integrate(ts, pmw, method = 'leapfrog')
In [167]:
matplotlib.rcParams['figure.figsize'] = (7, 7)
omwa2.plot(d1 = 'x', d2 = 'y',color = 'orange', overplot = False)
omwa.plot(d1 = 'x', d2 = 'y',color = '#00AF3F', overplot = True)
#omwa1.plot(d1 = 'x', d2 = 'y',color = 'crimson')#, overplot = True)
In [168]:
xai,yai = [],[]
for i in ts:
xai.append(omwa.x(i))
yai.append(omwa.y(i))
xar, yar = inrotframe(xai,yai,ts,pmw)
plt.plot(xar,yar)
In [171]:
#%%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,68,1000)
omwb = Orbit(vxvv=[0.2,0.238,0.7,0,0.0,0])
#omwb1 = Orbit(vxvv=[1,0.238,0.698,0,0.0,0])
#omwb2 = Orbit(vxvv=[1,0.238,0.697,0,0.0,0])
omwb.integrate(ts, pmw, method = 'leapfrog')
#omwb1.integrate(ts, pmw, method = 'leapfrog')
#omwb2.integrate(ts, pmw, method = 'leapfrog')
In [173]:
matplotlib.rcParams['figure.figsize'] = (7, 7)
#omwb2.plot(d1 = 'x', d2 = 'y',color = 'orange', overplot = False)
omwb.plot(d1 = 'x', d2 = 'y',color = '#00AF3F')#, overplot = True)
#omwb1.plot(d1 = 'x', d2 = 'y',color = 'crimson')#, overplot = True)
In [177]:
#%%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,68,1000)
omwc = Orbit(vxvv=[0.2,0.238,0.5,0,0.0,0])
#omwc1 = Orbit(vxvv=[1,0.238,0.698,0,0.0,0])
#omwc2 = Orbit(vxvv=[1,0.238,0.697,0,0.0,0])
omwc.integrate(ts, pmw, method = 'leapfrog')
#omwc1.integrate(ts, pmw, method = 'leapfrog')
#omwc2.integrate(ts, pmw, method = 'leapfrog')
In [178]:
matplotlib.rcParams['figure.figsize'] = (7, 7)
#omwc2.plot(d1 = 'x', d2 = 'y',color = 'orange', overplot = False)
omwc.plot(d1 = 'x', d2 = 'y',color = '#00AF3F')#, overplot = True)
#omwc1.plot(d1 = 'x', d2 = 'y',color = 'crimson')#, overplot = True)
In [179]:
xci,yci = [],[]
for i in ts:
xci.append(omwc.x(i))
yci.append(omwc.y(i))
xcr, ycr = inrotframe(xci,yci,ts,pmw)
plt.plot(xcr,ycr)
Out[179]:
In [184]:
omwc.plot(d1 = 'x', d2 = 'vx')
omwc.plot(d1 = 'R', d2 = 'vR')
In [186]:
omwd = Orbit(vxvv = [0.5,0.2,0.4,0.,0.,3*np.pi/10])
ts = np.linspace(0,68,1000)
omwd.integrate(ts, pmw, method ='leapfrog')
In [188]:
omwd.plot(d1 = 'R', d2 = 'vR')
In [199]:
vz = 0.03*16
omwe = Orbit(vxvv = [0.5,0.1,0.6,0.,vz,np.pi/4])
ts = np.linspace(0,500,10000)
omwe.integrate(ts, pmw, method = 'leapfrog')
In [200]:
omwe.plot(d1 = 'R', d2 = 'vR')
In [ ]: