diff --git a/Fig8_density_comparison.py b/Fig8_density_comparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..efd17b02331d24ab99070c209fbaa29f2f68d958
--- /dev/null
+++ b/Fig8_density_comparison.py
@@ -0,0 +1,174 @@
+import os
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+import numpy as np
+np.seterr(divide = 'ignore')
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'large',
+         'axes.labelsize': 'x-large',
+         'axes.titlesize':'x-large',
+         'figure.titlesize':'x-large',
+         'legend.title_fontsize':'large',
+         'legend.fontsize':'large',
+         'xtick.labelsize':'large',
+         'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+import sys
+sys.path.insert(0, "helpers/")
+
+from equations import density_critical_hmin
+
+from simulations_details import get_simulation_time_Myr
+from simulations_details import get_simulation_details_all
+
+from select_simulations_for_figures import Fig8_runnames as runnames
+from select_simulations_for_figures import Fig8_filebase_runs as filebase_runs
+from select_simulations_for_figures import Fig8_filebase_reruns as filebase_reruns
+
+from simulation_path import snapshotnumber
+
+from simulation_plotparameter import lognH_min
+from simulation_plotparameter import lognH_max
+from simulation_plotparameter import dlognH_ticks
+from simulation_plotparameter import logT_min
+from simulation_plotparameter import logT_max
+from simulation_plotparameter import massfrac_min
+from simulation_plotparameter import massfrac_max
+from simulation_plotparameter import nbins
+from simulation_plotparameter import colors
+from simulation_plotparameter import lognH_2Darr
+from simulation_plotparameter import logT_2Darr
+from simulation_plotparameter import nH_2Darr
+from simulation_plotparameter import T_2Darr
+from simulation_plotparameter import color_zone
+from simulation_plotparameter import eta_res
+
+from phasespace_plot import cumulative_mass_fraction
+from phasespace_plot import make_individual_phasespace_plot
+from phasespace_plot import add_zone_of_avoidance
+from phasespace_plot import add_Jeans_mass_contours
+
+snapshotnumber_rerun = 0
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig8_phasespace_comparison.png"
+
+
+# #################################################################
+# combines phase space plots and cumulative mass fraction plots
+# into one figure with:
+# row 1 (columns 0-2): phase-space plots for original snapshots for runs 
+#                      with 3 different minimum smoothing lengths
+# row 2 (columns 0-2): phase-space plots for real densities
+# final column: cumulative mass fractions for densities above x-axis value
+# #################################################################
+def make_plot():
+    mB_runs, eps_runs, h_min_runs, h_min_ratio_runs, XH_runs, esf_runs, \
+    mB_reruns, eps_reruns, h_min_reruns, h_min_ratio_reruns, XH_reruns, esf_reruns = \
+                                     get_simulation_details_all(filebase_runs, filebase_reruns, \
+                                                           snapshotnumber)
+
+    #get simulation time
+    time_Myr = get_simulation_time_Myr(filebase_runs[0] + '%4.4i'%(snapshotnumber) + '.hdf5')
+
+    fig = plt.figure(figsize = (12.3,5.0))
+    fig.subplots_adjust(bottom = 0.12 , left = 0.10, right = 0.97, top = 0.85)
+
+    titlestring  = "log m$_{\mathrm{B}}$ [M$_{\odot}$] = %.1f, "%(np.log10(mB_runs[0]))
+    titlestring += "$\epsilon$ = %i pc, "%(round(eps_runs[0]))
+    titlestring += "e$_{\mathrm{sf}}$ = %.2f"%(esf_runs[0]*100.)
+    titlestring += " per cent, "
+    titlestring += "t = %.2f Gyr"%(time_Myr/1000.)
+
+    fig.suptitle(titlestring)
+    gs = gridspec.GridSpec(2,5, width_ratios = [1.,1.,1.,0.32,1.], wspace = 0., hspace = 0.)
+
+    # add phasespace 2D histograms for original snapshot  
+    for iplot in range(len(filebase_runs)):
+        snapfilename = filebase_runs[iplot] + '%4.4i'%(snapshotnumber) + '.hdf5'
+        ax = plt.subplot(gs[iplot])
+        ax.set_title("h$_{\mathrm{min}}$ = %.2f pc"%(h_min_runs[iplot]))
+        ax.tick_params(axis = 'x', labelbottom = False, direction = 'inout')
+        mass_map = make_individual_phasespace_plot(ax, snapfilename)
+        add_zone_of_avoidance(ax, mB_runs[iplot], eps_runs[iplot], \
+                                  h_min_runs[iplot], XH_runs[iplot], mass_map)
+
+        add_Jeans_mass_contours(ax, mB_runs[iplot], eps_runs[iplot], \
+                                  h_min_runs[iplot], XH_runs[iplot])
+        if iplot == 0:
+            ax.set_ylabel('log T [K]')
+            ax.text(-0.34, 0.5, "SPH estimated\ngas densities", ha = 'center', va = 'center', transform = ax.transAxes, fontsize = 'x-large', rotation = 90)
+            ax.text(0.98, 0.5, "log M$_{\mathrm{C}}$", ha = 'right', va = 'baseline', transform = ax.transAxes, fontsize = 'x-large')
+        else:
+            ax.tick_params(axis = 'y', labelleft = False, direction = 'inout')
+
+        if iplot < 2:
+            ax.text(0.97, 0.95, '(C)', transform = ax.transAxes, ha = 'right', va = 'top', fontsize = 'x-large', zorder = 2000)
+        else:
+            ax.text(0.97, 0.05, '(C)', transform = ax.transAxes, ha = 'right', va = 'baseline', fontsize = 'x-large', zorder = 2000)
+        
+    # add phasespace 2D histograms for re-calculated (hmin -> 0) snapshot
+    for iplot in range(len(filebase_reruns)):
+        snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+        ax = plt.subplot(gs[iplot+5])
+        ax.set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+        mass_map = make_individual_phasespace_plot(ax, snapfilename)
+        add_zone_of_avoidance(ax, mB_runs[iplot], eps_runs[iplot], \
+                                  h_min_runs[iplot], XH_runs[iplot], mass_map)
+
+        add_Jeans_mass_contours(ax, mB_runs[iplot], eps_runs[iplot], \
+                                  h_min_runs[iplot], XH_runs[iplot])
+
+        if iplot == 0:
+            ax.set_ylabel('log T [K]')
+            ax.text(-0.34, 0.5, "Re-calculated\ngas densities", ha = 'center', va = 'center', transform = ax.transAxes, fontsize = 'x-large', rotation = 90)
+        else:
+            ax.tick_params(axis = 'y', labelleft = False, direction = 'inout')
+
+        if iplot < 2:
+            ax.text(0.97, 0.95, '(C)', transform = ax.transAxes, ha = 'right', va = 'top', fontsize = 'x-large', zorder = 2000)
+        else:
+            ax.text(0.97, 0.05, '(C)', transform = ax.transAxes, ha = 'right', va = 'baseline', fontsize = 'x-large', zorder = 2000)
+        
+    # mass fractions for >n 
+    ax_hi = plt.subplot(gs[4])
+    ax_hi.set_xlim(lognH_min, lognH_max)
+    ax_hi.xaxis.set_ticks(np.arange(lognH_min, lognH_max+dlognH_ticks, dlognH_ticks))
+    ax_hi.set_yscale('log')
+    ax_hi.set_ylim(massfrac_min, massfrac_max)
+    ax_hi.set_ylabel(r'M$_{ (> {\mathrm{n}}_{\mathrm{H}}) }$ / M$_{\mathrm{tot}}$')
+    ax_hi.tick_params(axis = 'x', labelbottom = False, direction = 'inout')
+    ax_lo = plt.subplot(gs[9])
+    ax_lo.set_xlim(lognH_min, lognH_max)
+    ax_lo.xaxis.set_ticks(np.arange(lognH_min, lognH_max+dlognH_ticks, dlognH_ticks))
+    ax_lo.set_yscale('log')
+    ax_lo.set_ylim(massfrac_min, massfrac_max)
+    ax_lo.set_ylabel(r'M$_{ (> {\mathrm{n}}_{\mathrm{H}}) }$ / M$_{\mathrm{tot}}$')
+    ax_lo.set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+
+    for iplot in range(len(filebase_runs)):
+        snapfilename = filebase_runs[iplot] + '%4.4i'%(snapshotnumber) + '.hdf5'
+        lognH, massfraction = cumulative_mass_fraction(snapfilename)       
+        ax_hi.plot(lognH, massfraction, color = colors[iplot], linestyle = 'solid', label = "%.2f"%(h_min_runs[iplot]))
+        nHcrit = density_critical_hmin(mB_runs[iplot], h_min_runs[iplot])
+        ax_hi.axvline(np.log10(nHcrit), ymin = 0.9, ymax = 1.0, color = colors[iplot], linestyle = 'dotted', lw = 2)
+        ax_lo.axvline(np.log10(nHcrit), ymin = 0.9, ymax = 1.0, color = colors[iplot], linestyle = 'dotted', lw = 2)
+
+    for iplot in range(len(filebase_reruns)):
+        snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+        lognH, massfraction = cumulative_mass_fraction(snapfilename)
+        ax_hi.plot(lognH, massfraction, color = colors[iplot], linestyle = 'dashed', lw = 1.0, alpha = 0.5)    
+        ax_lo.plot(lognH, massfraction, color = colors[iplot], linestyle = 'dashed')    
+
+    ax_hi.legend(loc = "lower left", title = 'h$_{\mathrm{min}}$ [pc]')
+
+    fig.savefig(outputname, dpi = 150)
+    plt.close()
+    print ("Figure saved: "+outputname)
+
+
+make_plot()
diff --git a/helpers/Fig6_runs_to_plot.py b/helpers/Fig6_runs_to_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..3faadab83927b8d28264ecd485793cd392f58334
--- /dev/null
+++ b/helpers/Fig6_runs_to_plot.py
@@ -0,0 +1,44 @@
+import glob
+from pathlib import Path
+import os
+
+########################################
+snapshotnumber = 100
+########################################
+snapshotnumber_rerun = 0
+folderbase = "/gpfs/data/fs71897/ploeckinger/SIMULATIONS/2022_11_IsolatedGalaxies_v2/"
+########################################
+# Simulation set run variations
+########################################
+mass_resolution_levels = ["M5"]
+softening_resolutions_in_kpc = [0.25]
+smoothing_resolutions_in_kpc = [0.0775, 0.0245, 0.00775]
+star_formation_efficiencies = [0.00316]
+########################################
+
+########################################
+# Set folder names and run names
+########################################
+runfolders = []
+rerunfolders = []
+foffolders = []
+for mass in mass_resolution_levels:
+    for soft in softening_resolutions_in_kpc:
+        for smooth in smoothing_resolutions_in_kpc:
+            for sfe in star_formation_efficiencies:
+                ##########################################
+                runname = "Galaxy" + mass + \
+                          "_soft%4.4ipc"%(int(soft*1000)) + \
+                          "_hmin%06.2fpc"%(smooth*1000.) + \
+                          "_sfe%06.3f"%(sfe)
+
+                runfolders.append( os.path.join(folderbase, mass, runname) )
+                rerunfolders.append( os.path.join(folderbase, mass+"_reruns_snap%4.4i"%(snapshotnumber), "Rerun" + runname) )
+                foffolders.append( os.path.join(folderbase, mass+"_fofruns_snap%4.4i"%(snapshotnumber), "FOF"+runname)  )
+
+filebase_runs   = [s + "/output_" for s in runfolders]
+filebase_reruns = [s + "/output_" for s in rerunfolders]
+filebase_fofruns = [s + "/fof_output_0000.hdf5" for s in foffolders]
+
+runnames = [os.path.split(s)[1] for s in runfolders]
+
diff --git a/helpers/phasespace_plot.py b/helpers/phasespace_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..67a79950065da8e2da407d617df3f2e5e29375ed
--- /dev/null
+++ b/helpers/phasespace_plot.py
@@ -0,0 +1,149 @@
+import os
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+import numpy as np
+np.seterr(divide = 'ignore') 
+
+from swiftsimio import load
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'large',
+         'axes.labelsize': 'x-large',
+         'axes.titlesize':'x-large',
+         'figure.titlesize':'x-large',
+         'legend.title_fontsize':'large',
+         'legend.fontsize':'large',
+         'xtick.labelsize':'large',
+         'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+from unyt import proton_mass_cgs as mH
+from unyt import cm, K, g
+from unyt import msun, pc
+
+from simulation_plotparameter import lognH_min
+from simulation_plotparameter import lognH_max
+from simulation_plotparameter import dlognH_ticks
+from simulation_plotparameter import logT_min
+from simulation_plotparameter import logT_max
+from simulation_plotparameter import massfrac_min
+from simulation_plotparameter import massfrac_max
+from simulation_plotparameter import nbins
+from simulation_plotparameter import colors
+from simulation_plotparameter import lognH_2Darr
+from simulation_plotparameter import logT_2Darr
+from simulation_plotparameter import nH_2Darr
+from simulation_plotparameter import T_2Darr
+from simulation_plotparameter import color_zone
+from simulation_plotparameter import color_zone_C
+from simulation_plotparameter import color_zone_D
+
+from equations import Jeans_length_0
+from equations import Jeans_length_W
+from equations import Jeans_mass_0
+from equations import Jeans_mass_W
+from equations import density_critical_hmin
+from equations import eta_res
+from equations import h_smooth
+
+
+Nneigh = 50
+
+# #################################################################
+# adds zone of induced collapse to axis
+# #################################################################
+def add_zone_of_induced_collapse(ax, mB, eps, XH):
+
+    softening_scale = 1.5 * eps
+    kernel_size = 1.936492 * h_smooth(mB, nH_2Darr, 0.)
+
+    all_zones = np.zeros_like(nH_2Darr)
+    all_zones[...] = 1.
+
+    # Zone D:
+    all_zones[ (Jeans_length_W(nH_2Darr, T_2Darr, eps) >  softening_scale) & \
+               (Jeans_length_W(nH_2Darr, T_2Darr, eps) <= kernel_size) ] = 0.
+
+    ax.contourf(lognH_2Darr, logT_2Darr, all_zones, levels = np.asarray([0., 0.99]), colors = color_zone_D, alpha = 0.8)
+    ax.contour (lognH_2Darr, logT_2Darr, all_zones, levels = np.asarray([0.99]), colors = color_zone_D, zorder = 100)    
+
+    return
+    
+
+# #################################################################
+# adds zone of avoidance to axis
+# #################################################################
+def add_zone_of_avoidance(ax, mB, eps, h_min, XH, mass_map):
+    zone_of_avoidance = np.zeros_like(lognH_2Darr)
+    zone_of_avoidance[:,:] = 1.
+    zone_of_avoidance[(Jeans_length_W(nH_2Darr, T_2Darr, eps) / h_min < 1.) & \
+                      (lognH_2Darr > np.log10(density_critical_hmin(mB, h_min)))] = 0.
+
+    mass_zone = np.sum(mass_map[zone_of_avoidance == 0.])
+
+    ax.contourf(lognH_2Darr, logT_2Darr, zone_of_avoidance, levels = np.asarray([0., 0.99]), colors = color_zone_C, alpha = 0.8)
+    ax.contour (lognH_2Darr, logT_2Darr, zone_of_avoidance, levels = np.asarray([0.99]), colors = color_zone_C, zorder = 100)
+    if mass_zone > 0.:
+        ax.text(0.98, 0.38, "%.2f"%(np.log10(mass_zone)), ha = 'right', va = 'baseline', transform = ax.transAxes, fontsize = 'x-large')
+
+    return
+
+# #################################################################
+# adds contours for softened Jeans mass equals 1 and 10 times 
+# the baryon particle mass
+# #################################################################
+def add_Jeans_mass_contours(ax, mB, eps, h_min, XH):
+    logMJ_W_2D_over_mB = np.log10(Jeans_mass_W(nH_2Darr, T_2Darr, eps) / (1. * mB))
+
+    ax.contour (lognH_2Darr, logT_2Darr, logMJ_W_2D_over_mB, levels = np.asarray([2.]), colors = 'black', linestyles = 'solid', zorder = 100, linewidths = 1)
+    ax.contour (lognH_2Darr, logT_2Darr, logMJ_W_2D_over_mB, levels = np.asarray([1.]), colors = 'black', linestyles = 'dashed', zorder = 100, linewidths = 1)
+    ax.contour (lognH_2Darr, logT_2Darr, logMJ_W_2D_over_mB, levels = np.asarray([0.]), colors = 'black', linestyles = 'dotted', zorder = 100, linewidths = 1)
+
+    return
+
+# #################################################################
+# plots 2D density temperature histogram into axis
+# for one snapshot
+# #################################################################
+def make_individual_phasespace_plot(ax, filename):
+    data = load(filename)
+    SPH_hydrogen_fraction = data.gas.element_mass_fractions.hydrogen
+    SPH_number_density    = (data.gas.densities / mH * SPH_hydrogen_fraction ).to(cm ** -3)
+    SPH_temperatures      = (data.gas.temperatures).to('K')
+    SPH_masses            = (data.gas.masses).to('msun')
+
+    SPH_hist, xedge, yedge = np.histogram2d(np.log10(SPH_number_density), np.log10(SPH_temperatures), bins = nbins, \
+                         range = [[lognH_min, lognH_max], [logT_min, logT_max]], weights = SPH_masses, density = False)  
+
+    SPH_hist = SPH_hist.T
+    logSPH_hist = np.log10(SPH_hist)
+
+    ax.imshow(logSPH_hist, interpolation='none', origin='lower',\
+              extent=[xedge[0], xedge[-1], yedge[0], yedge[-1]], \
+              aspect='auto', zorder = 10, vmin = 5., vmax = 7., cmap = 'viridis')
+
+    ax.set_xlim(lognH_min, lognH_max)
+    ax.set_ylim(logT_min, logT_max)
+    ax.xaxis.set_ticks(np.arange(lognH_min, lognH_max, dlognH_ticks))
+
+    return SPH_hist
+
+# #################################################################
+# calculates the mass fraction of particles with densities above n
+# for various n
+# #################################################################
+def cumulative_mass_fraction(filename):
+    data = load(filename)
+    data.gas.masses.convert_to_units('Msun')
+    data.gas.densities.convert_to_units(g * cm**-3)
+    nH = data.gas.densities / mH * data.gas.element_mass_fractions.hydrogen
+    nH.convert_to_units(cm**-3)
+
+    data.gas.masses[np.log10(nH.value) < lognH_min] = 0.    
+
+    histbins_dens = np.logspace(lognH_min, lognH_max, num = nbins, endpoint=True, base=10.0)
+    hist, edges = np.histogram(nH, bins = histbins_dens, weights = data.gas.masses.value, density = False)
+    cumhist  = np.cumsum(hist) / np.sum(data.gas.masses.value)
+
+    return np.log10(histbins_dens[:-1]), 1. - cumhist
+
diff --git a/helpers/select_simulations_for_figures.py b/helpers/select_simulations_for_figures.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ea8c5415af41c476a711ffc9526d5270e4afe95
--- /dev/null
+++ b/helpers/select_simulations_for_figures.py
@@ -0,0 +1,18 @@
+#################################################
+# Ploeckinger et al. (2023, subm)
+#################################################
+
+from simulations_details import get_runfolders
+
+#################################################
+# Simulations to plot for Figure 8
+# 3 different simulations; figure captions assume
+# identical values for epsilon, e_sf, and mB and
+# different values for h_min
+#################################################
+Fig8_mB   = ["M5", "M5", "M5"]   
+Fig8_eps  = [0.25, 0.25, 0.25]    # kpc
+Fig8_hmin = [0.0775, 0.0245, 0.00775]  # kpc
+Fig8_sfe  = [0.00316, 0.00316, 0.00316]
+
+Fig8_runnames, Fig8_filebase_runs, Fig8_filebase_reruns, Fig8_filebase_fofruns = get_runfolders(Fig8_mB, Fig8_eps, Fig8_hmin, Fig8_sfe)
diff --git a/helpers/simulation_path.py b/helpers/simulation_path.py
new file mode 100644
index 0000000000000000000000000000000000000000..3b0cddb57b77b154185a7badc557b35dc4d3cce7
--- /dev/null
+++ b/helpers/simulation_path.py
@@ -0,0 +1,12 @@
+
+#################################################
+# Set the path to the simulation directory
+#################################################
+simulation_path = "/gpfs/data/fs71897/ploeckinger/SIMULATIONS/2022_11_IsolatedGalaxies_v2/"
+
+#################################################
+# Select the snapshot number of original simulation
+# from which the simulation is restarted
+# (default: 100)
+#################################################
+snapshotnumber = 100
diff --git a/helpers/simulations_details.py b/helpers/simulations_details.py
new file mode 100644
index 0000000000000000000000000000000000000000..e2de310fa7fe61b8f9a668c63a0f0e4af1ca7bbf
--- /dev/null
+++ b/helpers/simulations_details.py
@@ -0,0 +1,119 @@
+import os
+import sys
+from swiftsimio import load
+import numpy as np
+
+from simulation_path import simulation_path as folderbase
+from simulation_path import snapshotnumber
+
+####################################################
+# Returns the full paths and runnames for a simulation with the
+# selected values for mB, eps, hmin, and the star
+# formation efficiency
+####################################################
+def get_runfolders(mB, eps, hmin, starforme):
+    if len(mB) == len(eps) == len(hmin) == len(starforme):
+
+        runfolders = []
+        rerunfolders = []
+        foffolders = []
+        for irun in range (len(mB)):
+            mass = mB[irun]
+            soft = eps[irun]
+            smooth = hmin[irun] 
+            sfe    = starforme[irun]
+            runname = "Galaxy" + mass + \
+                          "_soft%4.4ipc"%(int(soft*1000)) + \
+                          "_hmin%06.2fpc"%(smooth*1000.) + \
+                          "_sfe%06.3f"%(sfe)
+            runfolders.append( os.path.join(folderbase, mass, runname) )
+            rerunfolders.append( os.path.join(folderbase, mass+"_reruns_snap%4.4i"%(snapshotnumber), "Rerun" + runname) )
+            foffolders.append( os.path.join(folderbase, mass+"_fofruns_snap%4.4i"%(snapshotnumber), "FOF"+runname)  )
+
+        filebase_runs   = [s + "/output_" for s in runfolders]
+        filebase_reruns = [s + "/output_" for s in rerunfolders]
+        filebase_fofruns = [s + "/fof_output_0000.hdf5" for s in foffolders]
+
+        runnames = [os.path.split(s)[1] for s in runfolders]
+
+    else:
+        print("[get_runfolders] inconsistent simulation parameters")
+        sys.exit()
+
+    return runnames, filebase_runs, filebase_reruns, filebase_fofruns
+
+# #################################################################
+# returns the snapshot time in Myr
+# #################################################################
+def get_simulation_time_Myr(filename):
+    data = load(filename)
+    data.metadata.time.convert_to_units('Myr')
+    return data.metadata.time.value
+
+
+# #################################################################
+# returns some basic simulation information from one snapshot
+# #################################################################
+def get_simulation_details(filename):
+    data = load(filename)
+    eps         = float(data.metadata.parameters.get('Gravity:max_physical_baryon_softening')) * 1000.
+    h_min_ratio = float(data.metadata.parameters.get('SPH:h_min_ratio'))
+    h_min       = float(data.metadata.parameters.get('SPH:h_min_ratio')) * 1.55 * eps
+    XH_init     = float(data.metadata.parameters.get('SPH:H_mass_fraction'))
+    e_sf        = float(data.metadata.parameters.get('EAGLEStarFormation:star_formation_efficiency'))
+
+    Mgas_i = data.metadata.initial_mass_table.gas.to('Msun').value
+    if Mgas_i == 0:
+        mB = data.gas.masses.to('Msun')
+        Mgas_i = (np.sum(mB)/float(len(mB))).value
+    mB = Mgas_i
+
+    return mB, eps, h_min, h_min_ratio, XH_init, e_sf
+
+# #################################################################
+# returns some basic simulation information from list of runs and reruns
+# #################################################################
+def get_simulation_details_all(filebase_runs, filebase_reruns, \
+                                snapshotnumber):
+
+    snapshotnumber_rerun = 0
+    mB_runs  = []
+    eps_runs = []
+    h_min_runs = []
+    h_min_ratio_runs = []
+    XH_runs = []
+    e_sf_runs = []
+
+    for iplot in range(len(filebase_runs)):
+        snapfilename = filebase_runs[iplot] + '%4.4i'%(snapshotnumber) + '.hdf5'
+        mB, eps, h_min, h_min_ratio, XH_init, e_sf = get_simulation_details(snapfilename)
+        mB_runs.append(mB)
+        eps_runs.append(eps)
+        h_min_runs.append(h_min)
+        h_min_ratio_runs.append(h_min_ratio)
+        XH_runs.append(XH_init)
+        e_sf_runs.append(e_sf)
+
+    mB_reruns  = []
+    eps_reruns = []
+    h_min_reruns = []
+    h_min_ratio_reruns = []
+    XH_reruns = []
+    e_sf_reruns = []
+
+    for iplot in range(len(filebase_reruns)):
+        snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+        mB, eps, h_min, h_min_ratio, XH_init, e_sf = get_simulation_details(snapfilename)
+        mB_reruns.append(mB)
+        eps_reruns.append(eps)
+        h_min_reruns.append(h_min)
+        h_min_ratio_reruns.append(h_min_ratio)
+        XH_reruns.append(XH_init)
+        e_sf_reruns.append(e_sf)
+
+    return np.asarray(mB_runs), np.asarray(eps_runs), np.asarray(h_min_runs), \
+           np.asarray(h_min_ratio_runs), np.asarray(XH_runs), np.asarray(e_sf_runs),  \
+           np.asarray(mB_reruns), np.asarray(eps_reruns), np.asarray(h_min_reruns), \
+           np.asarray(h_min_ratio_reruns), np.asarray(XH_reruns), np.asarray(e_sf_reruns)
+
+