In [1]:
from __future__ import division, print_function
import os
# Third-party
import astropy.units as u
from astropy.constants import G
import matplotlib.pyplot as plt
import numpy as np
import yt
%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
G = G.decompose(galactic).value
In [2]:
figure_path = "/Users/adrian/projects/morphology-paper/figures/"
potential = gp.load("/Users/adrian/projects/morphology/potentials/triaxial-NFW.yml")
r_s = potential.parameters['r_s']
potential
Out[2]:
In [61]:
lim = 75
grid = np.linspace(-100,100,1000)
cmap = 'Blues_r'
vals = potential.value(np.array([[1.,0,0],[lim,0,0]]))
levels = np.linspace(vals.min(), vals.max(), 8)
# levels = np.logspace(vals.min(), vals.max(), 10)
# ------------------------------------------------
fig,axes = plt.subplots(2,2, figsize=(6,6), sharex=True, sharey=True)
fig.set_facecolor('w')
fig = potential.plot_contours(grid=(grid,grid,0), ax=axes[0,0], cmap=cmap, levels=levels)
fig = potential.plot_contours(grid=(grid,0,grid), ax=axes[1,0], cmap=cmap, levels=levels)
fig = potential.plot_contours(grid=(0,grid,grid), ax=axes[1,1], cmap=cmap, levels=levels)
axes[0,0].set_xlim(-75,75)
axes[0,0].set_ylim(*axes[0,0].get_xlim())
axes[0,0].set_ylabel(r"$y/r_s$")
axes[1,0].set_xlabel(r"$x/r_s$")
axes[1,0].set_ylabel(r"$z/r_s$")
axes[1,1].set_xlabel(r"$y/r_s$")
# Hack to hide hole
axes[0,0].plot(0,0,marker='o',zorder=-1000,c='k')
axes[1,0].plot(0,0,marker='o',zorder=-1000,c='k')
axes[1,1].plot(0,0,marker='o',zorder=-1000,c='k')
tk = np.arange(-3,3+1,1)
for ax in axes.flat:
ax.set_aspect('equal')
ax.xaxis.set_ticks(tk * r_s)
ax.yaxis.set_ticks(tk * r_s)
ax.xaxis.set_ticklabels(["{0:d}".format(x) for x in tk])
ax.yaxis.set_ticklabels(["{0:d}".format(x) for x in tk])
fig.text(0.77,0.75,"potential", fontsize=26, ha="center")
fig.tight_layout()
fig.subplots_adjust(wspace=0., hspace=0.)
axes[0,1].set_visible(False)
fig.savefig(os.path.join(figure_path, "potential.pdf"))
In [4]:
from scipy.misc import derivative
In [5]:
def func(x):
return potential.value([x,0.,0.])[0]
dens_grid = np.linspace(1,250,128)
num_dens = np.array([-derivative(func, x, 0.1, n=2) / (4*np.pi*potential.G) for x in dens_grid])
num_dens
Out[5]:
In [6]:
xyz = np.zeros((len(dens_grid),3))
xyz[:,0] = dens_grid
func_dens = potential.density(xyz)
func_dens
Out[6]:
In [7]:
func_dens / num_dens
Out[7]:
In [8]:
plt.loglog(dens_grid, func_dens)
plt.loglog(dens_grid, num_dens)
Out[8]:
In [9]:
levels
Out[9]:
In [71]:
lim = 75
grid = np.linspace(-100,100,1000)
cmap = 'Blues'
vals = potential.density(np.array([[1.,0,0],[lim,0,0]]))
# levels = np.linspace(vals.min(), vals.max(), 16)
# levels = np.logspace(np.log10(vals.min()), np.log10(vals.max()), 16)
levels = np.logspace(4, 7., 8)
# ------------------------------------------------
fig,axes = plt.subplots(2,2, figsize=(6,6), sharex=True, sharey=True)
fig.set_facecolor('w')
fig = potential.plot_densty_contours(grid=(grid,grid,0), ax=axes[0,0], cmap=cmap, levels=levels, extend='both')
fig = potential.plot_densty_contours(grid=(grid,0,grid), ax=axes[1,0], cmap=cmap, levels=levels, extend='both')
fig = potential.plot_densty_contours(grid=(0,grid,grid), ax=axes[1,1], cmap=cmap, levels=levels, extend='both')
axes[0,0].set_xlim(-75,75)
axes[0,0].set_ylim(*axes[0,0].get_xlim())
axes[0,0].set_ylabel(r"$y/r_s$")
axes[1,0].set_xlabel(r"$x/r_s$")
axes[1,0].set_ylabel(r"$z/r_s$")
axes[1,1].set_xlabel(r"$y/r_s$")
# Hack to hide hole
axes[0,0].plot(0,0,marker='o',zorder=-1000,c='k')
axes[1,0].plot(0,0,marker='o',zorder=-1000,c='k')
axes[1,1].plot(0,0,marker='o',zorder=-1000,c='k')
tk = np.arange(-3,3+1,1)
for ax in axes.flat:
ax.set_aspect('equal')
ax.xaxis.set_ticks(tk * r_s)
ax.yaxis.set_ticks(tk * r_s)
ax.xaxis.set_ticklabels(["{0:d}".format(x) for x in tk])
ax.yaxis.set_ticklabels(["{0:d}".format(x) for x in tk])
fig.text(0.77,0.75,"density", fontsize=26, ha="center")
fig.tight_layout()
fig.subplots_adjust(wspace=0., hspace=0.)
axes[0,1].set_visible(False)
fig.savefig(os.path.join(figure_path, "density-contours.pdf"))
In [8]:
maxx = 75.
grid1d = np.linspace(-maxx,maxx,128)
xx,yy,zz = np.meshgrid(grid1d,grid1d,grid1d)
In [9]:
xyz = np.vstack((xx.flat,yy.flat,zz.flat)).T
pot_val = np.abs(potential.value(xyz).reshape(xx.shape))
In [10]:
data = dict(potential=(pot_val, "kpc**2/Myr**2"))
bbox = np.array([[-maxx, maxx], [-maxx, maxx], [-maxx, maxx]])
ds = yt.load_uniform_grid(data, pot_val.shape, bbox=bbox, length_unit="kpc")
In [15]:
# Choose a field
field = 'potential'
# Find the bounds in log space of for your field
dd = ds.all_data()
mi, ma = dd.quantities.extrema(field)
# Instantiate the ColorTransferfunction.
tf = yt.ColorTransferFunction((mi, ma))
# Set up the camera parameters: center, looking direction, width, resolution
c = (ds.domain_right_edge + ds.domain_left_edge)/2.0
W = (100.,'kpc')
N = 512
for i,phi in enumerate(np.linspace(0,np.pi/2.,5)):
cam_dir = np.array([75.*np.cos(phi), 0., 75*np.sin(phi)])
# Create a camera object
cam = ds.camera(c, cam_dir, W, N, tf, fields=[field], log_fields=[False], no_ghost=False)
# tf.add_layers(5, mi=mi, ma=yt.YTQuantity(-0.15,mi.units), colormap='RdBu')
tf.add_layers(10, mi=mi, ma=ma, colormap='RdBu')
im = cam.snapshot()
im.write_png('/Users/adrian/projects/morphology/plots/triaxial-NFW/{0:03d}.png'.format(i))
# cam.show()
In [ ]: