In [1]:
from __future__ import division, print_function
# Standard library
import os
# Third-party
from astropy import log as logger
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 [2]:
potential = gp.LM10Potential()
In [5]:
fanning = np.array([8.312877511, 0.242593717, 16.811943627, -52.429087, -96.697363, -8.156130])
fanning[3:] = (fanning[3:]*u.km/u.s).to(u.kpc/u.Myr).value
fanning
Out[5]:
In [6]:
potential.total_energy(fanning[None,:3], fanning[None,3:])
Out[6]:
In [43]:
cend = np.loadtxt("/Users/adrian/projects/morphology/simulations/pal5/Pal5_triax_vr0006000all.txt")
cen_w0 = np.median(cend[:,1:7], axis=0)
cen_w0[:3] = cen_w0[:3] / 1000.
cen_w0[3:] = (cen_w0[3:]*u.km/u.s).to(u.kpc/u.Myr).value
t,cen_w = potential.integrate_orbit(cen_w0, dt=0.2, nsteps=50000, Integrator=gi.DOPRI853Integrator)
cen = np.squeeze(cen_w)
#cen_ac,cen_an,cen_fe = gd.find_actions(t, cen, N_max=8, units=galactic)
In [44]:
cen_ac,cen_an,cen_fe = gd.find_actions(t[::20], cen[::20], N_max=10, units=galactic)
In [53]:
Ns = [6,8,10,12,14,16,18]
acts,angs,fres = [],[],[]
for N_max in Ns:
if N_max < 10:
every = 20
else:
every = 5
cen_ac,cen_an,cen_fe = sd.find_actions(t[::every], cen[::every], N_max=N_max, units=galactic)
acts.append(cen_ac)
angs.append(cen_an)
fres.append(cen_fe)
In [65]:
acts = np.array(acts)
angs = np.array(angs)
fres = np.array(fres)
plt.figure(figsize=(8,6))
plt.plot(Ns, acts)
plt.xlabel(r"$N_{\rm max}$", fontsize=26)
plt.ylabel("$J$", fontsize=26)
Out[65]:
In [44]:
iso_p = sd.fit_toy_potential(cen, units=galactic)
aa = iso_p.action_angle(cen[:,:3],cen[:,3:])
In [72]:
plt.figure(figsize=(8,6))
plt.plot(t, aa[0][:,0], marker=None, alpha=0.8, label=r"$J_1$", lw=2.)
plt.plot(t, aa[0][:,1], marker=None, alpha=0.8, label=r"$J_2$", lw=2.)
plt.plot(t, aa[0][:,2], marker=None, alpha=0.8, label=r"$J_3$", lw=2.)
plt.legend(loc='lower right')
plt.xlabel("Time [Myr]")
plt.ylabel("Toy actions [kpc$^{2}$ Myr$^{-1}$]")
plt.xlim(0,10000)
Out[72]:
In [3]:
path = "/Users/adrian/projects/morphology/simulations/pal5/"
action_path = "/Users/adrian/projects/morphology/simulations/pal5/actions"
filename = "nbody6-6000.txt"
rdr = io.NBODY6Reader(path)
pal5 = io.tbl_to_w(rdr.read_snapshot(filename, units=galactic))[:,0]
In [4]:
pal5_2 = np.load(os.path.join(action_path, "w0.npy"))
In [5]:
debris = pal5_2[1:]
In [19]:
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.plot(debris[:,0], debris[:,1], linestyle='none', alpha=0.5, c='#74a9cf', ms=8)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_visible(False)
fig.savefig("/Users/adrian/projects/streams/output/talks/pal5_fan.pdf", transparent=True)
Compute the energy difference between particles and satellite
In [45]:
cen_E = potential.total_energy(cen[:,:3].copy(),cen[:,3:].copy())[0]
debris_E = potential.total_energy(debris[:,:3], debris[:,3:])
In [48]:
dE = debris_E - cen_E
plt.hist(dE);
In [56]:
absdE = np.abs(dE)
print(np.sqrt(np.std(absdE)))
In [20]:
freqs = np.load(os.path.join(action_path, "freqs_{}.npy".format(os.path.splitext(filename)[0])))
actions = np.load(os.path.join(action_path, "actions_{}.npy".format(os.path.splitext(filename)[0])))
In [21]:
actions = np.vstack((cen_ac, actions[1:]))
freqs = np.vstack((cen_fe, freqs[1:]))
In [22]:
actions.shape, dE.shape
Out[22]:
In [60]:
fig = sd.three_panel(actions[1:], c=dE, cmap='RdBu', symbol="J")
In [61]:
fig = sd.three_panel(freqs[1:], c=dE, cmap='RdBu', symbol="J")
In [148]:
ix = freqs[1:,1] < (-0.4*freqs[1:,0] + 0.0004)
ix2 = (freqs[1:,1] < (-0.4*freqs[1:,0] + 0.0015)) & (freqs[1:,1] > (-0.4*freqs[1:,0] + 0.0013))
ix3 = freqs[1:,1] > (-0.4*freqs[1:,0] + 0.0022)
ixs = [ix, ix3]
cs = ['r', 'b', 'g']
################################################
# Frequency figure
fig,axes = plt.subplots(1,3,figsize=(15,5))
axes[0].scatter(freqs[1:,0], freqs[1:,1], alpha=0.2, c='#777777')
axes[1].scatter(freqs[1:,0], freqs[1:,2], alpha=0.2, c='#777777')
axes[2].scatter(freqs[1:,1], freqs[1:,2], alpha=0.2, c='#777777')
for c,ix in zip(cs,ixs):
axes[0].scatter(freqs[1:][ix,0], freqs[1:][ix,1], alpha=0.5, c=c)
axes[1].scatter(freqs[1:][ix,0], freqs[1:][ix,2], alpha=0.5, c=c)
axes[2].scatter(freqs[1:][ix,1], freqs[1:][ix,2], alpha=0.5, c=c)
axes[0].set_xlim(0.02,0.0224)
axes[0].set_ylim(-0.01,-0.004)
axes[1].set_xlim(0.02,0.0224)
axes[1].set_ylim(0.015, 0.016)
axes[2].set_xlim(-0.01,-0.004)
axes[2].set_ylim(0.015, 0.016)
axes[0].set_xlabel(r"$\Omega_1$")
axes[0].set_ylabel(r"$\Omega_2$")
axes[1].set_xlabel(r"$\Omega_1$")
axes[1].set_ylabel(r"$\Omega_3$")
axes[2].set_xlabel(r"$\Omega_2$")
axes[2].set_ylabel(r"$\Omega_3$")
# Particle Positions
fig,axes = plt.subplots(1,3,figsize=(15,5))
axes[0].scatter(debris[:,0], debris[:,1], alpha=0.5, c='#777777')
axes[1].scatter(debris[:,0], debris[:,2], alpha=0.5, c='#777777')
axes[2].scatter(debris[:,1], debris[:,2], alpha=0.5, c='#777777')
for c,ix in zip(cs,ixs):
axes[0].scatter(debris[ix,0], debris[ix,1], alpha=0.4, c=c)
axes[1].scatter(debris[ix,0], debris[ix,2], alpha=0.4, c=c)
axes[2].scatter(debris[ix,1], debris[ix,2], alpha=0.4, c=c)
axes[0].set_xlabel(r"$x$ [kpc]")
axes[0].set_ylabel(r"$y$ [kpc]")
axes[1].set_xlabel(r"$x$ [kpc]")
axes[1].set_ylabel(r"$z$ [kpc]")
axes[2].set_xlabel(r"$y$ [kpc]")
axes[2].set_ylabel(r"$z$ [kpc]")
Out[148]:
Now try using NAFF
In [8]:
t,cen_w = potential.integrate_orbit(fanning, dt=0.5, nsteps=100000,
Integrator=gi.DOPRI853Integrator)
E = potential.total_energy(cen_w[:,0,:3],cen_w[:,0,3:])
plt.semilogy(np.abs(E[1:] - E[0]), marker=None)
Out[8]:
In [11]:
fig = gd.plot_orbits(cen_w, marker='.', linestyle='none', alpha=0.1)
In [13]:
In [14]:
nsteps = len(t)
naff = gd.NAFF(t[:nsteps//2+1])
fs = [cen_w[:nsteps//2+1,0,0]+1j*cen_w[:nsteps//2+1,0,0+i] for i in range(3)]
f1,d,f_ix = naff.find_fundamental_frequencies(fs)
fs = [cen_w[nsteps//2:,0,0]+1j*cen_w[nsteps//2:,0,0+i] for i in range(3)]
f2,d,f_ix = naff.find_fundamental_frequencies(fs)
f1, f2
Out[14]:
In [16]:
# frequency diffusion rate: percent change per Gyr
np.abs((np.abs(f1) - np.abs(f2)) / f1 / (t[20001]-t[0])) * 100. * 1000
Out[16]:
In [19]:
(np.abs(f1) - np.abs(f2)) / f1
Out[19]:
In [57]:
fanning_ball = np.random.normal(fanning, [0.02,0.02,0.02,0.0008,0.0008,0.0008], size=(100,6))
In [58]:
t,ball_w = potential.integrate_orbit(fanning_ball, dt=0.2, nsteps=40000,
Integrator=gi.DOPRI853Integrator)
E = potential.total_energy(ball_w[:,0,:3],ball_w[:,0,3:])
plt.semilogy(np.abs(E[1:] - E[0]), marker=None)
Out[58]:
In [59]:
ball_w.shape
Out[59]:
In [70]:
fig = gd.plot_orbits(ball_w[-1], linestyle='none', alpha=1., subplots_kwargs=dict(sharex=True, sharey=True))
fig.axes[0].set_xlim(-30,30)
fig.axes[0].set_ylim(-30,30)
Out[70]:
In [101]:
def in_hull(p, hull):
"""
Test if points in `p` are in `hull`
`p` should be a `NxK` coordinates of `N` points in `K` dimensions
`hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the
coordinates of `M` points in `K`dimensions for which Delaunay triangulation
will be computed
"""
return hull.find_simplex(p)>=0
In [103]:
from scipy.spatial import ConvexHull, Delaunay
# hull = ConvexHull(ball_w[-1,:,:3])
hull = Delaunay(ball_w[-1,:,:3])
In [99]:
testpts = np.random.uniform(hull.min_bound, hull.max_bound, size=(1000000,3))
In [102]:
init = [in_hull(testpt, hull) for testpt in testpts]
In [97]:
sum(init)
Out[97]:
In [ ]: