In [161]:
from __future__ import division, print_function

# Std lib
import glob
import os
import sys

morphology_path = "/Users/adrian/projects/morphology/"
try:
    sys.path.index(morphology_path)
except ValueError:
    sys.path.append(morphology_path)

# Third-party
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
%matplotlib inline
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

# project
from streammorphology.util import potential

In [299]:
# potential = gp.LeeSutoTriaxialNFWPotential(v_c=0.239225, r_s=30.,
#                                            a=1., b=0.75, c=0.55,
#                                            units=galactic)

In [552]:
path = "/Users/adrian/projects/morphology/output/E-0.210_LeeSutoTriaxialNFWPotential_loop"

In [553]:
w0 = np.load(os.path.join(path, 'w0.npy'))
allfreqs_filename = os.path.join(path, 'allfreqs.dat')
d = np.memmap(allfreqs_filename, shape=(len(w0),2,11), dtype='float64', mode='r')

done_ix = d[:,0,7] == 1.
print("{} total orbits".format(len(w0)))
print("{} box orbits".format((d[done_ix,0,8] == 0).sum()))
print("{} loop orbits".format((d[done_ix,0,8] == 1).sum()))


7849 total orbits
811 box orbits
6845 loop orbits

In [554]:
if path.endswith('box'):
    box_ix = done_ix & (d[:,0,8] == 0.) & np.all(np.isfinite(d[:,0,:3]), axis=1)
#     nperiods = d[box_ix,0,9]*d[box_ix,0,10] / (2*np.pi/np.abs(d[box_ix,0,:3]).max(axis=1))
    nperiods = d[:,0,9]*d[:,0,10] / (2*np.pi/np.abs(d[:,0,:3]).max(axis=1))
    
    # max_frac_diff = np.abs((d[box_ix,1,:3] - d[box_ix,0,:3]) / d[box_ix,0,:3]).max(axis=1)
    max_frac_diff = np.abs((d[:,1,:3] - d[:,0,:3]) / d[:,0,:3]).max(axis=1)
    all_freq_diff = np.log10(max_frac_diff / nperiods / 2.)
    box_freq_diff = all_freq_diff[box_ix]
    # np.log10(max_frac_diff[box_ix] / nperiods / 2.)

    box_ix = np.isfinite(box_freq_diff)
    box_freq_diff = box_freq_diff[box_ix]

    # color scaling
    delta = np.abs(box_freq_diff.max() - box_freq_diff.min())
    vmin = box_freq_diff.min() + delta/10.
    vmax = box_freq_diff.max() - delta/10.
else:
    all_freq_diff = np.zeros(len(loop_ix))
    
    loop_ix = done_ix & (d[:,0,8] == 1.) & np.all(np.isfinite(d[:,0,3:6]), axis=1) & np.all(np.isfinite(d[:,1,3:6]), axis=1)
    nperiods = d[loop_ix,0,9]*d[loop_ix,0,10] / (2*np.pi/np.abs(d[loop_ix,0,3:6]).max(axis=1))
    max_frac_diff = np.abs((d[loop_ix,1,3:6] - d[loop_ix,0,3:6]) / d[loop_ix,0,3:6]).max(axis=1)
    loop_freq_diff = np.log10(max_frac_diff / nperiods / 2.)
    all_freq_diff[loop_ix] = loop_freq_diff
    
    box_ix = done_ix & (d[:,0,8] == 0.) & np.all(np.isfinite(d[:,0,:3]), axis=1) & np.all(np.isfinite(d[:,1,:3]), axis=1)
    nperiods = d[box_ix,0,9]*d[box_ix,0,10] / (2*np.pi/np.abs(d[box_ix,0,3:6]).max(axis=1))
    max_frac_diff = np.abs((d[box_ix,1,3:6] - d[box_ix,0,3:6]) / d[box_ix,0,3:6]).max(axis=1)
    all_freq_diff[box_ix] = np.log10(max_frac_diff / nperiods / 2.)

    fin_loop_ix = np.isfinite(loop_freq_diff)
    bad_loop_ix = np.logical_not(fin_loop_ix)
    loop_freq_diff = loop_freq_diff[fin_loop_ix]

    # color scaling
    delta = np.abs(loop_freq_diff.max() - loop_freq_diff.min())
    vmin = loop_freq_diff.min() + delta/10.
    vmax = loop_freq_diff.max() - delta/10.

