In [ ]:
from __future__ import division, print_function
import os
import sys
from collections import OrderedDict
# Third-party
import astropy.coordinates as coord
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
pl.style.use('apw-notebook')
%matplotlib inline
# Custom
import gala.dynamics as gd
import gala.integrate as gi
import gala.potential as gp
from gala.units import galactic
from scipy.misc import factorial
# from ophiuchus import barred_mw, static_mw
import ophiuchus.potential as op
plotpath = "/Users/adrian/projects/ophiuchus-paper/figures/"
if not os.path.exists(plotpath):
os.mkdir(plotpath)
In [ ]:
barred_mw = op.load_potential("barred_mw_4")
static_mw = op.load_potential("static_mw")
In [ ]:
# transform from H&O 1992 coefficients to Lowing 2011 coefficients
nlms = np.array([[0,0,0],
[1,0,0],
[2,0,0],
[3,0,0],
[0,2,0],
[1,2,0],
[2,2,0],
[0,2,2],
[1,2,2],
[2,2,2],
[0,4,0],
[1,4,0],
[0,4,2],
[1,4,2],
[0,4,4],
[1,4,4],
[0,6,0],
[0,6,2],
[0,6,4],
[0,6,6]])
_Snlm = np.array([1.509,-0.086,-0.033,-0.02,-2.606,
-0.221,-0.001,0.665,0.129,0.006,6.406,
1.295,-0.66,-0.14,0.044,-0.012,-5.859,
0.984,-0.03,0.001])
NEW_S = _Snlm.copy()
for i,(n,l,m) in zip(range(len(_Snlm)), nlms):
if l != 0:
fac = np.sqrt(4*np.pi) * np.sqrt((2*l+1) / (4*np.pi) * factorial(l-m) / factorial(l+m))
NEW_S[i] /= fac
In [ ]:
nmax = 3
lmax = 6
Snlm = np.zeros((nmax+1,lmax+1,lmax+1))
for (n,l,m),A in zip(nlms,NEW_S):
Snlm[n,l,m] = A
In [ ]:
static_mw
In [ ]:
barred_mw
In [ ]:
# barpars = barred_mw.parameters.copy()
# barpars['halo']['q_z'] = 1.
# barpars['spheroid']['c'] = 0.2
# barpars['spheroid']['m'] = 5E9
# barpars['disk']['m'] = 4E10
# barpars['bar']['r_s'] = 1.2
# barpars['bar']['m'] = barpars['bar']['m']
# barred_mw = op.OphiuchusPotential(**barpars)
# stapars = static_mw.parameters.copy()
# stapars['halo']['q_z'] = 1.
# stapars['spheroid']['c'] = 0.3
# stapars['spheroid']['m'] = 1.2E10
# stapars['disk']['m'] = 6E10
# static_mw = op.OphiuchusPotential(**stapars)
In [ ]:
potential_classes = OrderedDict()
potential_classes['disk'] = gp.MiyamotoNagaiPotential
potential_classes['halo'] = gp.FlattenedNFWPotential
potential_classes['bar'] = op.WangZhaoBarPotential
potential_classes['spheroid'] = gp.HernquistPotential
In [ ]:
(0.19*u.kpc/u.Myr).to(u.km/u.s)
In [ ]:
ix = 0
xyz = np.zeros((3,128))
xyz[ix] = np.linspace(0.,10.,xyz.shape[1])
for pot in [static_mw, barred_mw]:
Menc = pot.mass_enclosed(xyz)
pl.loglog(xyz[ix], Menc, marker='')
pl.axvline(1)
pl.axhline(1E10)
In [ ]:
def density_on_grid(potential, t=0., grid_lim=(-15,15), ngrid=128):
grid = np.linspace(grid_lim[0], grid_lim[1], ngrid)
xyz = np.vstack(map(np.ravel, np.meshgrid(grid,grid,grid)))
# val = np.zeros((ngrid*ngrid*ngrid,))
val = potential.density(xyz, t=t).value
val[np.isnan(val)] = val[np.isfinite(val)].min()
val[val < 0] = 1.
gridx = xyz[0].reshape(ngrid,ngrid,ngrid)[:,:,0]
gridy = xyz[1].reshape(ngrid,ngrid,ngrid)[:,:,0]
return gridx, gridy, val
In [ ]:
ngrid = 128
xx,yy,barred_dens = density_on_grid(barred_mw, ngrid=ngrid)
xx,yy,static_dens = density_on_grid(static_mw, ngrid=ngrid)
In [ ]:
def side_by_side_surface_dens(xx, yy, dens):
ngrid = xx.shape[0]
fig,axes = pl.subplots(1, 2, figsize=(8,4),
sharex=True, sharey=True)
axes[0].pcolormesh(xx, yy, dens.reshape(ngrid,ngrid,ngrid).sum(axis=2),
cmap='Greys_r',
norm=mpl.colors.LogNorm(),
vmin=1E7, vmax=5E9)
axes[0].text(-8., 0, r"$\odot$", ha='center', va='center', fontsize=18, color='w')
axes[1].pcolormesh(xx, yy, dens.reshape(ngrid,ngrid,ngrid).sum(axis=0).T,
cmap='Greys_r',
norm=mpl.colors.LogNorm(),
vmin=1E7, vmax=5E9)
axes[0].set_xlim(xx.min(), xx.max())
axes[0].set_ylim(yy.min(), yy.max())
# TODO: fix the damn aspect ratio
# for ax in axes:
# ax.set_aspect('equal')
fig.tight_layout()
return fig
In [ ]:
fig = side_by_side_surface_dens(xx, yy, barred_dens)
fig = side_by_side_surface_dens(xx, yy, static_dens)
In [ ]:
def side_by_side_contour_plots(xx, yy, dens, levels=10**np.arange(7,12,0.25)):
ngrid = xx.shape[0]
fig,axes = pl.subplots(1,2,figsize=(7.8,4),sharex=True,sharey=True)
im = axes[0].contour(xx, yy, dens.reshape(ngrid,ngrid,ngrid).sum(axis=2),
colors='k',
levels=levels,
rasterized=True)
axes[0].text(-8., 0, r"$\odot$", ha='center', va='center', fontsize=18)
_ = axes[1].contour(xx, yy, dens.reshape(ngrid,ngrid,ngrid).sum(axis=1).T,
colors='k',
levels=levels,
rasterized=True)
# fig.subplots_adjust(bottom=0.2, right=0.85, wspace=0.25)
for ax in axes:
ax.xaxis.set_ticks([-10,0,10])
ax.yaxis.set_ticks([-10,0,10])
axes[0].set_xlabel("$x$ [kpc]")
axes[0].set_ylabel("$y$ [kpc]")
axes[1].set_xlabel("$y$ [kpc]")
axes[1].set_ylabel("$z$ [kpc]")
axes[0].set_xlim(xx.min(), xx.max())
axes[0].set_ylim(yy.min(), yy.max())
fig.tight_layout()
return fig
In [ ]:
barred_fig = side_by_side_contour_plots(xx, yy, barred_dens)
static_fig = side_by_side_contour_plots(xx, yy, static_dens)
# barred_fig.savefig(os.path.join(plotpath, "barred-surface-density-contour.pdf"), bbox_inches='tight')
# barred_fig.savefig(os.path.join(plotpath, "barred-surface-density-contour.png"), dpi=400, bbox_inches='tight')
# static_fig.savefig(os.path.join(plotpath, "static-surface-density-contour.pdf"), bbox_inches='tight')
# static_fig.savefig(os.path.join(plotpath, "static-surface-density-contour.png"), dpi=400, bbox_inches='tight')
In [ ]:
ngrid = 65
grid = np.linspace(-2,2,ngrid)
xyz = np.vstack(map(np.ravel, np.meshgrid(grid,grid,grid)))
val2 = np.zeros((ngrid*ngrid*ngrid,))
# for k in potentials.keys():
# val += potentials[k].density(xyz)
val2 += potentials['bar'].density(xyz)
val2[np.isnan(val2)] = val2[np.isfinite(val2)].max()
In [ ]:
surf_dens = (val2.reshape(ngrid,ngrid,ngrid).sum(axis=1).T*u.Msun/(u.kpc**2)/ngrid).to(u.Msun/u.pc**2)
pl.figure(figsize=(6,3))
pl.contourf(xyz[0].reshape(ngrid,ngrid,ngrid)[:,:,0],
xyz[1].reshape(ngrid,ngrid,ngrid)[:,:,0],
surf_dens.value,
norm=mpl.colors.LogNorm(),
levels=np.logspace(1., 4, 8),
cmap='Blues')
# cmap='Greys_r',
# norm=mpl.colors.LogNorm(),
# vmin=5E8, vmax=5E10)
pl.xlim(-2,2)
pl.ylim(-1.1,1.1)
pl.colorbar()
pl.tight_layout()
In [ ]:
def circ_vel_plot(potential, name):
""" name = barred, static """
rr = np.linspace(0.1, 20., 1024)
xyz = np.zeros((3, len(rr)))
xyz[0] = rr
potentials = OrderedDict()
for k,P in potential_classes.items():
potentials[k] = P(units=galactic, **potential.parameters[k])
# vcirc = (np.sqrt(potential.G * potential.mass_enclosed(xyz) / rr)*u.kpc/u.Myr).to(u.km/u.s).value
vcirc = (np.sqrt(potential.G * np.sum([p.mass_enclosed(xyz) for p in potentials.values()], axis=0) / rr)*u.kpc/u.Myr).to(u.km/u.s).value
fig,ax = pl.subplots(1,1,figsize=(6,5))
ax.plot(rr, vcirc, marker='', lw=3.)
styles = dict(
halo=dict(lw=2, ls='-.'),
bar=dict(lw=3., ls=':'),
spheroid=dict(lw=3., ls=':'),
disk=dict(lw=2., ls='--')
)
for k,p in potentials.items():
if k != 'halo' and potential.parameters[k]['m'] == 0:
continue
if k == 'bar':
continue
if name == 'static':
disk_other = 'Spher'
elif name == 'barred':
disk_other = 'Bar+Spher'
vc = (np.sqrt(potential.G * p.mass_enclosed(xyz).value / rr)*u.kpc/u.Myr).to(u.km/u.s).value
if name == 'barred' and k == 'spheroid':
menc_sph = p.mass_enclosed(xyz)
p = potentials['bar']
vc = (np.sqrt(potential.G * (menc_sph + p.mass_enclosed(xyz)).value / rr)*u.kpc/u.Myr).to(u.km/u.s).value
label = 'Bar+Spheroid'
else:
label = k.capitalize()
ax.plot(rr, vc, marker='', label=label, **styles[k])
if name == 'barred':
vc = (np.sqrt(potential.G * (potentials['spheroid'].mass_enclosed(xyz)+potentials['bar'].mass_enclosed(xyz)+potentials['disk'].mass_enclosed(xyz)).value / rr)*u.kpc/u.Myr).to(u.km/u.s).value
ax.plot(rr, vc, marker='', label='Disk+Bar+Spher', lw=2.)
else:
vc = (np.sqrt(potential.G * (potentials['spheroid'].mass_enclosed(xyz)+potentials['disk'].mass_enclosed(xyz)).value / rr)*u.kpc/u.Myr).to(u.km/u.s).value
ax.set_xlabel("$R$ [kpc]")
ax.set_ylabel(r"$v_c$ [${\rm km}\,{\rm s}^{-1}$]")
ax.legend(loc='upper right', fontsize=12)
ax.set_ylim(0,300)
# ax.set_ylim(150,300)
# ax.axhline(220, alpha=0.2, lw=1.)
# ax.axvline(8., color='#cccccc', lw=2., zorder=-100)
rcolor = '#dddddd'
rect = mpl.patches.Rectangle((0.,215), rr.max(), 20., zorder=-100, color=rcolor)
ax.add_patch(rect)
rect2 = mpl.patches.Rectangle((8.,0), 0.3, ax.get_ylim()[1], zorder=-100, color=rcolor)
ax.add_patch(rect2)
fig.tight_layout()
return fig
In [ ]:
fig = circ_vel_plot(barred_mw, 'barred')
# fig.savefig(os.path.join(plotpath, "barred-circ-vel.pdf"))
# fig.savefig(os.path.join(plotpath, "barred-circ-vel.png"), dpi=400)
In [ ]:
fig = circ_vel_plot(static_mw, name='static')
# fig.savefig(os.path.join(plotpath, "static-circ-vel.pdf"))
# fig.savefig(os.path.join(plotpath, "static-circ-vel.png"), dpi=400)
In [ ]:
fig,axes = pl.subplots(2,2,figsize=(9,8.5),sharex='col')
# Circular velocity
styles = dict(
halo=dict(lw=2, ls='-.'),
bar=dict(lw=3., ls=':'),
spheroid=dict(lw=3., ls=':'),
disk=dict(lw=2., ls='--')
)
# Contour
levels = 10**np.arange(7,12,0.25)
rr = np.linspace(0.1, 22., 1024)
fac = static_mw.G / rr
xyz = np.zeros((3, len(rr)))
xyz[0] = rr
for i,(name,pot,dens) in enumerate(zip(['barred','static'], [barred_mw, static_mw],[barred_dens, static_dens])):
# Circular velocity
ax = axes[i,0]
potentials = OrderedDict()
for k,P in potential_classes.items():
potentials[k] = P(units=galactic, **pot.parameters[k])
# vcirc = (np.sqrt(potential.G * potential.mass_enclosed(xyz) / rr)*u.kpc/u.Myr).to(u.km/u.s).value
vcirc = (np.sqrt(pot.G * np.sum([p.mass_enclosed(xyz) for p in potentials.values()], axis=0) / rr)*u.kpc/u.Myr)\
.to(u.km/u.s).value
ax.plot(rr, vcirc, marker='', lw=3.)
menc = dict()
for k,p in potentials.items():
menc[k] = p.mass_enclosed(xyz)
# Halo
vc = np.sqrt(fac * menc['halo'].value)
ax.plot(rr, (vc*u.kpc/u.Myr).to(u.km/u.s),
marker='', label='Halo', **styles['halo'])
# disk, etc.
if name == 'static':
vc = np.sqrt(fac * (menc['disk']+menc['spheroid']).value)
ax.plot(rr, (vc*u.kpc/u.Myr).to(u.km/u.s),
marker='', label='Disk+Sph', **styles['disk'])
elif name == 'barred':
vc = np.sqrt(fac * (menc['disk']+menc['spheroid']+menc['bar']).value)
ax.plot(rr, (vc*u.kpc/u.Myr).to(u.km/u.s),
marker='', label='Disk+Sph+Bar', **styles['disk'])
ax.legend(loc='upper right', fontsize=12)
ax.set_ylim(0,300)
# ax.set_ylim(150,300)
# ax.axhline(220, alpha=0.2, lw=1.)
# ax.axvline(8., color='#cccccc', lw=2., zorder=-100)
rcolor = '#dddddd'
rect = mpl.patches.Rectangle((0.,215), rr.max(), 22., zorder=-100, color=rcolor)
ax.add_patch(rect)
rect2 = mpl.patches.Rectangle((8.,0), 0.3, ax.get_ylim()[1], zorder=-100, color=rcolor)
ax.add_patch(rect2)
# Surface density
ngrid = xx.shape[0]
ax = axes[i,1]
im = ax.contour(xx, yy, dens.reshape(ngrid,ngrid,ngrid).sum(axis=2),
colors='k', levels=levels, rasterized=True)
ax.text(-8., 0, r"$\odot$", ha='center', va='center', fontsize=18)
ax.xaxis.set_ticks([-10,0,10])
ax.yaxis.set_ticks([-10,0,10])
ax.set_xlim(xx.min(), xx.max())
ax.set_ylim(yy.min(), yy.max())
if i == 0:
ax = axes[0,0]
ax.text(8.4, 40, r'$R_\odot$', fontsize=18, color='#666666')
# ax.annotate(r'$R_\odot$', xy=(8.3, 50), xytext=(12, 75.),
# fontsize=18,
# xycoords='data', textcoords='data',
# arrowprops=dict(arrowstyle="fancy",
# fc="0.6", ec="none",
# patchB=rect2,
# connectionstyle="angle3,angleA=0,angleB=90"),
# )
axes[0,0].text(1, 260, "Barred", fontsize=24, fontstyle='italic', ha='left')
axes[1,0].text(1, 260, "Static", fontsize=24, fontstyle='italic', ha='left')
axes[1,0].set_xlabel("$R$ [kpc]")
axes[1,1].set_xlabel("$x$ [kpc]")
axes[0,0].set_ylabel(r"$v_c$ [${\rm km}\,{\rm s}^{-1}$]")
axes[1,0].set_ylabel(r"$v_c$ [${\rm km}\,{\rm s}^{-1}$]")
axes[0,0].set_xlim(0,22)
axes[0,1].set_ylabel("$y$ [kpc]")
axes[1,1].set_ylabel("$y$ [kpc]")
axes[0,1].yaxis.set_label_position('right')
axes[1,1].yaxis.set_label_position('right')
axes[0,1].yaxis.tick_right()
axes[1,1].yaxis.tick_right()
axes[1,1].set_aspect('equal')
fig.tight_layout()
# fig.savefig(os.path.join(plotpath, "potentials-four.pdf"))
# fig.savefig(os.path.join(plotpath, "potentials-four.png"), dpi=400)
In [ ]:
pot = op.WangZhaoBarPotential(**barred_mw.parameters['bar'])
T = (2*np.pi/(60*u.km/u.s/u.kpc)).to(u.Myr).value
for time in np.linspace(0.,T/4,4):
xx,yy,_dens = density_on_grid(pot, t=time, ngrid=64)
fig = side_by_side_surface_dens(xx, yy, _dens)
In [ ]:
pars = barred_mw.parameters['bar'].copy()
pars['alpha'] = 0.
pot = op.WangZhaoBarPotential(**pars)
X = np.linspace(-15,15,256)
_xyz = np.zeros((X.size,3))
_xyz[:,0] = X
along_x = pot.acceleration(_xyz)[:,0]
_xyz = np.zeros((X.size,3))
_xyz[:,1] = X
along_y = pot.acceleration(_xyz)[:,1]
In [ ]:
pl.plot(X, np.abs(along_x))
pl.plot(X, np.abs(along_y))
In [ ]:
engrid = 32
derp = np.linspace(-15,15,engrid)
xy = np.vstack(map(np.ravel, np.meshgrid(derp,derp))).T
xyz = np.zeros((len(xy),3))
xyz[:,[0,2]] = xy
dens = pot.density(xyz, t=0)
dens[np.isnan(dens)] = dens[np.isfinite(dens)].max()
xx = xyz[:,0].reshape(engrid,engrid)
yy = xyz[:,2].reshape(engrid,engrid)
pl.figure(figsize=(5,5))
pl.contour(xx, yy, dens.reshape(engrid,engrid),
colors='k', rasterized=True)
In [ ]: