In [41]:
from __future__ import division, print_function
# Third-party
from IPython.html.widgets import interactive
from IPython.html import widgets as wdg
from IPython.display import display
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# Custom
import gary.dynamics as gd
import gary.integrate as gi
import gary.io as io
import gary.potential as gp
from gary.units import galactic
In [76]:
potential = gp.LeeSutoTriaxialNFWPotential(v_c=0.239225, r_s=30.,
a=1., b=0.75, c=0.5,
units=galactic)
In [95]:
w0 = [19.57343622,0.,31.10792321,0.,0.13777349,0.]
# w0 = [20.,1.,1.,0.,0.26,-0.12]
dt = 2.5
T = 100000.
nsteps = int(T/dt)*2
t,w = potential.integrate_orbit(w0, dt=dt, nsteps=nsteps,
Integrator=gi.DOPRI853Integrator, Integrator_kwargs=dict(nsteps=10))
In [96]:
fig = gd.plot_orbits(w, linestyle='none', alpha=0.1, subplots_kwargs=dict(sharex=True, sharey=True))
In [98]:
E = potential.total_energy(w[:,0,:3].copy(), w[:,0,3:].copy())
dE = np.abs(E[1:] - E[0])
dE.max()
Out[98]:
In [99]:
nintvec = 50
loop = gd.classify_orbit(w)
new_ws = gd.align_circulation_with_z(w, loop)
In [80]:
naff = gd.NAFF(t[:])
fs = gd.poincare_polar(new_ws[:,0])
freqs,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=nintvec)
In [82]:
fprimes = []
for i in range(3):
_d = d[d['n'] == i]
fp = _d['A'][:,np.newaxis] * np.exp(1j * _d['freq'][:,np.newaxis] * t[np.newaxis])
fp = np.sum(fp[:], axis=0)
fprimes.append(fp)
# plotting...
fig,axes = plt.subplots(1,4,figsize=(18,5))
alpha = 0.1
axes[0].plot(fs[i].real, fs[i].imag, linestyle='none', alpha=alpha)
axes[1].plot(fp.real, fp.imag, linestyle='none', alpha=alpha)
axes[2].plot(fs[i].real, marker=None)
axes[2].plot(fp.real, marker=None)
axes[2].set_xlim(0,2000)
axes[3].plot(fs[i].imag, marker=None)
axes[3].plot(fp.imag, marker=None)
axes[3].set_xlim(0,2000)
In [62]:
def dat_interact(ix, n):
_d = d[d['n'] == ix]
fp = _d['A'][:,np.newaxis] * np.exp(1j * _d['freq'][:,np.newaxis] * t[np.newaxis])
fp = np.sum(fp[:n], axis=0)
# plotting...
fig,axes = plt.subplots(1,2,figsize=(10,5))
axes[0].plot(fs[ix].real, marker=None)
axes[0].plot(fp.real, marker=None)
axes[0].set_xlim(0,2000)
axes[1].plot(fs[ix].imag, marker=None)
axes[1].plot(fp.imag, marker=None)
axes[1].set_xlim(0,2000)
In [63]:
v = interactive(dat_interact,
ix=wdg.RadioButtonsWidget(values=[0,1,2],value=0),
n=wdg.FloatSliderWidget(min=1, max=nintvec, step=1, value=1))
display(v)
In [83]:
# dp = double prime
fdp,ddp,ixesdp = naff.find_fundamental_frequencies(fprimes, nintvec=nintvec, break_condition=None)
In [84]:
fdp - freqs
Out[84]:
In [100]:
naff = gd.NAFF(t[:nsteps//2+1])
fs1 = gd.poincare_polar(new_ws[:nsteps//2+1,0])
freqs1,d,ixes = naff.find_fundamental_frequencies(fs1, nintvec=nintvec)
fs2 = gd.poincare_polar(new_ws[nsteps//2:,0])
freqs2,d,ixes = naff.find_fundamental_frequencies(fs2, nintvec=nintvec)
In [101]:
(freqs1 - freqs2) / freqs1
Out[101]:
In [94]:
(freqs1 - freqs2) / freqs1
Out[94]:
In [ ]:
In [ ]:
In [115]:
from scipy.signal import argrelmax,argrelmin
def ptp_freqs(t, *args):
freqs = []
for x in args:
ix = argrelmax(x, mode='wrap')[0]
f_max = np.mean(2*np.pi / (t[ix[1:]] - t[ix[:-1]]))
ix = argrelmin(x, mode='wrap')[0]
f_min = np.mean(2*np.pi / (t[ix[1:]] - t[ix[:-1]]))
freqs.append(np.mean([f_max, f_min]))
return np.array(freqs)
def estimate_max_period(t, w):
if w.ndim < 3:
w = w[:,np.newaxis]
norbits = w.shape[1]
periods = []
for i in range(norbits):
loop = gd.classify_orbit(w[:,i])
if np.any(loop):
# flip coords
new_w = gd.align_circulation_with_z(w[:,i], loop[0])[:,0]
# convert to cylindrical
R = np.sqrt(new_w[:,0]**2 + new_w[:,1]**2)
phi = np.arctan2(new_w[:,1], new_w[:,0])
z = new_w[:,2]
T = 2*np.pi / ptp_freqs(t, R, phi, z)
else:
T = 2*np.pi / ptp_freqs(t, *w[:,i,:3].T)
periods.append(T)
return np.array(periods)
In [119]:
max_T = round(estimate_max_period(t, w).max() * 200, -4)
dt = round(max_T * 1.E-5, 0)
nsteps = int(max_T / dt)
In [ ]:
In [ ]:
In [151]:
# _t,_w = potential.integrate_orbit([0.1,0.,8.62452752,0.,0.30335661,0.], dt=2., nsteps=10000,
_t,_w = potential.integrate_orbit([0.1,0.,0.1,0.,0.38919491,0.], dt=2., nsteps=85000,
Integrator=gi.DOPRI853Integrator,
Integrator_kwargs=dict(nsteps=8192,atol=1E-14,rtol=1E-9))
In [154]:
E = potential.total_energy(_w[:,0,:3],_w[:,0,3:])
np.abs(E[1:] - E[0]).max() / E[0]
Out[154]:
In [153]:
fig = gd.plot_orbits(_w, marker='.', linestyle='none', alpha=0.1)
In [ ]: