diff --git a/Fig9a_gallery_gas.py b/Fig9a_gallery_gas.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b09c77599bc2e18f498e753be9e24cd055e7e1a
--- /dev/null
+++ b/Fig9a_gallery_gas.py
@@ -0,0 +1,146 @@
+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 plotannotation import add_arrows
+
+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 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 simulation_path import snapshotnumber
+
+
+snapshotnumber_rerun = 0
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig9a_gallery_gas.png"
+
+############################################################################
+# Some extra plotting settings
+############################################################################
+
+axis_xsize = 0.355
+axis_ysize = 0.28
+axis_dy = 0.010
+axis_dx_left = 0.07
+axis_dx = 0.2
+axis_dy_bottom = 0.11
+
+xticks = np.arange(-4,8,2)
+yticks = np.arange(1,5,1)
+
+# #################################################################
+# Creates 5 phasespace plots of re-calculated density vs. temperature
+# of the 5 runs identified in helpers/Fig7_runs_to_plot.py
+# Labels need to be updated if different runs are selected
+# #################################################################
+
+def make_phasespace_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 = (10*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)
+
+
+ #######################################################
+ # Setup axes
+ #######################################################
+ ax = []
+ # 0: top left panel
+ ax.append(fig.add_axes([axis_dx_left, axis_dy_bottom + 2.* axis_dy + 2.*axis_ysize, axis_xsize, axis_ysize]))
+ # 1: top right panel
+ ax.append(fig.add_axes([axis_dx_left + axis_dx + axis_xsize, axis_dy_bottom + 2.* axis_dy + 2.*axis_ysize, axis_xsize, axis_ysize]))
+ # 2: middle panel
+ ax.append(fig.add_axes([0.5 * axis_dx_left + 0.5 - 0.5 * axis_xsize, axis_dy_bottom + 1.* axis_dy + 1.*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([axis_dx_left + axis_xsize + axis_dx, 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 = 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].set_xticks(xticks)
+ ax[iplot].set_yticks(yticks)
+
+ if iplot == 0:
+ ax[iplot].set_ylabel('log T [K]')
+ ax[iplot].yaxis.set_tick_params(labelsize=16)
+ ax[iplot].set_xticklabels([])
+ elif iplot == 3:
+ ax[iplot].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ ax[iplot].set_ylabel('log T [K]')
+ ax[iplot].xaxis.set_tick_params(labelsize=16)
+ ax[iplot].yaxis.set_tick_params(labelsize=16)
+ elif iplot == 4:
+ ax[iplot].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ ax[iplot].xaxis.set_tick_params(labelsize=16)
+ ax[iplot].set_yticklabels([])
+ else:
+ ax[iplot].set_xticklabels([])
+ ax[iplot].set_yticklabels([])
+
+
+ for iax in [0,1,2,4]:
+ ax[iax].text(0.98, 0.5, "log M$_{\mathrm{zone}}$", ha = 'right', va = 'baseline', \
+ transform = ax[iax].transAxes, fontsize = 'x-large')
+
+ 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.5, "Increase $\epsilon$\n (x 2)", transform = ax[2].transAxes, fontsize = 'x-large', va = 'bottom', ha = 'center')
+ ax[2].text(0.5,-0.5, "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_phasespace_gallery_plot(outputname)
diff --git a/helpers/plotannotation.py b/helpers/plotannotation.py
new file mode 100644
index 0000000000000000000000000000000000000000..29e2c741560b41d54d246548f754378bc6970904
--- /dev/null
+++ b/helpers/plotannotation.py
@@ -0,0 +1,63 @@
+from matplotlib import patches
+
+
+############################################################################
+# Add arrows to connect subplots
+############################################################################
+def add_arrows(fig,ax):
+ # from center to top right
+ arrow0 = patches.ConnectionPatch(
+ [0.5, 1.0],
+ [0., 0.5],
+ coordsA=ax[2].transAxes,
+ coordsB=ax[1].transAxes,
+ color="black",
+ arrowstyle="-|>", # "normal" arrow
+ mutation_scale=30, # controls arrow head size
+ linewidth=3,
+ )
+
+ # from center to bottom right
+ arrow1 = patches.ConnectionPatch(
+ [1.0, 0.5],
+ [0.5, 1.0],
+ coordsA=ax[2].transAxes,
+ coordsB=ax[4].transAxes,
+ color="black",
+ arrowstyle="-|>", # "normal" arrow
+ mutation_scale=30, # controls arrow head size
+ linewidth=3,
+ )
+
+ # from center to bottom left
+ arrow2 = patches.ConnectionPatch(
+ [0.5, 0.0],
+ [1.0, 0.5],
+ coordsA=ax[2].transAxes,
+ coordsB=ax[3].transAxes,
+ color="black",
+ arrowstyle="-|>", # "normal" arrow
+ mutation_scale=30, # controls arrow head size
+ linewidth=3,
+ )
+
+ # from center to top left
+ arrow3 = patches.ConnectionPatch(
+ [0.0, 0.5],
+ [0.5, 0.0],
+ coordsA=ax[2].transAxes,
+ coordsB=ax[0].transAxes,
+ color="black",
+ arrowstyle="-|>", # "normal" arrow
+ mutation_scale=30, # controls arrow head size
+ linewidth=3,
+ )
+
+ fig.patches.append(arrow0)
+ fig.patches.append(arrow1)
+ fig.patches.append(arrow2)
+ fig.patches.append(arrow3)
+
+ return
+
+
diff --git a/helpers/select_simulations_for_figures.py b/helpers/select_simulations_for_figures.py
index 3ea8c5415af41c476a711ffc9526d5270e4afe95..913b365afb56ed41a0482ae55effc7002c700b8c 100644
--- a/helpers/select_simulations_for_figures.py
+++ b/helpers/select_simulations_for_figures.py
@@ -16,3 +16,18 @@ 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)
+
+#################################################
+# Simulations to plot for Figure 9
+# 5 different simulations, if different
+# simulations are selected, some annotations might
+# have to be manually updated
+#################################################
+
+Fig9_mB = ["M5", "M5", "M5", "M5", "M6"]
+Fig9_eps = [0.25, 0.5, 0.25, 0.25, 0.25]
+Fig9_hmin = [0.0775, 0.0775, 0.0775, 0.00775, 0.0775]
+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)
+