diff --git a/Fig9b_gallery_stars.py b/Fig9b_gallery_stars.py
new file mode 100644
index 0000000000000000000000000000000000000000..f48337ac8d73a74440fedf1eeb6f43e1ae68ab2a
--- /dev/null
+++ b/Fig9b_gallery_stars.py
@@ -0,0 +1,133 @@
+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':'large',
+ 'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+import sys
+sys.path.insert(0, "helpers/")
+
+from plotannotation import add_arrows
+
+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 Fig9_runnames as runnames
+from select_simulations_for_figures import Fig9_filebase_runs as filebase_runs
+from select_simulations_for_figures import Fig9_filebase_reruns as filebase_reruns
+from select_simulations_for_figures import Fig9_filebase_fofruns as filebase_fofruns
+
+from simulation_path import snapshotnumber
+
+
+snapshotnumber_rerun = 0
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig7_gallery_stars.png"
+
+############################################################################
+# Some extra plotting settings
+############################################################################
+
+axis_xsize = 0.40
+axis_ysize = axis_xsize
+axis_dy = 0.015
+axis_dx_left = 0.02
+axis_dx = 0.2
+axis_dy_bottom = 0.02
+axis_overlap = 0.1
+axis_overlap_x = 0.02
+
+# #################################################################
+# Creates 5 stellar mass density plots
+# of the 5 runs identified in helpers/Fig7_runs_to_plot.py
+# Labels need to be updated if different runs are selected
+# #################################################################
+
+def make_stars_gallery_plot(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 = (7.5*0.9,7*0.9))
+ axframe = fig.add_axes([0,0,1,1])
+ for axis in ['top','bottom','left','right']:
+ axframe.spines[axis].set_linewidth(10)
+ axframe.spines[axis].set_color("black")
+ axframe.spines[axis].set_zorder(0)
+
+ ax = []
+
+ #######################################################
+ # Setup axes
+ #######################################################
+ # 0: top left panel
+ ax.append(fig.add_axes([axis_dx_left, \
+ 1. - axis_dy_bottom - axis_ysize, \
+ axis_xsize, \
+ axis_ysize]))
+ # 1: top right panel
+ ax.append(fig.add_axes([1. - axis_dx_left - axis_xsize , \
+ 1. - axis_dy_bottom - axis_ysize, \
+ axis_xsize, \
+ axis_ysize]))
+ # 2: middle panel
+ ax.append(fig.add_axes([0.5 - 0.5 * axis_xsize, \
+ 0.5 - 0.5 * axis_ysize, \
+ axis_xsize, \
+ axis_ysize]))
+ # 3: bottom left panel
+ ax.append(fig.add_axes([axis_dx_left, \
+ axis_dy_bottom, \
+ axis_xsize, \
+ axis_ysize]))
+ # 4: bottom right panel
+ ax.append(fig.add_axes([1. - axis_dx_left - axis_xsize , \
+ axis_dy_bottom, \
+ axis_xsize, \
+ axis_ysize]))
+
+ 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], snapfilename)
+
+ foffilename = filebase_fofruns[iplot]
+ add_fof_groups(ax[iplot], foffilename)
+
+ ax[2].text(1.2, 0.5, "Increase m$_{\mathrm{B}}$\n (x 8)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'center', ha = 'left')
+ ax[2].text(0.5, 1.35, "Increase $\epsilon$\n (x 2)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'bottom', ha = 'center')
+ ax[2].text(0.5,-0.35, "Decrease h$_{\mathrm{min}}$\n (x 0.1)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'top', ha = 'center')
+ ax[2].text(-0.2,0.5, "Increase e$_{\mathrm{sf}}$\n (x 3)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'center', ha = 'right')
+
+ add_arrows(fig, ax)
+
+ plt.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+ return
+
+make_stars_gallery_plot(outputname)
diff --git a/helpers/stellar_surface_density_plot.py b/helpers/stellar_surface_density_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..9b39e1bc2f9e1b180e843e3c65219bd6d1a7e028
--- /dev/null
+++ b/helpers/stellar_surface_density_plot.py
@@ -0,0 +1,174 @@
+import os
+import numpy as np
+np.seterr(divide = 'ignore')
+
+import h5py
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.pyplot as plt
+from matplotlib import gridspec
+from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
+
+from unyt import unyt_array, unyt_quantity
+from unyt import msun, kpc, g, cm, K, pc
+from unyt import boltzmann_constant_cgs as kB
+from unyt import proton_mass_cgs as mH
+
+import swiftsimio as sw
+from swiftsimio import load
+from swiftsimio.visualisation.projection import project_pixel_grid
+from swiftsimio.visualisation.slice import kernel_gamma
+from swiftsimio.visualisation.smoothing_length_generation import (
+ generate_smoothing_lengths,
+)
+# ignore cosmology warnings
+import warnings
+warnings.filterwarnings('ignore', module='swiftsimio')
+
+from simulations_details import get_simulation_time_Myr
+from simulations_details import get_simulation_details_all
+
+from simulation_plotparameter import MaxAgeYoungStarsMyr
+from simulation_plotparameter import r_img_kpc
+from simulation_plotparameter import slice_thickness_in_kpc
+from simulation_plotparameter import npix
+from simulation_plotparameter import nbins
+from simulation_plotparameter import vmin
+from simulation_plotparameter import vmax
+from simulation_plotparameter import vmax_sigma_line
+from simulation_plotparameter import cmap_stars
+from simulation_plotparameter import colors
+from simulation_plotparameter import sigma_min
+from simulation_plotparameter import eta_res
+from simulation_plotparameter import mask_central_kpc
+
+# #################################################################
+# Calculates the fraction of pixels in a mass map above a
+# stellar surface density Sigma for various Sigma
+# #################################################################
+def get_cumulative_pixel_fraction(mass_map):
+
+ surface_densities = np.power(10., mass_map.ravel())
+ surface_densities[surface_densities < sigma_min] = 0.
+ indx = np.where(surface_densities >= sigma_min)[0]
+ histbins_sigma = np.logspace(vmin, vmax_sigma_line, num = nbins, endpoint=True, base=10.0)
+ hist, edges = np.histogram(surface_densities, bins = histbins_sigma, density = False)
+ cumhist = 1. - np.cumsum(hist) / float(len(indx))
+
+ return np.log10(histbins_sigma[:-1]), cumhist
+
+
+# #################################################################
+# Adds stellar surface density image to subplot
+# #################################################################
+def add_stellar_surface_density(ax, filename):
+ data = load(filename)
+ t_Myr = data.metadata.time.to('Myr')
+ data.metadata.boxsize.convert_to_units("kpc")
+ boxsize = data.metadata.boxsize
+
+ data.stars.birth_times.convert_to_units('Myr')
+ data.stars.coordinates.convert_to_units("kpc")
+ data.stars.masses.convert_to_units("Msun")
+ Mstar = np.sum(data.stars.masses)
+
+ data.metadata.time.convert_to_units('Myr')
+ tmin = np.maximum(data.metadata.time.value - MaxAgeYoungStarsMyr, 0.)
+
+ data.stars.masses = data.stars.masses[data.stars.birth_times >= tmin]
+ data.stars.coordinates = data.stars.coordinates[data.stars.birth_times >= tmin, :]
+
+ total_usermass = np.sum(data.stars.masses)
+
+ data.stars.smoothing_lengths = generate_smoothing_lengths(
+ coordinates=data.stars.coordinates,
+ boxsize=data.metadata.boxsize,
+ kernel_gamma=kernel_gamma,
+ neighbours=11,
+ speedup_fac=1,
+ dimension=3,
+ )
+
+ deltax = unyt_quantity(r_img_kpc, 'kpc')
+ deltay = unyt_quantity(r_img_kpc, 'kpc')
+
+ range_plot = [0.5 * boxsize[0] - deltax, 0.5 * boxsize[0] + deltax , 0.5 * boxsize[1] - deltay, 0.5 * boxsize[1] + deltay]
+ extent = [range_plot[0].value-0.5 * boxsize[0].value, \
+ range_plot[1].value-0.5 * boxsize[0].value, \
+ range_plot[2].value-0.5 * boxsize[1].value, \
+ range_plot[3].value-0.5 * boxsize[1].value ]
+
+ pixel_grid = project_pixel_grid(
+ data.stars,
+ boxsize=data.metadata.boxsize,
+ resolution=int(npix),
+ project="masses",
+ region=range_plot,
+ )
+ x_range = range_plot[1] - range_plot[0]
+ y_range = range_plot[3] - range_plot[2]
+ units = 1.0 / (x_range * y_range)
+ # Unfortunately this is required to prevent us from {over,under}flowing
+ # the units...
+ units.convert_to_units(1.0 / (x_range.units * y_range.units))
+ units *= getattr(data.stars, 'masses').units
+
+ mass_map = unyt_array(pixel_grid, units=units)
+ mass_map.convert_to_units(msun / kpc**2)
+
+ mass_map_plot = np.log10(mass_map.value)
+
+ nonzero = mass_map_plot >= 1e-2
+ min_nonzero = np.min(mass_map_plot[nonzero])
+ mass_map_plot[~nonzero] = min_nonzero
+
+ cblab = 'Star'
+ cmap_loc = matplotlib.cm.get_cmap(cmap_stars)
+ tcolor = 'white'
+ ccolor = 'black'
+
+
+ # mask with circle
+ lx, ly = mass_map_plot.shape
+ X, Y = np.ogrid[0:lx, 0:ly]
+ mask = (X - lx / 2) ** 2 + (Y - ly / 2) ** 2 > lx * ly / 4
+ mass_map_plot[mask] = np.nan
+
+ # make plot
+ ax.tick_params(labelleft = False, labelbottom = False)
+ ax.set_xticks([])
+ ax.set_yticks([])
+ ax.set_axis_off()
+
+ im = ax.imshow(mass_map_plot, cmap=cmap_stars, extent = extent, vmin = vmin, vmax = vmax, origin = 'lower')
+ circle = plt.Circle((0, 0), 0.99 * r_img_kpc, color=ccolor, fill=False, linewidth = 2)
+ ax.add_artist(circle)
+
+ return mass_map_plot
+
+def add_fof_groups(ax, foffile):
+
+ with h5py.File(foffile, 'r') as f:
+ group_centers_kpc = f["Groups/Centres"][:]
+ group_masses_msun = f["Groups/Masses"][:] * f["Groups/Masses"].attrs["Conversion factor to physical CGS (including cosmological corrections)"] / msun
+ boxsize_kpc = f["Header"].attrs["BoxSize"]
+
+ group_masses_msun.convert_to_units('1/g')
+ group_masses_msun = group_masses_msun.value
+
+ group_centers_kpc[:,0] = group_centers_kpc[:,0] - 0.5 * boxsize_kpc[0]
+ group_centers_kpc[:,1] = group_centers_kpc[:,1] - 0.5 * boxsize_kpc[1]
+ group_centers_kpc[:,2] = group_centers_kpc[:,2] - 0.5 * boxsize_kpc[2]
+
+ r = np.sqrt( group_centers_kpc[:,1]**2 + group_centers_kpc[:,0]**2 )
+
+ ax.scatter(group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun < 1.e8), 1], \
+ group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun < 1.e8), 0], marker = 'o', edgecolor = 'yellow', facecolor = 'none',
+ s = 110, linewidths = 2.5)
+
+ ax.scatter(group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun > 1.e8), 1], \
+ group_centers_kpc[(r>mask_central_kpc) & (group_masses_msun > 1.e8), 0], marker = 'D', edgecolor = 'yellow', facecolor = 'none', \
+ s = 220, linewidths = 2.5)
+ return
+
+