diff --git a/Fig10_alternative.py b/Fig10_alternative.py
new file mode 100644
index 0000000000000000000000000000000000000000..9c2a77c94cb47090e767bdbec5d2d15a1cdfa027
--- /dev/null
+++ b/Fig10_alternative.py
@@ -0,0 +1,168 @@
+import os
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib import patches
+
+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':'xx-large',
+ 'ytick.labelsize':'xx-large'}
+plt.rcParams.update(params)
+
+import sys
+sys.path.insert(0, "helpers/")
+
+from phasespace_plot import make_individual_phasespace_plot
+from phasespace_plot import add_zone_of_avoidance
+from phasespace_plot import add_Jeans_mass_contours
+from phasespace_plot import add_zone_of_induced_collapse
+
+from stellar_surface_density_plot import add_stellar_surface_density
+from stellar_surface_density_plot import add_fof_groups
+
+from simulations_details import get_simulation_time_Myr
+from simulations_details import get_simulation_details_all
+
+from select_simulations_for_figures import Fig10a_runnames as runnames_M5
+from select_simulations_for_figures import Fig10a_filebase_runs as filebase_runs_M5
+from select_simulations_for_figures import Fig10a_filebase_reruns as filebase_reruns_M5
+from select_simulations_for_figures import Fig10a_filebase_fofruns as filebase_fofruns_M5
+from simulation_path import snapshotnumber as snapshotnumber_M5
+
+from select_simulations_for_figures import Fig10b_runnames as runnames_M6
+from select_simulations_for_figures import Fig10b_filebase_runs as filebase_runs_M6
+from select_simulations_for_figures import Fig10b_filebase_reruns as filebase_reruns_M6
+from select_simulations_for_figures import Fig10b_filebase_fofruns as filebase_fofruns_M6
+from simulation_path import snapshotnumber as snapshotnumber_M6
+
+snapshotnumber_rerun = 0
+############################################################################
+# Set output filename
+############################################################################
+outputname_M5 = "Fig10a_gallery_alternative.png"
+outputname_M6 = "Fig10b_gallery_alternative.png"
+
+############################################################################
+# Some extra plotting settings
+############################################################################
+
+axis_dx_left = 0.1
+axis_dy_bottom = 0.1
+axis_xsize = 0.44
+axis_ysize = 0.3
+
+figsizex = 8.
+figsizey = 8.
+
+# #################################################################
+# Creates phasespace and stellar mass surface density plots
+# of 2 runs identified in helpers/Fig8_runs_to_plot.py
+# and of 2 runs identified in helpers/Fig9_runs_to_plot.py
+# #################################################################
+def make_gallery(filebase_runs, filebase_reruns, filebase_fofruns, snapshotnumber, outputname):
+ 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 = (figsizex,figsizex))
+ if (round(mB_runs[0] / 1.e5) == 1):
+ suptitlestring = "Mass resolution: 10$^{5}$ M$_{\odot}$, star formation efficiency: %.2f"%(100.*esf_runs[0])
+ else:
+ suptitlestring = r"Mass resolution: %i $\times$ 10$^{5}$ M$_{\odot}$, star formation efficiency: %.2f"%(round(mB_runs[0]/1.e5),100.*esf_runs[0])
+ suptitlestring += " per cent"
+ fig.suptitle(suptitlestring, fontsize = 16)
+
+ #######################################################
+ # Setup axes
+ #######################################################
+ ax = []
+ # 0: bottom left panel
+ ax.append(fig.add_axes([axis_dx_left, \
+ axis_dy_bottom, \
+ axis_xsize, \
+ axis_ysize]))
+
+ # 1: bottom right panel
+ ax.append(fig.add_axes([axis_dx_left + axis_xsize, \
+ axis_dy_bottom, \
+ axis_xsize, \
+ axis_ysize]))
+
+ # 2: top left panel
+ ax.append(fig.add_axes([axis_dx_left, \
+ axis_dy_bottom + axis_ysize, \
+ axis_xsize, \
+ axis_xsize]))
+
+ # 3: top right panel
+ ax.append(fig.add_axes([axis_dx_left + axis_xsize, \
+ axis_dy_bottom + axis_ysize, \
+ axis_xsize, \
+ axis_xsize]))
+
+ for iplot in range(len(filebase_reruns)):
+ snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+ mass_map = make_individual_phasespace_plot(ax[iplot], snapfilename)
+ add_zone_of_avoidance(ax[iplot], mB_runs[iplot], eps_runs[iplot], \
+ h_min_runs[iplot], XH_runs[iplot], mass_map)
+ add_Jeans_mass_contours(ax[iplot], mB_runs[iplot], eps_runs[iplot], \
+ h_min_runs[iplot], XH_runs[iplot])
+
+ ax[iplot].xaxis.set_tick_params(labelsize=16)
+ ax[iplot].yaxis.set_tick_params(labelsize=16)
+ ax[iplot].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]', fontsize = 16)
+
+ if iplot > 0:
+ ax[iplot].set_yticklabels([])
+ else:
+ ax[iplot].set_ylabel('log T [K]', fontsize = 16)
+
+
+ for iplot in range(len(filebase_reruns)):
+ snapfilename = filebase_reruns[iplot] + '%4.4i'%(snapshotnumber_rerun) + '.hdf5'
+ mass_map_plot = add_stellar_surface_density(ax[iplot+2], snapfilename)
+
+ foffilename = filebase_fofruns[iplot]
+ add_fof_groups(ax[iplot+2], foffilename)
+ ax[iplot+2].set_axis_on()
+
+
+ if iplot == 0:
+ titlestring = "Fiducial:\n$\epsilon$ = %i pc, h$_{\mathrm{min}}$ = %.2f pc"%(eps_runs[iplot],\
+ h_min_runs[iplot])
+ ax[iplot+2].set_title(titlestring, fontsize = 16, loc = 'left')
+ else:
+ titlestring = "Alternative:\n$\epsilon$ = %i pc, h$_{\mathrm{min}}$ = %.2f pc"%(eps_runs[iplot],\
+ h_min_runs[iplot])
+ ax[iplot+2].set_title(titlestring, fontsize = 16, loc = 'right')
+
+ ax[0].text(0.98, 0.5, "log M$_{\mathrm{zone}}$", ha = 'right', va = 'baseline', \
+ transform = ax[0].transAxes, fontsize = 'x-large')
+
+ #ax[0].text(0.97, 0.95, '(C)', transform = ax[0].transAxes, ha = 'right', va = 'top', fontsize = 'x-large', zorder = 2000)
+ #ax[0].text(0.01, 0.15, '(D)', transform = ax[0].transAxes, ha = 'left', va = 'bottom', fontsize = 'x-large', zorder = 2000)
+ #ax[1].text(0.01, 0.15, '(D)', transform = ax[1].transAxes, ha = 'left', va = 'bottom', fontsize = 'x-large', zorder = 2000)
+
+
+ plt.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+ return
+
+make_gallery(filebase_runs_M5, filebase_reruns_M5, filebase_fofruns_M5, snapshotnumber_M5, outputname_M5)
+make_gallery(filebase_runs_M6, filebase_reruns_M6, filebase_fofruns_M6, snapshotnumber_M6, outputname_M6)
diff --git a/helpers/select_simulations_for_figures.py b/helpers/select_simulations_for_figures.py
index 913b365afb56ed41a0482ae55effc7002c700b8c..14fd42d702c8af680781159069d0aa52a7997a16 100644
--- a/helpers/select_simulations_for_figures.py
+++ b/helpers/select_simulations_for_figures.py
@@ -31,3 +31,23 @@ Fig9_sfe = [0.01, 0.00316, 0.00316, 0.00316, 0.00316]
Fig9_runnames, Fig9_filebase_runs, Fig9_filebase_reruns, Fig9_filebase_fofruns = get_runfolders(Fig9_mB, Fig9_eps, Fig9_hmin, Fig9_sfe)
+#################################################
+# Simulations to plot for Figure 10
+# 4 different simulations in total
+# 2 for each mass resolution level (10a, 10b)
+#################################################
+
+Fig10a_mB = ["M5", "M5"]
+Fig10a_eps = [0.25, 0.5]
+Fig10a_hmin = [0.0775, 0.00775]
+Fig10a_sfe = [0.00316, 0.00316]
+
+Fig10a_runnames, Fig10a_filebase_runs, Fig10a_filebase_reruns, Fig10a_filebase_fofruns = get_runfolders(Fig10a_mB, Fig10a_eps, Fig10a_hmin, Fig10a_sfe)
+
+Fig10b_mB = ["M6", "M6"]
+Fig10b_eps = [0.25, 0.5]
+Fig10b_hmin = [0.0775, 0.00775]
+Fig10b_sfe = [0.00316, 0.00316]
+
+Fig10b_runnames, Fig10b_filebase_runs, Fig10b_filebase_reruns, Fig10b_filebase_fofruns = get_runfolders(Fig10b_mB, Fig10b_eps, Fig10b_hmin, Fig10b_sfe)
+