Checking failed orbits


In [555]:
bad_ix = (d[:,0,7] != 1.) | ((d[:,0,8] == 0.) & np.all(np.isnan(d[:,0,:3]), axis=1)) | ((d[:,0,8] == 1.) & np.all(np.isnan(d[:,0,3:6]), axis=1))

In [556]:
if path.endswith('box'):
    # plot initial condition grid, colored by fractional diffusion rate
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot3D(w0[bad_ix,0], w0[bad_ix,1], w0[bad_ix,2], linestyle='none')
#     ax.plot3D(w0[~bad_ix,0], w0[~bad_ix,1], w0[~bad_ix,2], linestyle='none')
    ax.set_xlim(-1, w0[:,0].max()+5)
    ax.set_ylim(-1, w0[:,1].max()+5)
    ax.set_zlim(-1, w0[:,2].max()+5)
    ax.set_xlabel("$x$")
    ax.set_ylabel("$y$")
    ax.set_zlabel("$z$")

    ax.elev = 45
    ax.azim = 45
else:
    # plot initial condition grid, colored by fractional diffusion rate
    fig,ax = plt.subplots(1,1,figsize=(9,8))
    c = ax.scatter(w0[bad_ix,0], w0[bad_ix,2])
    ax.set_xlim(-1, w0[:,0].max()+5)
    ax.set_ylim(-1, w0[:,2].max()+5)
    ax.set_xlabel(r'$x_0$ $[{\rm kpc}]$')
    ax.set_ylabel(r'$z_0$ $[{\rm kpc}]$')



In [557]:
derp = w0[bad_ix][5]
nsteps = 200000
dt = 2.
t,w = potential.integrate_orbit(derp, dt=dt, nsteps=nsteps, 
                                Integrator=gi.DOPRI853Integrator)
fig = gd.plot_orbits(w, linestyle='none', alpha=0.02)
E = potential.total_energy(w[:,0,:3],w[:,0,3:])
dE = np.abs(E[1:] - E[0])
print(dE.max())


ERROR: KeyboardInterrupt [numpy.core.shape_base]
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-557-079833357b90> in <module>()
      2 nsteps = 200000
      3 t,w = potential.integrate_orbit(derp, dt=2., nsteps=nsteps, 
----> 4                                 Integrator=gi.DOPRI853Integrator)
      5 fig = gd.plot_orbits(w, linestyle='none', alpha=0.02)
      6 E = potential.total_energy(w[:,0,:3],w[:,0,3:])

/Users/adrian/projects/gary/gary/potential/core.pyc in integrate_orbit(self, w0, Integrator, Integrator_kwargs, **time_spec)
    293 
    294         integrator = Integrator(acc, **Integrator_kwargs)
--> 295         return integrator.run(w0, **time_spec)
    296 
    297 class CartesianPotential(Potential):

/Users/adrian/projects/gary/gary/integrate/dopri853.pyc in run(self, w0, mmap, **time_spec)
    154         k = 1
    155         while self._ode.successful() and k < (nsteps+1):
--> 156             self._ode.integrate(times[k])
    157             outy = self._ode.y
    158             ws[k] = outy.reshape(nparticles,ndim)

/Users/adrian/anaconda/lib/python2.7/site-packages/scipy/integrate/_ode.pyc in integrate(self, t, step, relax)
    386             self._y, self.t = mth(self.f, self.jac or (lambda: None),
    387                                 self._y, self.t, t,
--> 388                                 self.f_params, self.jac_params)
    389         except SystemError:
    390             # f2py issue with tuple returns, see ticket 1187.

