In [1]:

    

from __future__ import division, print_function

import sys
sys.path.append("/Users/adrian/projects/morphology")

# Third-party
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
from streammorphology.initialconditions import loop_grid


    

In [127]:

    

((0.239225/np.sqrt(np.log(2.)-0.5))*u.kpc/u.Myr).to(u.km/u.s)


    

    
Out[127]:





$532.242 \; \mathrm{\frac{km}{s}}$

In [23]:

    

potential = gp.LeeSutoTriaxialNFWPotential(v_c=0.239225, r_s=30.,
                                           a=1., b=0.75, c=0.55,
                                           units=galactic)


    

In [25]:

    

def xz_to_w(E, potential, x, z):
    x = np.atleast_1d(x)
    z = np.atleast_1d(z)
    
    xyz = np.zeros((len(x),3))
    xyz[:,0] = x
    xyz[:,2] = z

    # now, for each grid point, compute the y velocity
    vxyz = np.zeros_like(xyz)
    vxyz[:,1] = np.sqrt(2*(E - potential.value(xyz)))

    return np.hstack((xyz,vxyz))


    

In [28]:

    

w0 = np.load('/Users/adrian/projects/morphology/output/E-0.210_LeeSutoTriaxialNFWPotential_loop/w0.npy')
d = np.memmap('/Users/adrian/projects/morphology/output/E-0.210_LeeSutoTriaxialNFWPotential_loop/allfreqs.dat',
              mode='r', dtype='float64', shape=(len(w0),2,11))


    

In [135]:

    

all_diff = np.abs(d[:,1,3:6] - d[:,0,3:6])
ixx = ((w0[:,0] > 33) & (w0[:,0] < 35.5) & (w0[:,2] > 23.) & (w0[:,2] < 28) & 
       np.all(np.isfinite(all_diff),axis=-1) & 
       np.all(d[:,0,3:6] != 0., axis=-1))

nperiods = d[ixx,0,9]*d[ixx,0,10] / (2*np.pi/np.abs(d[ixx,0,3:6]).max())
frac_diff_rate = np.abs((d[ixx,1,3:6] - d[ixx,0,3:6]) / d[ixx,0,3:6]).max(axis=-1) / nperiods / 2.
sorted_ix = frac_diff_rate.argsort()

region_d = d[ixx][sorted_ix]
frac_diff_rate = frac_diff_rate[sorted_ix]
index = np.where(ixx)[0][sorted_ix]


    

In [136]:

    

plt.plot(w0[ixx,0],w0[ixx,2],linestyle='none')


    

    
Out[136]:





[<matplotlib.lines.Line2D at 0x10c724b10>]






    

In [142]:

    

ixes = [1,int(3*len(frac_diff_rate)/4.),len(frac_diff_rate)-1]

for ix in ixes:
    T = region_d[ix,0,-2]*region_d[ix,0,-1]
    assert region_d[ix,0,8] == 1. # loop 
    
    print(frac_diff_rate[ix])
    print(region_d[ix,0,9], region_d[ix,0,10])
    print()


    

    




8.0008114152e-08
2.0 32500.0

3.97324781559e-05
2.0 32500.0

0.00157051904792
2.0 25000.0

In [138]:

    

for ix in ixes:
    t,w = potential.integrate_orbit(w0[index[ix]], dt=2., nsteps=25000, 
                                    Integrator=gi.DOPRI853Integrator)
    fig = gd.plot_orbits(w, marker=None)


    

In [139]:

    

for ix in ixes:
    ball_w0 = np.random.normal(w0[index[ix]], [0.02,0.02,0.02,0.005,0.005,0.005], size=(1000,6))
    t,w = potential.integrate_orbit(ball_w0, dt=1., nsteps=10000, 
                                    Integrator=gi.DOPRI853Integrator)
    
    fig = gd.plot_orbits(w[-1], marker='.', linestyle='none', subplots_kwargs=dict(sharex=True,sharey=True))
    fig.axes[0].set_xlim(-40,40)
    fig.axes[0].set_ylim(-40,40)


    

In [ ]:

    




    

In [ ]:

    




    

In [149]:

    

for ix,name in zip(ixes,['slow','med','fast']):
    blerg = w0[index[ix]]
    pos = "{},{},{} kpc".format(*blerg[:3])
    vel = "{},{},{} kpc/Myr".format(*blerg[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()


    

    




python scripts/makerun.py --name=slow_freq_diff --pos='34.8135167323,0.0,26.6290009664 kpc' --vel='0.0,0.103136838247,0.0 kpc/Myr' --dt=0.2 --nsteps=50000 --mass=2.5e5 --nsnap=1000 -o

python scripts/makerun.py --name=med_freq_diff --pos='33.5435100226,0.0,26.9735334465 kpc' --vel='0.0,0.108198556941,0.0 kpc/Myr' --dt=0.2 --nsteps=50000 --mass=2.5e5 --nsnap=1000 -o

python scripts/makerun.py --name=fast_freq_diff --pos='34.3901811624,0.0,24.2172736058 kpc' --vel='0.0,0.120832289335,0.0 kpc/Myr' --dt=0.2 --nsteps=50000 --mass=2.5e5 --nsnap=1000 -o

In [ ]: