diff --git a/Fig3_zone_of_avoidance.py b/Fig3_zone_of_avoidance.py
new file mode 100644
index 0000000000000000000000000000000000000000..f986f839a1fd20d065a59f10ef5606bc2fccd029
--- /dev/null
+++ b/Fig3_zone_of_avoidance.py
@@ -0,0 +1,238 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from unyt import msun, pc
+from unyt import proton_mass_cgs as mH
+import matplotlib.pylab as pl
+
+import sys
+sys.path.insert(0, "helpers/")
+
+from equations import Jeans_length_0
+from equations import Jeans_length_W
+from equations import density_critical_hmin
+
+
+eta_res = 1.2348
+XH = 0.74
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "Fig3_zone_of_avoidance.png"
+
+############################################################################
+# set resolution properties to plot
+# epsarr ... gravitational softening scale in pc
+# mBarr ... gas particle mass in Msol
+# h_min_ratio_arr ... minimum smoothing length,
+# defined as h_min = 1.55 * h_min_ratio * eps
+############################################################################
+epsarr = np.asarray([31.6, 100., 316., 700., 1000.])
+mBarr = np.asarray([1.81e5, 1.81e6, 1.81e6, 1.81e7])
+h_min_ratio_arr = np.asarray([0.01, 0.0316, 0.1, 0.316])
+
+lwfat = 4
+lwthin = 2
+
+############################################################################
+# equations from paper
+############################################################################
+# mB in units of solar mass
+# hmin in units of pc
+def zone_min_density(mB, hmin):
+ return 58. * (mB / 1.e6) * np.power(hmin / 100., -3)
+
+# hmin in units of pc
+# eps in units of pc
+# nH in units of cm-3
+def zone_max_temperature(hmin, eps, nH):
+ return 27. * np.power(hmin / 100., 5) * np.power(eps / 700., -3) * (nH / 100.)
+
+############################################################################
+# set plot limits and define 2D arrays for contour plots
+############################################################################
+idefault_eps = 3
+idefault_mB = 2
+idefault_hmin = 2
+
+ncolors = np.maximum(len(epsarr), len(mBarr))
+ncolors = np.maximum(ncolors, len(h_min_ratio_arr))
+c = pl.cm.copper(np.linspace(0,0.9,ncolors))
+ls= ['dotted', 'dashed', 'dotted', 'dashed', 'dotted', 'dashed', 'dotted', 'dashed']
+
+lognH_min = -6.
+lognH_max = 8.
+dlognH = 0.01
+
+logT_min = 1.
+logT_max = 6.
+dlogT = 0.1
+
+xticks = np.arange(-6., 10., 2)
+
+lognH_arr = np.arange(lognH_min, lognH_max + dlognH, dlognH)
+logT_arr = np.arange(logT_min , logT_max + dlogT, dlogT)
+logT_arr_large = np.arange(logT_min , logT_max + 4. + dlogT, dlogT)
+
+lognH_2Darr = np.tile(lognH_arr, (len(logT_arr), 1))
+logT_2Darr = (np.tile(logT_arr, (len(lognH_arr), 1))).T
+
+nH_2Darr = np.power(10., lognH_2Darr)
+T_2Darr = np.power(10., logT_2Darr)
+
+eplot = [lognH_arr[0], lognH_arr[-1], logT_arr[0], logT_arr[-1]]
+
+def logTeff_EAGLE(lognH):
+ return 1./3. * (lognH +1.) + np.log10(8000.)
+
+############################################################################
+# make plot
+############################################################################
+def make_plot(outputname):
+
+ fig, ax = plt.subplots(nrows=3, ncols=1, sharex=True, sharey=True, squeeze=True, figsize = (3*1.3,7.4))
+ fig.subplots_adjust(top = 0.915, bottom = 0.07 , left = 0.15, right = 0.95)
+ fig.suptitle('Runaway collapse zone')
+
+ for axloc in ax:
+ axloc.set_xlim(lognH_min, lognH_max)
+ axloc.set_xlim(logT_min, logT_max)
+
+ ################### Plot 0 ############################
+ ax[0].set_title('Vary baryon particle mass')
+ #######################################################
+ ax[0].set_ylabel('log T [K]')
+ ax[0].xaxis.set_ticks(xticks)
+ ax[0].set_ylim(logT_min, logT_max)
+ ieps = idefault_eps
+ ihmin = idefault_hmin
+ h = []
+ l = []
+ for imB in range(len(mBarr)):
+ eps = epsarr[ieps]
+ h_min = 1.55 * eps * h_min_ratio_arr[ihmin]
+ mB = mBarr[imB]
+
+
+ if imB == idefault_mB:
+ currentcolor = "black"
+ currentlinestyle = "solid"
+ lw = lwfat
+ else:
+ currentcolor = c[imB]
+ currentlinestyle = ls[imB]
+ lw = lwthin
+
+
+ ax[0].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+ np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '10$^{%.1f}$'%((np.log10(mBarr[imB]))))
+
+ ax[0].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle)
+
+
+ # add EAGLE EOS:
+ ax[0].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+ # with rotated annotation:
+ x = [ 0., 4. ]
+ y = [logTeff_EAGLE(0.), logTeff_EAGLE(4.)]
+ dy = y[1] - y[0]
+ dx = x[1] - x[0]
+ angle = np.rad2deg(np.arctan2(dy, dx))
+
+ ax[0].text(x[0], y[0]+0.1, "EAGLE T$_{\mathrm{eff}}$", ha = 'left', va = 'bottom',\
+ transform_rotates_text=True, rotation=angle, rotation_mode='anchor')
+
+ ax[0].text(0.05,0.05, '$\epsilon$ = %i pc'%(round(eps))+'\n'+\
+ 'h$_{\mathrm{min}}$ = %.1f pc'%(h_min),\
+ va = 'bottom', ha = 'left', transform=ax[0].transAxes)
+ ax[0].text(0.97, 0.05, "runaway\n collapse zone", va = 'bottom', ha = 'right', transform=ax[0].transAxes)
+ ax[0].legend(loc = 'upper left', title = 'm$_{\mathrm{B}}$ [M$_{\odot}$]')
+
+ ################### Plot 1 ############################
+ ax[1].set_title('Vary softening scale')
+ #######################################################
+ ax[1].set_ylabel('log T [K]')
+ ax[1].xaxis.set_ticks(xticks)
+ imB = idefault_mB
+ ihmin = idefault_hmin
+ h = []
+ l = []
+ for ieps in range(len(epsarr)):
+ eps = epsarr[ieps]
+ h_min = 1.55 * epsarr[idefault_eps] * h_min_ratio_arr[idefault_hmin]
+ mB = mBarr[imB]
+
+ if ieps == idefault_eps:
+ currentcolor = "black"
+ currentlinestyle = "solid"
+ lw = lwfat
+ else:
+ currentcolor = c[ieps]
+ currentlinestyle = ls[ieps]
+ lw = lwthin
+
+ ax[1].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+ np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '%i'%(round(epsarr[ieps])))
+
+ ax[1].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle)
+
+ # add EAGLE EOS:
+ ax[1].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+
+ ax[1].text(0.05,0.05, 'm$_{\mathrm{B}}$ = 10$^{%.1f}$ M$_{\odot}$'%(np.log10(mB))+'\n'+\
+ 'h$_{\mathrm{min}}$ = %.1f pc'%(h_min),\
+ va = 'bottom', ha = 'left', transform=ax[1].transAxes)
+ ax[1].legend(loc = 'upper left', title = '$\epsilon$ [pc]')
+
+
+ ################### Plot 2 ############################
+ ax[2].set_title('Vary minimum smoothing length')
+ #######################################################
+ ax[2].set_ylabel('log T [K]')
+ ax[2].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ ax[2].xaxis.set_ticks(xticks)
+ imB = idefault_mB
+ ieps = idefault_eps
+ h = []
+ l = []
+ for ihmin in range(len(h_min_ratio_arr)):
+ eps = epsarr[ieps]
+ h_min = 1.55 * eps * h_min_ratio_arr[ihmin]
+ mB = mBarr[imB]
+
+ if ihmin == idefault_hmin:
+ currentcolor = "black"
+ currentlinestyle = "solid"
+ lw = lwfat
+ else:
+ currentcolor = c[ihmin]
+ currentlinestyle = ls[ihmin]
+ lw = lwthin
+
+ ax[2].plot(lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))], \
+ np.log10(zone_max_temperature(h_min, eps, np.power(10., lognH_arr[lognH_arr>=np.log10(zone_min_density(mB, h_min))]))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle, label = '%.1f'%(h_min))
+
+ ax[2].vlines(np.log10(zone_min_density(mB, h_min)), ymin = 1., ymax = np.log10(zone_max_temperature(h_min, eps, zone_min_density(mB, h_min))),\
+ lw = lw, color = currentcolor, linestyle = currentlinestyle)
+
+ h, l = ax[2].get_legend_handles_labels()
+
+ # add EAGLE EOS:
+ ax[2].plot(lognH_arr[lognH_arr>=-1.], logTeff_EAGLE(lognH_arr[lognH_arr>=-1.]), color = 'grey', lw = 5, alpha = 0.4)
+
+ ax[2].text(0.05,0.05, 'm$_{\mathrm{B}}$ = 10$^{%.1f}$ M$_{\odot}$'%(np.log10(mB))+'\n' +\
+ '$\epsilon$ = %i pc'%(round(eps)),\
+ va = 'bottom', ha = 'left', transform=ax[2].transAxes)
+ ax[2].legend(h[::-1],l[::-1], loc = 'upper left', title = 'h$_{\mathrm{min}}$ [pc]')
+
+ fig.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+if __name__ == "__main__":
+ make_plot(outputname)