This notebook generates all simulation results figures for my current paper. All simulations run on stampede.tacc.utexas.edu
%matplotlib inline
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
mpl.rc('font', size=20.)
mpl.rc('font', family='serif')
mpl.rc('text', usetex=True)
import pyGadget
Use the final snapshot of the vanilla simulation.
sims = ['vanilla', 'xr_tau_J0', 'xr_tau_J1', 'xr_tau2_J2', 'xr_tau_J3', 'XR_sfr_1e-1', 'XR_sfr_1e-2', 'XR_sfr_1e-3']
n0 = [355, 306, 327, 227, 235, 200, 201, 269]
n5k = [1857, 1546, 1852, 1758, 1687, 1616, 1900]
t0 = '_t0'
t5k = '_t5k'
n, tag = n0, t0
i = 0
sim = pyGadget.sim.Simulation('stampede/'+sims[i])
sim.refine_by_mass(False)
sim.set_coordinate_system('physical')
snap = sim.load_snapshot(n[i])
imzoom = []
for scale in ['5376pc', '1000pc', '10pc', '100pc']:
imzoom.append(pyGadget.visualize.project(snap, 'ndensity', scale, 'xz', centering='avg'))
#snap.close()
from mpl_toolkits.axes_grid1 import ImageGrid
scales = ['140 kpc (comoving)', '1 kpc (physical)', '10 pc (physical)', '100 pc (physical)']
ratio = [.1788, .1, None, .1]
zoom = ['right', 'down', None, 'left']
clims = [(-2.5,1.5), (-2.,2.), (1.5,7.5), (-0.5,5.)]
ticks = [(-2,-1,0,1), (-1,0,1), (2,3,4,5,6,7), (0,1,2,3,4)]
cpad = [-17, -17, -15, -16]
clabel = [False, True, False, True]
bbox_props = dict(boxstyle="round", fc="k", ec="k", alpha=0.5)
zc = 'w'
zls = '--'
zlw = 1.5
fig = plt.figure(1, (12., 12.), dpi=600)
grid = ImageGrid(fig, 111, # similar to subplot(111)
nrows_ncols = (2, 2), # creates 2x2 grid of axes
axes_pad=0.0, # pad between axes in inch.
cbar_mode = 'each', cbar_size='7%', cbar_pad=0.
)
for i in range(4):
x = imzoom[i][0]
y = imzoom[i][1]
im = imzoom[i][2]
ax = grid[i]
img = ax.imshow(im, cmap=plt.cm.bone, origin='lower')
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
img.set_clim(clims[i])
cb = plt.colorbar(img, cax=grid.cbar_axes[i])
cb.set_ticks(ticks[i])
cb.ax.tick_params(left='on', pad=cpad[i],
labelsize=15, labelcolor='k', labelleft='on', labelright='off')
if clabel[i]: cb.set_label('Log Number Density [cm$^{-3}$]')
ax.text(0.5, 0.025, scales[i], color='w', ha='center', va='bottom', size=12,
transform=grid[i].transAxes, bbox=bbox_props)
if ratio[i]:
axmin, axmax = ax.get_xlim()
axlength = axmax - axmin
mid = axlength/2
s = ratio[i] * axlength
s00 = [mid - s/2, mid - s/2]
s01 = [mid - s/2, mid + s/2]
s11 = [mid + s/2, mid + s/2]
ax.add_line(plt.Line2D(s00, s01, c=zc, lw=zlw))
ax.add_line(plt.Line2D(s11, s01, c=zc, lw=zlw))
ax.add_line(plt.Line2D(s01, s00, c=zc, lw=zlw))
ax.add_line(plt.Line2D(s01, s11, c=zc, lw=zlw))
if zoom[i] == 'right':
ax.add_line(plt.Line2D([mid+s/2, axmax], [mid+s/2, axmax], c=zc, lw=zlw, ls=zls))
ax.add_line(plt.Line2D([mid+s/2, axmax], [mid-s/2, axmin], c=zc, lw=zlw, ls=zls))
elif zoom[i] == 'down':
ax.add_line(plt.Line2D([mid-s/2, axmin], [mid-s/2, axmin], c=zc, lw=zlw, ls=zls))
ax.add_line(plt.Line2D([mid+s/2, axmax], [mid-s/2, axmin], c=zc, lw=zlw, ls=zls))
elif zoom[i] == 'left':
ax.add_line(plt.Line2D([mid-s/2, axmin], [mid+s/2, axmax], c=zc, lw=zlw, ls=zls))
ax.add_line(plt.Line2D([mid-s/2, axmin], [mid-s/2, axmin], c=zc, lw=zlw, ls=zls))
plt.show()
#fig.savefig('figures/structure/structure-'+sim.name.split('/')[-1]+tag+'.png', bbox_inches='tight', dpi=100)
fig.savefig('figures/structure/structure-'+sim.name.split('/')[-1]+tag+'.pdf', bbox_inches='tight', dpi=150)
For each sim, picking the snapshot just prior to the formation of the first sink.
simV = pyGadget.sim.Simulation('stampede/vanilla', track_sinks=True)
sim1 = pyGadget.sim.Simulation('stampede/xr_tau2_J0', track_sinks=True)
sim2 = pyGadget.sim.Simulation('stampede/xr_tau2_J1', track_sinks=True)
sim3 = pyGadget.sim.Simulation('stampede/xr_tau2_J2', track_sinks=True)
snapV = simV.load_snapshot(1900)
snap1 = sim1.load_snapshot(1778)
snap2 = sim2.load_snapshot(1669)
snap3 = sim3.load_snapshot(1727)
import copy
snaplist = [snapV, snap1, snap2, snap3]
imlist = []
sinklist = []
scale = '15000AU'
#shifty = [None, None, 5000, -2000, None, None, None, None]
shifty = [None, None, None, None, None, None, None, None]
face = [[('x', 0.29518), ('z', 0.825), ('x',np.pi/2)],
[('y', 0.6346), ('z', 2.03), ('x',np.pi/2)],
[('x', 1.865), ('z', 2.919), ('x',np.pi/2)],
[('x', 1.7136), ('z', 0.18), ('x',np.pi/2)]]
edge = [[('x', 0.29518), ('z', 0.825)],
[('y', 0.6346), ('z', 2.03)],
[('x', 1.865), ('z', 2.919)],
[('x', 1.7136), ('z', 0.18)]
]
ocount = 0
for view in [face, edge]:
count = 0
for snap in snaplist:
imlist.append(pyGadget.visualize.project(snap, 'ndensity', scale, view[count], centering='avg',
depth=2., shifty=shifty[ocount], dens_lim=None))
sinklist.append(copy.deepcopy(snap.sinks))
count += 1
ocount += 1
# snap.close()
from mpl_toolkits.axes_grid1 import ImageGrid
sim = ['J = 0', 'J = J$_{0}$', 'J = 10 J$_{0}$', 'J = 100 J$_{0}$']
bbox_props = dict(boxstyle="round", fc="k", ec="k", alpha=0.5)
ticks = [(7,8,9,10,11),(6,7,8,9,10)]
fig = plt.figure(1, (20, 8), dpi=600)
grid = ImageGrid(fig, 111, # similar to subplot(111)
nrows_ncols = (2, 4), # creates 4x2 grid of axes
axes_pad=0.0, # pad between axes in inch.
cbar_mode = 'edge', cbar_location = 'right', cbar_size='7%', cbar_pad=0.0
)
for i in range(8):
x = imlist[i][0]
y = imlist[i][1]
im = imlist[i][2]
ax = grid[i]
img = ax.imshow(im, extent=[x.min(),x.max(),y.min(),y.max()], cmap=plt.cm.bone, origin='lower')
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
img.set_clim((6.5,10.75))
cb = plt.colorbar(img, cax=grid.cbar_axes[i])
#cb.set_ticks(ticks[i/4])
cb.set_ticks((7,8,9,10))
cb.ax.tick_params(left='on', labelsize=15, labelcolor='k')
#if clabel[i]:
cb.set_label('Log Number Density \n [cm$^{-3}$]')
for sink in sinklist[i]:
#mscale = sink.mass*6./27. + .33
mscale = np.log(sink.mass) +1
ax.plot(sink.x, sink.y, 'ko', ms=mscale, mew=1)
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
if i > 3:
cb.set_ticks(ticks[i/4])
ax.text(0.5, 0.025, sim[i-4], color='w', ha='center', va='bottom', size=18,
transform=grid[i].transAxes, bbox=bbox_props)
plt.show()
fig.savefig('figures/structure/disks.png', bbox_inches='tight', dpi=100)
reload(pyGadget.visualize)
import copy
snaplist = [snapV, snap1, snap2, snap3]
imlist = []
sinklist = []
scale = '15000AU'
imscale=['log','log','linear']
shifty = [None, None, 5000, -2000]
face = [[('x', 0.29518), ('z', 0.825), ('x',np.pi/2)],
[('y', 0.6346), ('z', 2.03), ('x',np.pi/2)],
[('x', 1.865), ('z', 2.919), ('x',np.pi/2)],
[('x', 1.7136), ('z', 0.18), ('x',np.pi/2)]]
for property in ['ndensity', 'temp', 'h2frac']:
count = 0
for snap in snaplist:
imlist.append(pyGadget.visualize.project(snap, property, scale, face[count], centering='avg', depth=.05,
shifty=shifty[count], imscale=imscale[count/4]))
sinklist.append(copy.deepcopy(snap.sinks))
count += 1
# snap.close()
from mpl_toolkits.axes_grid1 import ImageGrid
sim = ['J = 0', 'J = J$_{0}$', 'J = 10 J$_{0}$', 'J = 100 J$_{0}$']
colormap = [plt.cm.RdGy_r, plt.cm.afmhot, plt.cm.Blues]
color_lims = [(6.5,11), (1.8,3.8), (-3,.5)]
labels = ['Log Number Density \n [cm$^{-3}$]', 'Log Gas Temperature\n [K]', 'H$_2$ Fraction']
bbox_props = dict(boxstyle="round", fc="k", ec="k", alpha=0.5)
ticks = [(7,8,9,10,11),(2.0,2.4,2.8,3.2,3.6,4),(-3,-2,-1,0)]
fig = plt.figure(1, (20, 12), dpi=600)
grid = ImageGrid(fig, 111, # similar to subplot(111)
nrows_ncols = (3, 4), # creates 4x2 grid of axes
axes_pad=0.0, # pad between axes in inch.
cbar_mode = 'each', cbar_location = 'right', cbar_size='7%', cbar_pad=0.0
)
for i in range(12):
x = imlist[i][0]
y = imlist[i][1]
im = imlist[i][2]
ax = grid[i]
img = ax.imshow(im, extent=[x.min(),x.max(),y.min(),y.max()], cmap=colormap[i/4], origin='lower')
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
img.set_clim(color_lims[i/4])
cb = plt.colorbar(img, cax=grid.cbar_axes[i])
cb.set_ticks(ticks[i/4])
cb.ax.tick_params(left='on', labelsize=15, labelcolor='k')
if (i+1) % 4 == 0:
cb.set_label(labels[i/4])
else:
plt.setp(cb.ax.get_yticklabels(), visible=False)
#cb.ax.set_axis_off()
for sink in sinklist[i]:
mscale = np.log(sink.mass) +1
ax.plot(sink.x, sink.y, 'ko', ms=mscale, mew=1)
ax.set_xlim(x.min(), x.max())
ax.set_ylim(y.min(), y.max())
if i > 7:
ax.text(0.5, 0.025, sim[i/3], color='w', ha='center', va='bottom', size=18,
transform=grid[i].transAxes, bbox=bbox_props)
plt.show()
fig.savefig('figures/diskprops.png', bbox_inches='tight', dpi=300)