

In [1]:

    
import os
import glob
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

import streamteam.io as io
from streamteam.units import galactic
import streamteam.dynamics as sd

Intro

I ran a few simulations over the past few days in spherical, oblate, and triaxial halos with and without a disk and bulge. Below I just summarize two simulations -- one for an oblate halo + disk + bulge and one for a triaxial, NFW halo + bulge + disk where the bulge and disk parameters are the same for the two simulations, $q=0.8$ for the oblate case, and the triaxial halo is sort of prolate with a long axis $q_z/q_y = 1.3$ and short axis $q_x/q_y = 0.8$. In both cases, I use a cluster mass of $2.5\times10^5~\mathrm{M}_\odot$.



In [8]:

    
# Ignore all of this -- just some plotting code...
def panel_plot(x, symbol, lim, relative=True):
    if relative:
        x1,x2,x3 = ((x - x[0])/x[0]).T
    else:
        x1,x2,x3 = x.T
    
    fig,axes = plt.subplots(1,3,figsize=(15,5),sharex=True,sharey=True)
    with mpl.rc_context({'lines.marker': '.', 'lines.linestyle': 'none'}):
        axes[0].plot(x1, x3, alpha=0.25)
        axes[1].plot(x2, x3, alpha=0.25)
        axes[2].plot(x1, x2, alpha=0.25)
    
    if relative:
        axes[0].set_xlabel(r"$({sym}_1 - {sym}_{{1,{{\rm sat}}}})/{sym}_{{1,{{\rm sat}}}}$".format(sym=symbol))
        axes[1].set_xlabel(r"$({sym}_2 - {sym}_{{2,{{\rm sat}}}})/{sym}_{{2,{{\rm sat}}}}$".format(sym=symbol))
        axes[0].set_ylabel(r"$({sym}_3 - {sym}_{{3,{{\rm sat}}}})/{sym}_{{3,{{\rm sat}}}}$".format(sym=symbol))
        axes[2].set_ylabel(r"$({sym}_2 - {sym}_{{2,{{\rm sat}}}})/{sym}_{{2,{{\rm sat}}}}$".format(sym=symbol))
        axes[2].set_xlabel(r"$({sym}_1 - {sym}_{{1,{{\rm sat}}}})/{sym}_{{1,{{\rm sat}}}}$".format(sym=symbol))
    else:
        axes[0].set_xlabel(r"${sym}_1$".format(sym=symbol))
        axes[1].set_xlabel(r"${sym}_2$".format(sym=symbol))
        axes[0].set_ylabel(r"${sym}_3$".format(sym=symbol))
        axes[2].set_ylabel(r"${sym}_2$".format(sym=symbol))
        axes[2].set_xlabel(r"${sym}_1$".format(sym=symbol))
        
    with mpl.rc_context({'lines.marker' : 'o', 'lines.markersize' : 10}):
        axes[0].plot(x1[0], x3[0], alpha=0.75, color='r')
        axes[1].plot(x2[0], x3[0], alpha=0.75, color='r')
        axes[2].plot(x1[0], x2[0], alpha=0.75, color='r')
    
    plt.locator_params(nbins=5)
    
    axes[0].set_xlim(*lim)
    axes[0].set_ylim(*lim)
    fig.tight_layout()
    
    return fig

def snap_plot(path):
    scf = io.APWSCFReader(path)
    w = io.tbl_to_w(scf.read_snap("SNAP015", units=galactic))
    fig = sd.plot_orbits(w, marker='.', alpha=0.1, linestyle='none')
    fig = sd.plot_orbits((w[...,3:]*u.kpc/u.Myr).to(u.km/u.s).value, marker='.', alpha=0.1, linestyle='none')
    fig.axes[0].set_xlabel("$v_x$")
    fig.axes[0].set_ylabel("$v_y$")
    fig.axes[1].set_xlabel("$v_x$")
    fig.axes[1].set_ylabel("$v_z$")
    fig.axes[2].set_xlabel("$v_y$")
    fig.axes[2].set_ylabel("$v_z$")
    
def make_action_plots(path, action_lim=0.25, freq_lim=0.05):
    action_path = os.path.join(path, "actions")   
    actions = np.load(glob.glob(os.path.join(action_path, "actions_*.npy"))[0])
    angles = np.load(glob.glob(os.path.join(action_path, "angles_*.npy"))[0])
    freqs = np.load(glob.glob(os.path.join(action_path, "freqs_*.npy"))[0])
    
    fig_actions = panel_plot(actions, symbol="J", lim=(-action_lim,action_lim))
    fig_angles = panel_plot(angles % (2*np.pi), symbol=r"\theta", lim=(0.,2*np.pi), relative=False)
    fig_freqs = panel_plot(freqs, symbol=r"\Omega", lim=(-freq_lim,freq_lim))

Oblate potential - thin stream



In [9]:

    
oblate_path = "/Users/adrian/projects/morphology/simulations/runs/oblate/"

Here is a plot of the final snapshot of the simulation at $t = 6~\mathrm{Gyr}$. The units of the axis are kpc.



In [10]:

    
snap_plot(oblate_path)

Here are plots of the actions, angles, and frequencies in all three projections, but plotted relative to the progenitor. So, on these axes, 0.01 = 1%, 0.2 = 20% and etc. Red dots mark the position of the progenitor.

The spread in all of these plots matches what I would expect from comparing to Sanders et al. (2013) and Bovy (2014). Especially Figure 4 of Sanders, which shows the equivalent frequency plot for a $10^4~\mathrm{M}_\odot$ cluster (vs. my $2.5\times10^5~\mathrm{M}_\odot$ run). Also see Bovy's Figure 2 to compare all spread in all three quantities.



In [93]:

    
make_action_plots(oblate_path)

Triaxial halo, thin stream



In [89]:

    
triaxial_path = "/Users/adrian/projects/morphology/simulations/runs/thin2/"

Here now is the final snapshot of a thin stream in a triaxial potential, again after $6$ Gyr of evolution.



In [91]:

    
snap_plot(triaxial_path)

Here's where things get funky. Again, here are plots of the actions, angles, and frequencies relative to the progenitor, but notice now I've changed the axes so they go out to 500%! The spread in all quantities is strange, and erratic. Some quantities get as large as $10^8$. Interestingly, the angles seem to bunch up around $\theta = \pi$. Obviously there is some bug...



In [90]:

    
make_plots(triaxial_path, action_lim=5., freq_lim=5.)

Final notes

I compared the actions my code computes vs. Sanders' code for a single random orbit in spherical, oblate, and triaxial potential and they match at the <1% level. But when I compare the values computed for the above orbit, they are wayyy off. I need to debug that and get back to you. For now, hold of on playing with the Hessian until I can figure this out!