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()))
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.
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())
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.)
In [30]:
d[bad_ix][5][0]
Out[30]:
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]:
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]:
In [518]:
_w0 = w0[_ix]
_freqdiff = all_freq_diff[_ix]
In [519]:
list(_w0[_freqdiff.argmin()])
Out[519]:
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]:
In [562]:
fast_w0.tolist(),slow_w0.tolist(),
Out[562]:
In [349]:
fast_w0 = _w0[7].copy()
print(_freqdiff[7])
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())
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)
In [353]:
(freqs1[1] - freqs1[0]) / freqs1[0]
Out[353]:
In [526]:
T1 = (2*np.pi/freqs1)
2*nsteps / T1
Out[526]:
In [572]:
2*28000 / T1
Out[572]:
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]:
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())
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)
In [576]:
print((freqs1[1] - freqs1[0]) / (dt*nsteps/2.))
print((freqs2[1] - freqs2[0]) / (dt*nsteps/2.))
In [528]:
T2 = (2*np.pi/freqs2)
2*nsteps / T2
Out[528]:
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]:
In [359]:
print(slow_w0)
print(fast_w0)
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]:
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()
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]:
In [ ]: