In [18]:
from __future__ import division, print_function
import sys
# Third-party
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
morphology_path = "/Users/adrian/projects/morphology/"
try:
sys.path.index(morphology_path)
except ValueError:
sys.path.append(morphology_path)
# 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
# project
from streammorphology.potential import potential_registry
potential = potential_registry['triaxial-NFW']
In [45]:
kwargs = dict(marker=None,color='#1b75bb')
In [46]:
# pretzel
t,w = potential.integrate_orbit([31.712141115952516, 0.0, 15.551957526800706, 0.0, 0.15403197509440325, 0.0],
dt=3., nsteps=10000, Integrator=gi.DOPRI853Integrator)
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.set_xlim(-40,40)
ax.set_ylim(-40,40)
ax.plot(w[:,0,0],w[:,0,2],**kwargs)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_facecolor('None')
fig.patch.set_facecolor('None')
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/orbit1.pdf", transparent=True)
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.set_xlim(-40,40)
ax.set_ylim(-40,40)
ax.plot(w[:,0,1],w[:,0,2],**kwargs)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_facecolor('None')
fig.patch.set_facecolor('None')
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/orbit2.pdf", transparent=True)
In [47]:
# closed trefoil
t,w = potential.integrate_orbit([24.417031627655781, 0.0, 0.1, 0.0, 0.23727061350345613, 0.0],
dt=3., nsteps=10000, Integrator=gi.DOPRI853Integrator)
fig = gd.plot_orbits(w, marker=None)
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.set_xlim(-40,40)
ax.set_ylim(-40,40)
ax.plot(w[:,0,0],w[:,0,1],**kwargs)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_facecolor('None')
fig.patch.set_facecolor('None')
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/orbit3.pdf", transparent=True)
In [431]:
# higher order
t,w = potential.integrate_orbit([20.039965934677742, 0.0, 16.7718489104955, 0.0, 0.21218860087477895, 0.0],
dt=3., nsteps=10000, Integrator=gi.DOPRI853Integrator)
fig = gd.plot_orbits(w, marker=None)
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.set_xlim(-40,40)
ax.set_ylim(-40,40)
ax.plot(w[:,0,0],w[:,0,1],**kwargs)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_facecolor('None')
fig.patch.set_facecolor('None')
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/orbit4.pdf", transparent=True)
fig,ax = plt.subplots(1,1,figsize=(8,8))
ax.set_xlim(-40,40)
ax.set_ylim(-40,40)
ax.plot(w[:,0,1],w[:,0,2],**kwargs)
ax.set_frame_on(False)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
ax.patch.set_facecolor('None')
fig.patch.set_facecolor('None')
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/orbit5.pdf", transparent=True)
In [441]:
# fast, slow
w0s = np.array([[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 [442]:
t,ws = potential.integrate_orbit(w0s, dt=2., nsteps=28000, Integrator=gi.DOPRI853Integrator)
In [451]:
loop = gd.classify_orbit(ws)
ws2 = gd.align_circulation_with_z(ws, loop)
In [469]:
fig,axes = plt.subplots(2,2,figsize=(10,10),sharex=True,sharey=True)
axes[0,0].plot(ws2[:,1,0], ws2[:,1,1], marker=None)
axes[1,0].plot(ws2[:,1,0], ws2[:,1,2], marker=None)
axes[0,1].plot(ws2[:,0,0], ws2[:,0,1], marker=None)
axes[1,1].plot(ws2[:,0,0], ws2[:,0,2], marker=None)
axes[0,0].set_xticks(range(-30,31,15))
axes[0,0].set_yticks(range(-30,31,15))
axes[0,0].set_xlim(-40,40)
axes[0,0].set_ylim(-40,40)
fig.tight_layout()
fig.savefig("/Users/adrian/papers/posters/AAS225/figures/figure2.pdf", transparent=True)
In [473]:
dt = 2.
nsteps = 150000
t,ws = potential.integrate_orbit(w0s, dt=dt, nsteps=nsteps,
Integrator=gi.DOPRI853Integrator)
In [478]:
fast_w = ws[:,0:1]
slow_w = ws[:,1:]
In [480]:
naff = gd.NAFF(t[:nsteps//2+1])
fast_freqs = []
for slc in [slice(None,nsteps//2+1,None),slice(nsteps//2,None,None)]:
loop = gd.classify_orbit(fast_w[slc,])
ww = gd.align_circulation_with_z(fast_w[slc], loop)
fs = gd.poincare_polar(ww[:,0])
f,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
fast_freqs.append(f)
fast_freqs = np.array(fast_freqs)
In [481]:
naff = gd.NAFF(t[:nsteps//2+1])
slow_freqs = []
for slc in [slice(None,nsteps//2+1,None),slice(nsteps//2,None,None)]:
loop = gd.classify_orbit(slow_w[slc,])
ww = gd.align_circulation_with_z(slow_w[slc], loop)
fs = gd.poincare_polar(ww[:,0])
f,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
slow_freqs.append(f)
slow_freqs = np.array(slow_freqs)
In [521]:
0.5*50000 / np.abs(2*np.pi/slow_freqs).max()
Out[521]:
In [490]:
(slow_freqs[1] - slow_freqs[0]) / slow_freqs[0] / (dt*nsteps/2.) * 1000.
Out[490]:
In [489]:
(fast_freqs[1] - fast_freqs[0]) / fast_freqs[0] / (dt*nsteps/2.) * 1000.
Out[489]:
In [ ]:
In [ ]:
In [ ]:
In [21]:
slow_w0 = np.array([-17.504234723,-17.2283745157,-9.07711397761,0.0721992194851,0.021293129758,0.199775306493])
fast_w0 = np.array([-1.69295332221,-13.78418595,15.6309115075,0.0580704842,0.228735516722,0.0307028904261])
In [22]:
name = 'slow'
In [49]:
scf = io.APWSCFReader("/Users/adrian/projects/morphology/simulations/runs/fast_freq_diff_2.5e4/")
tbl = scf.read_snap("SNAP050",units=galactic)
tbl.meta
Out[49]:
In [23]:
if name == 'fast':
scf = io.APWSCFReader("/Users/adrian/projects/morphology/simulations/runs/fast_freq_diff_2.5e4/")
w0 = fast_w0
extra = 0*u.deg
fast_cen_w = io.tbl_to_w(scf.read_cen(units=galactic))
fast_w = io.tbl_to_w(scf.read_snap("SNAP050",units=galactic))
elif name == 'slow':
scf = io.APWSCFReader("/Users/adrian/projects/morphology/simulations/runs/slow_freq_diff_2.5e4/")
w0 = slow_w0
# extra = 180*u.deg
extra = 0*u.deg
slow_cen_w = io.tbl_to_w(scf.read_cen(units=galactic))
slow_w = io.tbl_to_w(scf.read_snap("SNAP050",units=galactic))
else:
raise ValueError("Dumby!")
In [24]:
fig = gd.plot_orbits(slow_w[:,3:], marker='.', linestyle='none', alpha=0.2)
fig = gd.plot_orbits(fast_w[:,3:], marker='.', linestyle='none', alpha=0.2)
In [25]:
t,ws = potential.integrate_orbit(w0, dt=0.5, nsteps=50000)
In [26]:
fig = gd.plot_orbits(ws, marker=None)
fig = gd.plot_orbits(cen_w, marker=None, axes=fig.axes)
fig = gd.plot_orbits(cen_w[-1:], marker='s', axes=fig.axes)
fig = gd.plot_orbits(cen_w[0:1], marker='o', axes=fig.axes)
Realign point so that it lies along the x=0 (y) axis
In [27]:
from astropy.coordinates.angles import rotation_matrix
In [28]:
new_cen_x = cen_w[:,:3].copy()
new_cen_v = cen_w[:,3:].copy()
# put endpoint on x axis in x-z plane
theta = np.arctan2(new_cen_x[-1,2],new_cen_x[-1,0]) * u.radian
R = rotation_matrix(-theta, 'y')
new_cen_x = R.dot(new_cen_x.T).T
new_cen_v = R.dot(new_cen_v.T).T
full_R = R
# put endpoint on x axis in x-y plane
theta = np.arctan2(new_cen_x[-1,1], new_cen_x[-1,0]) * u.radian
print(theta.to(u.deg))
R = rotation_matrix(theta, 'z')
new_cen_x = R.dot(new_cen_x.T).T
new_cen_v = R.dot(new_cen_v.T).T
full_R = R*full_R
# now align L
L = np.cross(new_cen_x[-1], new_cen_v[-1])[0]
theta = np.arccos(L[2]/np.sqrt(np.sum(L**2)))*u.radian
R = rotation_matrix(-theta, 'x')
full_R = R*full_R
new_cen_x = np.array(R.dot(new_cen_x.T).T)
new_cen_v = np.array(R.dot(new_cen_v.T).T)
In [29]:
fig = gd.plot_orbits(new_cen_x, marker=None)
fig = gd.plot_orbits(new_cen_x[-1:], marker='s', axes=fig.axes)
fig = gd.plot_orbits(new_cen_x[0:1], marker='o', axes=fig.axes)
In [30]:
# slow_x = np.array(full_R.dot(w[:,:3].T).T)
# slow_cen_x = new_cen_x.copy()
fig = gd.plot_orbits(slow_x, marker='.', linestyle='none', subplots_kwargs=dict(sharex=True, sharey=True), alpha=0.4)
fig.axes[0].set_xlim(-50,50)
fig.axes[0].set_ylim(-50,50)
fig.suptitle('slow', fontsize=22, y=1.02)
In [31]:
# fast_x = np.array(full_R.dot(w[:,:3].T).T)
# fast_cen_x = new_cen_x.copy()
fig = gd.plot_orbits(fast_x, marker='.', linestyle='none', subplots_kwargs=dict(sharex=True, sharey=True), alpha=0.4)
fig.axes[0].set_xlim(-50,50)
fig.axes[0].set_ylim(-50,50)
fig.suptitle('fast', fontsize=22, y=1.02)
In [641]:
fig,axes = plt.subplots(1,2,figsize=(12,5),sharex=True,sharey=True)
for i,(x,cen_x) in enumerate(zip([fast_x[:,:3], slow_x[:,:3]],[fast_cen_x,slow_cen_x])):
r = np.sqrt(np.sum(x[:,:3]**2, axis=-1))
phi = np.arctan2(x[:,1],x[:,0])
axes[i].plot(np.degrees(phi), r, linestyle='none', alpha=0.25)
cen_r = np.sqrt(np.sum(cen_x[:,:3]**2, axis=-1))
cen_phi = np.arctan2(cen_x[:,1], cen_x[:,0])
axes[i].plot(np.degrees(cen_phi[-1]), cen_r[-1], linestyle='none', marker='o')
axes[i].set_xlim(-50,50.)
axes[i].set_ylim(25,40.)
axes[i].set_xlabel("Stream longitude [deg]")
axes[0].set_ylabel(r"Distance [kpc]")
fig.tight_layout()
In [43]:
fig,axes = plt.subplots(1,2,figsize=(12,5),sharex=True,sharey=True)
for i,(x,cen_x) in enumerate(zip([fast_w[:,:3], slow_w[:,:3]],[fast_cen_w[:,:3],slow_cen_w[:,:3]])):
r = np.sqrt(np.sum(x[:,:3]**2, axis=-1))
phi = np.arctan2(x[:,1],x[:,0])
theta = np.arccos(x[:,2]/r)
axes[i].plot(np.degrees(phi), np.degrees(theta)-90, linestyle='none', alpha=0.1)
cen_r = np.sqrt(np.sum(cen_x[:,:3]**2, axis=-1))
cen_phi = np.arctan2(cen_x[:,1], cen_x[:,0])
cen_theta = np.arccos(cen_x[:,2]/cen_r)
axes[i].plot(np.degrees(cen_phi[-1]), np.degrees(cen_theta[-1])-90, linestyle='none', marker='o')
axes[i].set_xlim(-90,90.)
axes[i].set_ylim(-10,10.)
axes[i].set_xlabel("Stream longitude [deg]")
axes[0].set_ylabel(r"Latitude [deg]")
fig.tight_layout()
In [35]:
fig,axes = plt.subplots(2,2,figsize=(10,10),sharex='row',sharey='row')
kwargs = dict(marker='o', linestyle='none', alpha=0.1)
for i,(x,cen_x) in enumerate(zip([slow_w[:,:3],fast_w[:,:3]],[slow_cen_w[:,:3],fast_cen_w[:,:3]])):
ax = axes[0]
ax[i].plot(x[:,0], x[:,1], **kwargs)
ax[i].set_xlim(-25,45)
ax[i].set_ylim(-35,35)
ax[i].set_xticks(range(-20,41,20))
ax[i].set_yticks(range(-20,21,20))
# -----------------------------------------------------------------------
ax = axes[1]
r = np.sqrt(np.sum(x[:,:3]**2, axis=-1))
phi = np.arctan2(x[:,1],x[:,0])
ax[i].plot(np.degrees(phi), r, **kwargs)
cen_r = np.sqrt(np.sum(cen_x[:,:3]**2, axis=-1))
cen_phi = np.arctan2(cen_x[:,1], cen_x[:,0])
ax[i].plot(np.degrees(cen_phi[-1]), cen_r[-1], marker='o', markersize=10, color='r')
ax[i].set_xlim(-180,180.)
ax[i].set_ylim(5,45)
ax[i].set_xticks(range(-180,181,90))
ax[i].set_yticks(range(10,41,10))
fig.tight_layout()
fig.subplots_adjust(wspace=0.15, hspace=0.2)
# fig.savefig("/Users/adrian/papers/posters/AAS225/figures/figure4.pdf", transparent=True)
In [37]:
fig,axes = plt.subplots(1,2,figsize=(11,5.5),sharex=True,sharey=True)
kwargs = dict(marker='o', linestyle='none', alpha=0.05, markersize=7)
colors = ["#1b75bb","#b11f23"]
for i,(x,cen_x) in enumerate(zip([slow_w[:,:3],fast_w[:,:3]],[slow_cen_w[:,:3],fast_cen_w[:,:3]])):
axes[i].plot(x[:,0], x[:,1], **kwargs)
axes[i].plot(cen_x[-100:,0], cen_x[-100:,1], marker=None, linestyle='-', linewidth=4, color=colors[i], zorder=-1)
# axes[i].plot(cen_x[-1,0], cen_x[-1,1], marker='x', markersize=10, markeredgewidth=2, color='k')
axes[i].set_xlim(-20,40)
axes[i].set_ylim(-30,30)
axes[i].set_xticks(range(-15,31,15))
axes[i].set_yticks(range(-15,16,15))
for tick in axes[i].xaxis.get_major_ticks() + axes[i].yaxis.get_major_ticks():
tick.label.set_fontsize(40)
fig.tight_layout()
# fig.savefig("/Users/adrian/papers/posters/AAS225/figures/figure4.pdf", transparent=True)
In [ ]:
In [575]:
from sklearn.neighbors import KernelDensity
In [599]:
kde = KernelDensity(bandwidth=0.5, kernel='exponential', algorithm='ball_tree')
In [617]:
Hs = []
In [38]:
fig,axes = plt.subplots(1,2,figsize=(11,5.75),sharex=True,sharey=True)
kwargs = dict(marker='o', linestyle='none', alpha=0.1)
colors = ["#1b75bb","#b11f23"]
for i,(x,cen_x) in enumerate(zip([slow_w[:,:3],fast_w[:,:3]],[slow_cen_w[:,:3],fast_cen_w[:,:3]])):
# axes[i].plot(x[:,0], x[:,1], **kwargs)
xbins = np.linspace(-20,40,250)
ybins = np.linspace(-30,30,250)
sh = (len(xbins),len(ybins))
if len(Hs) < 2:
kde.fit(x[:,:2])
xgrid,ygrid = map(np.ravel, np.meshgrid(xbins,ybins))
grid = np.vstack((xgrid,ygrid)).T
H = kde.score_samples(grid)
Hs.append(H)
HH = Hs[i].reshape(sh)
# H,xedges,yedges = np.histogram2d(x[:,0], x[:,1], bins=(xbins,ybins))
# HH = np.log(H - H.min())
# HH[~np.isfinite(HH)] = -9999
axes[i].pcolormesh(xgrid.reshape(sh), ygrid.reshape(sh), HH, cmap='Greys',
vmin=-25, vmax=-3)
axes[i].plot(cen_x[-100:,0], cen_x[-100:,1], marker=None, linestyle='-', linewidth=3, color=colors[i], alpha=0.5)
axes[i].set_xlim(xbins.min(),xbins.max())
axes[i].set_ylim(ybins.min(),ybins.max())
axes[i].set_xticks(range(-15,31,15))
axes[i].set_yticks(range(-15,16,15))
for tick in axes[i].xaxis.get_major_ticks() + axes[i].yaxis.get_major_ticks():
tick.label.set_fontsize(20)
fig.tight_layout()
# fig.savefig("/Users/adrian/papers/posters/AAS225/figures/figure4.pdf", transparent=True)
In [42]:
fig,axes = plt.subplots(1,2,figsize=(11,5.5),sharex=True,sharey=True)
kwargs = dict(marker='o', linestyle='none', alpha=0.05, markersize=7)
colors = ["#1b75bb","#b11f23"]
for i,(x,cen_x) in enumerate(zip([slow_w[:,:3],fast_w[:,:3]],[slow_cen_w[:,:3],fast_cen_w[:,:3]])):
axes[i].plot(x[:,0], x[:,1], **kwargs)
axes[i].plot(cen_x[-100:,0], cen_x[-100:,1], marker=None, linestyle='-', linewidth=4, color=colors[i], zorder=-1)
# axes[i].plot(cen_x[-1,0], cen_x[-1,1], marker='x', markersize=10, markeredgewidth=2, color='k')
axes[i].set_xlim(15,25)
axes[i].set_ylim(-10,10)
fig.tight_layout()
# fig.savefig("/Users/adrian/papers/posters/AAS225/figures/figure4.pdf", transparent=True)
In [ ]: