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]:
array([  8.31287751e+00,   2.42593717e-01,   1.68119436e+01,
        -5.36198651e-02,  -9.88935695e-02,  -8.34137337e-03])

In [6]:
potential.total_energy(fanning[None,:3], fanning[None,3:])


Out[6]:
array([ 0.07198427])

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]:
<matplotlib.text.Text at 0x10d2ff250>


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]:
(0, 10000)

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


0.0119469922702

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]:
((1001, 3), (1000,))

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]:
<matplotlib.text.Text at 0x1212b1450>

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]:
[<matplotlib.lines.Line2D at 0x10aca4f50>]

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]:
(array([-0.01571391, -0.01467847,  0.01571203]),
 array([-0.01570794,  0.01468284,  0.01570922]))

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]:
array([ 0.00379698,  0.00297488,  0.00178777])

In [19]:
(np.abs(f1) - np.abs(f2)) / f1


Out[19]:
array([-0.00037972,  0.0002975 ,  0.00017879])

Integrate a ball


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]:
[<matplotlib.lines.Line2D at 0x12ef6d750>]

In [59]:
ball_w.shape


Out[59]:
(40001, 100, 6)

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]:
(-30, 30)

Convex hull bullshit


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]


ERROR: AttributeError: 'ConvexHull' object has no attribute 'find_simplex' [IPython.core.interactiveshell]
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-102-e17c79afd72c> in <module>()
----> 1 init = [in_hull(testpt, hull) for testpt in testpts]

<ipython-input-101-e45215595376> in in_hull(p, hull)
      8     will be computed
      9     """
---> 10     return hull.find_simplex(p)>=0

AttributeError: 'ConvexHull' object has no attribute 'find_simplex'

In [97]:
sum(init)


Out[97]:
12763

In [ ]: