diff --git a/examples/ZeldovichPancake_3D/plotSolution.py b/examples/ZeldovichPancake_3D/plotSolution.py
index 599283d683be8f7ca7722523fd6b8fc571487ecd..d9e78b68fc39f99d030ecca058367021365d31c2 100644
--- a/examples/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/ZeldovichPancake_3D/plotSolution.py
@@ -33,6 +33,7 @@ import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
+import os.path
# Plot parameters
params = {'axes.labelsize': 10,
@@ -85,28 +86,48 @@ phi = sim["/PartType0/Potential"][:]
x -= 0.5 * boxSize
-sim_g = h5py.File("snapshot_%03d.hdf5"%snap, "r")
-x_g = sim_g["/PartType0/Coordinates"][:,0]
-v_g = sim_g["/PartType0/Velocities"][:,0]
-u_g = sim_g["/PartType0/InternalEnergy"][:]
-rho_g = sim_g["/PartType0/Density"][:]
-phi_g = sim_g["/PartType0/Potential"][:]
-
-x_g -= 0.5 * boxSize
+filename_g = "snapshot_%03d.hdf5"%snap
+print filename_g
+if os.path.exists(filename_g):
+ print "exists!"
+ sim_g = h5py.File(filename_g, "r")
+ x_g = sim_g["/PartType0/Coordinates"][:,0]
+ v_g = sim_g["/PartType0/Velocities"][:,0]
+ u_g = sim_g["/PartType0/InternalEnergy"][:]
+ rho_g = sim_g["/PartType0/Density"][:]
+ phi_g = sim_g["/PartType0/Potential"][:]
+
+ x_g -= 0.5 * boxSize
+else:
+ x_g = np.zeros(1)
+ v_g = np.zeros(1)
+ u_g = np.zeros(1)
+ rho_g = np.zeros(1)
+ phi_g = np.zeros(1)
# Derived parameters
rho_0 = m.sum() / (boxSize**3) # critical density of the box
lambda_i = boxSize # wavelength of the perturbation
-x_s = linspace(-0.5 * lambda_i, 0.5 * lambda_i, 1000)
+#x_s = linspace(-0.5 * lambda_i, 0.5 * lambda_i, 256)
+#k_i = 2. * pi / lambda_i
+#zfac = (1. + z_c) / (1. + redshift)
+#rho_s = rho_0 / (1 - zfac * cos(k_i * x_s))
+#v_s = -H_0 * (1. + z_c) / sqrt(1. + redshift) * sin(k_i * x_s) / k_i
+#T_s = T_i * (((1. + redshift) / (1. + z_i))**3 / (1. - zfac * cos(k_i * x_s)))**(2. / 3.)
+#P_s = zeros(256)
+#u_s = zeros(256)
+#s_s = zeros(256)
+
+# Solution taken from Springel 2010. Eqs. 127 - 130
+q = linspace(-0.5 * lambda_i, 0.5 * lambda_i, 256)
k_i = 2. * pi / lambda_i
zfac = (1. + z_c) / (1. + redshift)
-rho_s = rho_0 / (1 - zfac * cos(k_i * x_s))
-v_s = -H_0 * (1. + z_c) / sqrt(1. + redshift) * sin(k_i * x_s) / k_i
-T_s = T_i * (((1. + redshift) / (1. + z_i))**3 / (1. - zfac * cos(k_i * x_s)))**(2. / 3.)
-P_s = zeros(1000)
-u_s = zeros(1000)
-s_s = zeros(1000)
+x_s = q - zfac * sin(k_i * q) / k_i
+rho_s = rho_0 / (1. - zfac * cos(k_i * q))
+v_s = -H_0 * (1. + z_c) / sqrt(1. + redshift) * sin(k_i * q) / k_i
+T_s = T_i * (((1. + redshift) / ( 1. + z_i))**3. * rho_s / rho_0)**(2. / 3.)
+
unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
@@ -121,7 +142,8 @@ figure()
# Velocity profile --------------------------------
subplot(231)
-plot(x_g, v_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
+if np.size(x_g) > 1:
+ plot(x_g, v_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
plot(x, v, '.', color='r', ms=4.0)
plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
@@ -130,7 +152,8 @@ ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
# Density profile --------------------------------
subplot(232)
-plot(x_g, rho_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
+if np.size(x_g) > 1:
+ plot(x_g, rho_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
plot(x, rho, '.', color='r', ms=4.0)
plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
@@ -138,16 +161,14 @@ ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
# Potential profile --------------------------------
subplot(233)
-plot(x_g, phi_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
+if np.size(x_g) > 1:
+ plot(x_g, phi_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
plot(x, phi, '.', color='r', ms=4.0)
-#plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Potential}}~\\phi$", labelpad=0)
# Internal energy profile -------------------------
-subplot(234)
-#plot(x, u, '.', color='r', ms=4.0)
-#plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+subplot(234, yscale="log")
u *= (unit_length_in_si**2 / unit_time_in_si**2)
u_g *= (unit_length_in_si**2 / unit_time_in_si**2)
u /= a**(3 * (gas_gamma - 1.))
@@ -155,7 +176,8 @@ u_g /= a**(3 * (gas_gamma - 1.))
T = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
T_g = (gas_gamma - 1.) * u_g * mH_in_kg / k_in_J_K
print "z = {0:.2f}, T_avg = {1:.2f}".format(redshift, T.mean())
-plot(x_g, T_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
+if np.size(x_g) > 1:
+ plot(x_g, T_g, 's', color='g', alpha=0.8, lw=1.2, ms=4)
plot(x, T, '.', color='r', ms=4.0)
plot(x_s, T_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
@@ -176,4 +198,4 @@ ylim(0, 1)
xticks([])
yticks([])
-savefig("ZeldovichPancake.png", dpi=200)
+savefig("ZeldovichPancake_%.3d.png"%snap, dpi=200)
diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index aadc982e2fb4463194b277db3702518b0b66eb7b..ec20a7121fc7c10212a0b7a8c4a73a3e908a6c12 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -15,7 +15,7 @@ TimeIntegration:
Snapshots:
basename: zeldovichPancake # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 1.01 # Time difference between consecutive outputs (in internal units)
+ delta_time: 1.04 # Time difference between consecutive outputs (in internal units)
scale_factor_first: 0.00991
# Parameters governing the conserved quantities statistics