/Users/adrian/anaconda/lib/python2.7/site-packages/scipy/integrate/_ode.pyc in run(self, f, jac, y0, t0, t1, f_params, jac_params)
    944     def run(self, f, jac, y0, t0, t1, f_params, jac_params):
    945         x, y, iwork, idid = self.runner(*((f, t0, y0, t1) +
--> 946                                           tuple(self.call_args) + (f_params,)))
    947         if idid < 0:
    948             warnings.warn(self.name + ': ' +

/Users/adrian/projects/gary/gary/integrate/dopri853.pyc in func_wrapper(t, x)
    137         def func_wrapper(t,x):
    138             _x = x.reshape((nparticles,ndim))
--> 139             return self.func(t,_x,*self._func_args).reshape((nparticles*ndim,))
    140 
    141         self._ode = ode(func_wrapper, jac=None)

/Users/adrian/projects/gary/gary/potential/core.pyc in <lambda>(t, w)
    290             acc = lambda t,w: self.acceleration(w)
    291         else:
--> 292             acc = lambda t,w: np.hstack((w[...,3:],self.acceleration(w[...,:3])))
    293 
    294         integrator = Integrator(acc, **Integrator_kwargs)

/Users/adrian/anaconda/lib/python2.7/site-packages/numpy/core/shape_base.pyc in hstack(tup)
    270 
    271     """
--> 272     arrs = [atleast_1d(_m) for _m in tup]
    273     # As a special case, dimension 0 of 1-dimensional arrays is "horizontal"
    274     if arrs[0].ndim == 1:

/Users/adrian/anaconda/lib/python2.7/site-packages/numpy/core/shape_base.pyc in atleast_1d(*arys)
     53             result = ary
     54         res.append(result)
---> 55     if len(res) == 1:
     56         return res[0]
     57     else:

KeyboardInterrupt: 

In [ ]:
loop = gd.classify_orbit(w)
ww = gd.align_circulation_with_z(w, loop)

In [ ]:
fig = gd.plot_orbits(ww[:nsteps//2], linestyle='none', alpha=0.02)
fig = gd.plot_orbits(ww[nsteps//2:], linestyle='none', alpha=0.02)

In [124]:
naff = gd.NAFF(t[:nsteps//2+1])

if np.any(loop):
    fs1 = gd.poincare_polar(ww[:nsteps//2+1,0])
    fs2 = gd.poincare_polar(ww[nsteps//2:,0])
else:
    fs1 = [(w[:nsteps//2+1,0,i] + 1j*w[:nsteps//2+1,0,i+3]) for i in range(3)]
    fs2 = [(w[nsteps//2:,0,i] + 1j*w[nsteps//2:,0,i+3]) for i in range(3)]
    
ff1,dd1,ixes = naff.find_fundamental_frequencies(fs1, nintvec=15)
ff2,dd2,ixes = naff.find_fundamental_frequencies(fs2, nintvec=15)

In [129]:
print(ff1)
print(ff2)
print(np.abs((ff2-ff1)/ff1)*100.)


[-0.00493796 -0.00610292 -0.00704913]
[-0.00491576 -0.00609525 -0.00700765]
[ 0.44963651  0.12566835  0.58835493]

In [30]:
d[bad_ix][5][0]


Out[30]:
memmap([ -6.14463013e-03,  -6.13618307e-03,  -6.52241937e-03,
                    nan,              nan,              nan,
         2.00500991e-13,   1.00000000e+00,   1.00000000e+00,
         2.00000000e+00,   5.25000000e+04])


In [558]:
if path.endswith('box'):
    # plot initial condition grid, colored by fractional diffusion rate
    fig = plt.figure(figsize=(14,10))
    ax = fig.add_subplot(111, projection='3d')
#     c = ax.scatter3D(w0[box_ix,0], w0[box_ix,1], w0[box_ix,2], c=box_freq_diff,
#                      vmin=vmin, vmax=vmax, cmap='Greys', s=20, marker='o')
    c = ax.scatter3D(w0[~bad_ix,0], w0[~bad_ix,1], w0[~bad_ix,2], c=all_freq_diff[~bad_ix],
                     vmin=vmin, vmax=vmax, cmap='Greys', s=12, marker='o')
    ax.scatter3D(w0[bad_ix,0], w0[bad_ix,1], w0[bad_ix,2], c='k', s=12, marker='o')
    ax.set_xlim(-1, w0[:,0].max()+5)
    ax.set_ylim(-1, w0[:,1].max()+5)
    ax.set_zlim(-1, w0[:,2].max()+5)
    ax.set_xlabel("$x$")
    ax.set_ylabel("$y$")
    ax.set_zlabel("$z$")

    ax.elev = 45
    ax.azim = 45
    fig.colorbar(c)
else:
    # plot initial condition grid, colored by fractional diffusion rate
    fig,ax = plt.subplots(1,1,figsize=(9,8))
#     c = ax.scatter(w0[loop_ix,0][fin_loop_ix], w0[loop_ix,2][fin_loop_ix], c=loop_freq_diff,
#                    vmin=vmin, vmax=vmax, cmap='jet', s=20, marker='s')
    c = ax.scatter(w0[:,0], w0[:,2], c=all_freq_diff,
                   vmin=vmin, vmax=vmax, cmap='jet', s=20, marker='s')
    ax.set_xlim(-1, w0[:,0].max()+5)
    ax.set_ylim(-1, w0[:,2].max()+5)
    ax.set_xlabel(r'$x_0$ $[{\rm kpc}]$')
    ax.set_ylabel(r'$z_0$ $[{\rm kpc}]$')
    fig.colorbar(c)
    
#     c.set_sizes([500])
#     ax.set_xlim(28,37)
#     ax.set_ylim(11,20)



In [126]:
# plot histograms of diffusion rates
fig,ax = plt.subplots(1,1,figsize=(8,6))
if path.endswith('box'):
    n,bins,pa = ax.hist(box_freq_diff, alpha=0.4, normed=True, label='box')
else:
    n,bins,pa = ax.hist(loop_freq_diff, alpha=0.4, normed=True, label='box')
ax.legend(loc='upper right')
ax.axvline(vmin, alpha=0.1, c='k')
ax.axvline(vmax, alpha=0.1, c='k')


Out[126]:
<matplotlib.lines.Line2D at 0x10c483990>

In [127]:
# plot frequency maps
fig,ax = plt.subplots(1,1,figsize=(8,8))

if path.endswith('box'):
    box_freqs = d[box_ix,:,:3].mean(axis=1)
    ax.plot(box_freqs[:,0]/box_freqs[:,2], box_freqs[:,1]/box_freqs[:,2],
            linestyle='none', marker='.', alpha=0.3)
    ax.set_xlabel(r'$\Omega_x/\Omega_z$')
    ax.set_ylabel(r'$\Omega_y/\Omega_z$')
    ax.set_xlim(0.6,1)
    ax.set_ylim(0.7,1.1)
else:
    loop_freqs = d[loop_ix,:,3:6].mean(axis=1)
    ax.plot(loop_freqs[:,1]/loop_freqs[:,0], loop_freqs[:,2]/loop_freqs[:,0],
            linestyle='none', marker='.', alpha=0.4)
    ax.set_xlabel(r'$\Omega_\phi/\Omega_R$')
    ax.set_ylabel(r'$\Omega_z/\Omega_R$')
    ax.set_xlim(0.45,0.8)
    ax.set_ylim(0.45,0.8)



Pick out two orbits: one with a large frequency difference, one with a small


In [516]:
# _ix = (w0[:,0] > 20.) & (w0[:,0] < 21) & (w0[:,2] > 30) & (w0[:,2] < 32.5)
_ix = (w0[:,0] > 19.9) & (w0[:,0] < 20.5) & (w0[:,2] > 16.5) & (w0[:,2] < 17)

In [517]:
w0[_ix].shape


Out[517]:
(1, 6)

In [518]:
_w0 = w0[_ix]
_freqdiff = all_freq_diff[_ix]

In [519]:
list(_w0[_freqdiff.argmin()])


Out[519]:
[20.039965934677742, 0.0, 16.7718489104955, 0.0, 0.21218860087477895, 0.0]

In [515]:
plt.figure(figsize=(8,8))
plt.scatter(_w0[:,0], _w0[:,2], c=_freqdiff, s=300, vmin=vmin, vmax=vmax)
# plt.xlim(30.5,33.5)
# plt.ylim(14.5,17.5)


Out[515]:
<matplotlib.collections.PathCollection at 0x11031a410>

In [562]:
fast_w0.tolist(),slow_w0.tolist(),


Out[562]:
([20.526306567230858, 0.0, 30.597284592369814, 0.0, 0.11200968228030868, 0.0],
 [20.526306567230858, 0.0, 32.223806437296204, 0.0, 0.09553080349445207, 0.0])

In [349]:
fast_w0 = _w0[7].copy()
print(_freqdiff[7])


-4.19988695266

In [523]:
nsteps = 150000
dt = 2.
t1,w1 = potential.integrate_orbit(fast_w0, dt=dt, nsteps=nsteps, 
                                  Integrator=gi.DOPRI853Integrator)

In [351]:
E = potential.total_energy(w1[:,0,:3],w1[:,0,3:])
dE = np.abs(E[1:] - E[0])
print(dE.max())


1.21569421196e-14

In [352]:
naff = gd.NAFF(t1[:nsteps//2+1])

freqs1 = []
for slc in [slice(None,nsteps//2+1,None),slice(nsteps//2,None,None)]:
    loop = gd.classify_orbit(w1[slc])
    ww1 = gd.align_circulation_with_z(w1[slc], loop)
    fs1 = gd.poincare_polar(ww1[:,0])
    f1,d1,ixes1 = naff.find_fundamental_frequencies(fs1, nintvec=15)
    freqs1.append(f1)
    print(f1)
freqs1 = np.array(freqs1)


[-0.0114901  -0.00754239 -0.00562795]
[-0.01150306 -0.00753852 -0.0056129 ]

In [353]:
(freqs1[1] - freqs1[0]) / freqs1[0]


Out[353]:
array([ 0.00112829, -0.0005139 , -0.00267345])

In [526]:
T1 = (2*np.pi/freqs1)
2*nsteps / T1


Out[526]:
array([[-548.61184982, -360.12273194, -268.71479654],
       [-549.23084042, -359.93766511, -267.99640005]])

In [572]:
2*28000 / T1


Out[572]:
array([[-102.4075453 ,  -67.22290996,  -50.16009535],
       [-102.52309021,  -67.18836415,  -50.02599468]])

In [527]:
fig = gd.plot_orbits(w1, linestyle='none', alpha=0.1)
fig.suptitle("~500 periods", fontsize=20, y=1.05)

fig = gd.plot_orbits(w1[:nsteps//50], linestyle='none', alpha=0.1)
fig.suptitle("~10 periods", fontsize=20, y=1.05)


Out[527]:
<matplotlib.text.Text at 0x113db1ad0>


In [354]:
slow_w0 = _w0[-1].copy()

In [356]:
E = potential.total_energy(w2[:,0,:3],w2[:,0,3:])
dE = np.abs(E[1:] - E[0])
print(dE.max())


8.74300631892e-15

In [357]:
naff = gd.NAFF(t2[:nsteps//2+1])

freqs2 = []
for slc in [slice(None,nsteps//2+1,None),slice(nsteps//2,None,None)]:
    loop = gd.classify_orbit(w2[slc])
    ww2 = gd.align_circulation_with_z(w2[slc], loop)
    fs2 = gd.poincare_polar(ww2[:,0])
    f2,d2,ixes2 = naff.find_fundamental_frequencies(fs2, nintvec=15)
    freqs2.append(f2)
    print(f2)
freqs2 = np.array(freqs2)


[-0.01124301 -0.00759634 -0.00562151]
[-0.01124301 -0.00759633 -0.00562151]

In [576]:
print((freqs1[1] - freqs1[0]) / (dt*nsteps/2.))
print((freqs2[1] - freqs2[0]) / (dt*nsteps/2.))


[ -6.48205451e-11   1.93801532e-11   7.52303043e-11]
[ -5.74738496e-15   1.05386994e-14   3.84619306e-15]

In [528]:
T2 = (2*np.pi/freqs2)
2*nsteps / T2


Out[528]:
array([[-536.81426835, -362.69828635, -268.40713115],
       [-536.81432324, -362.69818572, -268.40709442]])

In [529]:
fig = gd.plot_orbits(w2, linestyle='none', alpha=0.1)
fig.suptitle("~500 periods", fontsize=20, y=1.05)

fig = gd.plot_orbits(w2[:nsteps//50], linestyle='none', alpha=0.1)
fig.suptitle("~10 periods", fontsize=20, y=1.05)


Out[529]:
<matplotlib.text.Text at 0x117ab4950>


In [359]:
print(slow_w0)
print(fast_w0)


[ 20.52630657   0.          32.22380644   0.           0.0955308    0.        ]
[ 20.52630657   0.          30.59728459   0.           0.11200968   0.        ]

In [330]:
# for name,diff_w0 in zip(['slow','fast'],[slow_w0,fast_w0]):
#     pos = "{},{},{} kpc".format(*diff_w0[:3])
#     vel = "{},{},{} kpc/Myr".format(*diff_w0[3:])
    
#     cmd = ("python scripts/makerun.py --name={name}_freq_diff --pos='{pos}' "
#            "--vel='{vel}' --dt=0.2 --nsteps=50000 --mass=2.5e5 --nsnap=1000 -o")
#     print(cmd.format(pos=pos, vel=vel, name=name))
    
#     print()

In [267]:
(potential.parameters['v_c']*u.kpc/u.Myr).to(u.km/u.s).value / np.sqrt(np.log(2.)-0.5)


Out[267]:
532.24177023360915

In [370]:
scf_dt = 0.5
scf_nsteps = 50000

In [371]:
new_w0 = dict()
for name,www0 in zip(['slow','fast'],[slow_w0, fast_w0]):
    t,w = potential.integrate_orbit(www0, dt=1., nsteps=1000, Integrator=gi.DOPRI853Integrator)
    R = np.sqrt(np.sum(w[:,0,:3]**2,axis=-1))

    fig = gd.plot_orbits(w, marker=None)
    fig = gd.plot_orbits(w[0], marker='o', axes=fig.axes)
    fig = gd.plot_orbits(w[R.argmax()], marker='o', axes=fig.axes)

    t,w = potential.integrate_orbit(w[R.argmax()].copy(), dt=-scf_dt, nsteps=scf_nsteps, 
                                    Integrator=gi.DOPRI853Integrator)
    
    new_w0[name] = w[-1,0]
    pos = "{},{},{} kpc".format(*w[-1,0][:3])
    vel = "{},{},{} kpc/Myr".format(*w[-1,0][3:])
    
    cmd = ("python scripts/makerun.py --name={name}_freq_diff_2.5e4 --pos='{pos}' "
           "--vel='{vel}' --dt={dt} --nsteps={nsteps} --mass=2.5e4 --nsnap=1000 -o")
    print(cmd.format(pos=pos, vel=vel, name=name, dt=scf_dt, nsteps=scf_nsteps))
    
    print()


python scripts/makerun.py --name=slow_freq_diff_2.5e4 --pos='-17.504234723,-17.2283745157,-9.07711397761 kpc' --vel='0.0721992194851,0.021293129758,0.199775306493 kpc/Myr' --dt=0.2 --nsteps=50000 --mass=2.5e4 --nsnap=1000 -o

python scripts/makerun.py --name=fast_freq_diff_2.5e4 --pos='-1.69295332221,-13.78418595,15.6309115075 kpc' --vel='0.0580704842,0.228735516722,0.0307028904261 kpc/Myr' --dt=0.2 --nsteps=50000 --mass=2.5e4 --nsnap=1000 -o


In [368]:
t,w = potential.integrate_orbit(new_w0['slow'], dt=scf_dt, nsteps=scf_nsteps, Integrator=gi.DOPRI853Integrator)
fig = gd.plot_orbits(w, marker='.', linestyle='none', alpha=0.1)
fig = gd.plot_orbits(w[-1], marker='o', axes=fig.axes)



In [369]:
t,w = potential.integrate_orbit(new_w0['fast'], dt=scf_dt, nsteps=scf_nsteps, Integrator=gi.DOPRI853Integrator)
fig = gd.plot_orbits(w, marker='.', linestyle='none', alpha=0.1)
fig = gd.plot_orbits(w[-1], marker='o', axes=fig.axes)



In [372]:
slow_w0


Out[372]:
array([ 20.52630657,   0.        ,  32.22380644,   0.        ,
         0.0955308 ,   0.        ])

In [ ]: