diff --git a/Fig1_Jeans_mass_comparison.py b/Fig1_Jeans_mass_comparison.py
new file mode 100644
index 0000000000000000000000000000000000000000..8fd8705db4bebd1313c4195f77b608d19e2edd0a
--- /dev/null
+++ b/Fig1_Jeans_mass_comparison.py
@@ -0,0 +1,122 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.lines import Line2D
+
+import matplotlib.pylab as pylab
+params = {'legend.fontsize': 'large',
+ 'axes.labelsize': 'x-large',
+ 'axes.titlesize':'x-large',
+ 'xtick.labelsize':'large',
+ 'ytick.labelsize':'large'}
+pylab.rcParams.update(params)
+
+
+import sys
+sys.path.insert(0, "helpers/")
+from equations import Jeans_length_0
+from equations import Jeans_mass_0
+from equations import Jeans_mass_P
+from equations import Jeans_mass_W
+
+
+############################################################################
+# Set Plummer softening scale in pc
+############################################################################
+eps = 100. # in pc
+H = 3. * eps
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = 'Fig1_Jeans_mass_softening_%4.4ipc.png'%(round(eps))
+
+############################################################################
+# set plot limits and define 2D arrays for contour plots
+############################################################################
+lognH_min = -6.
+lognH_max = 6.
+dlognH = 0.1
+
+logT_min = 1.
+logT_max = 6.
+dlogT = 0.1
+
+LJmin = np.log10(1.)
+LJmax = np.log10(10. * 1000.)
+
+MJmin = np.log10(1.e1)
+MJmax = np.log10(1.e10)
+
+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]]
+
+# manual locations for contour labels
+manual_locations = [ (4.5,1.5), \
+ (3.0,1.5), \
+ (1.0,1.5), \
+ (-1.,1.5), \
+ (-2.0,1.5), \
+ (-4.5,1.5), \
+ (-2.0,3.), \
+ (-4.5,3.), \
+ (-2.0,4.5), \
+ (-4.5,4.) ]
+
+ls = ['solid', 'dotted']
+lines = [Line2D([0], [0], color='black', linewidth=2, linestyle=l) for l in ls]
+labels = ['Plummer', 'Wendland C2']
+
+linesN = [Line2D([0], [0], color='black', linewidth=2, linestyle='solid')]
+labelsN = ['Newton']
+
+linesL = [Line2D([0], [0], color='#CC0066', linewidth=2, linestyle='dashed')]
+labelsL = ['$\lambda_{\mathrm{J,N}} \equiv \epsilon = %i$ pc'%(round(eps))]
+
+############################################################################
+# get Jeans masses and Newtonian Jeans lengths for full phase-space
+############################################################################
+logMJ_0_2D = np.log10(Jeans_mass_0(nH_2Darr, T_2Darr))
+logMJ_P_2D = np.log10(Jeans_mass_P(nH_2Darr, T_2Darr, eps))
+logMJ_W_2D = np.log10(Jeans_mass_W(nH_2Darr, T_2Darr, H))
+logLJ_0_2D = np.log10(Jeans_length_0(nH_2Darr, T_2Darr))
+
+############################################################################
+# make plot
+############################################################################
+
+def make_plot(outputname):
+ fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, squeeze=True, figsize = (5*1.1,4*1.1))
+ fig.subplots_adjust(top = 0.9, bottom = 0.15)
+ fig.suptitle('Jeans mass (log M$_{\mathrm{J}}$ [M$_{\odot}$])', fontsize = 14)
+ ax[0].set_ylabel('log T [K]')
+ CF_0 = ax[0].contour(lognH_2Darr, logT_2Darr, logMJ_0_2D, levels = np.arange(MJmin, MJmax + 1, 1), colors = 'black')
+ ax[0].clabel(CF_0, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations)
+ ax[0].legend(linesN,labelsN, loc = 'upper left', framealpha = 1.0)
+
+ ax[1].set_ylabel('log T [K]')
+ ax[1].set_xlabel('log n$_{\mathrm{H}}$ [cm$^{-3}$]')
+ CF_P = ax[1].contour(lognH_2Darr, logT_2Darr, logMJ_P_2D, levels = np.arange(MJmin, MJmax + 1, 1), colors = 'black')
+ CF_W = ax[1].contour(lognH_2Darr, logT_2Darr, logMJ_W_2D, levels = np.arange(MJmin, MJmax + 1, 1), linestyles = 'dotted', colors = 'black')
+ ax[1].clabel(CF_P, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations[5:], colors = 'white')
+ ax[1].clabel(CF_W, inline=1, fontsize=10, fmt='%1.1f',manual=manual_locations[5:])
+ legend1 = ax[1].legend(lines,labels, loc = 'upper left', framealpha = 1.0)
+
+ CF_extra = ax[1].contour(lognH_2Darr, logT_2Darr, logLJ_0_2D-np.log10(eps), levels = np.asarray([0.]), colors ='#CC0066', linestyles = 'dashed')
+ legend2 = ax[1].legend(linesL,labelsL, loc = 'lower right', framealpha = 1.0)
+ ax[1].add_artist(legend1)
+ ax[1].add_artist(legend2)
+
+ fig.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+make_plot(outputname)
diff --git a/helpers/equations.py b/helpers/equations.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ab6fea85a07a1e67d1c95e16e980bb7a6347db2
--- /dev/null
+++ b/helpers/equations.py
@@ -0,0 +1,78 @@
+import numpy as np
+from unyt import msun, pc
+from unyt import proton_mass_cgs as mH
+
+XH = 0.74
+eta_res = 1.2348
+
+############################################################################
+# nH ... gas density in cm-3
+# T ... gas temperature in K
+# returns the Newtonian Jeans length in pc
+############################################################################
+def Jeans_length_0(nH,T):
+ return 5.4 * np.power(T, 1./2.) * np.power(nH, -1./2.)
+
+############################################################################
+# nH ... gas density in cm-3
+# T ... gas temperature in K
+# returns the Newtonian Jeans mass in Msol
+############################################################################
+def Jeans_mass_0(nH, T):
+ return 21.2 * np.power(T, 3./2.) * np.power(nH, -1./2.)
+
+############################################################################
+# nH ... gas density in cm-3
+# T ... gas temperature in K
+# returns the Plummer softened Jeans mass in Msol
+############################################################################
+def Jeans_mass_P(nH, T, eps):
+ return Jeans_mass_0(nH, T) * np.power(1. + 1.42 * np.power(eps/Jeans_length_0(nH,T), 3./2.), 6./5.)
+
+############################################################################
+# nH ... gas density in cm-3
+# T ... gas temperature in K
+# returns the softened Jeans length (for Wendland C2) in pc
+############################################################################
+def Jeans_length_W(nH,T,eps):
+ H = 3. * eps
+ return Jeans_length_0(nH,T) * np.power(1. + 0.27 * np.power(H / Jeans_length_0(nH,T), 2.), 0.3)
+
+############################################################################
+# nH ... gas density in cm-3
+# T ... gas temperature in K
+# returns the Wendland C2 softened Jeans mass in Msol
+############################################################################
+def Jeans_mass_W(nH, T, H):
+ return Jeans_mass_0(nH, T) * np.power(1. + 0.27 * np.power(H/Jeans_length_0(nH,T), 2.), 9./10.)
+
+
+############################################################################
+# mB ... particle mass resolution in Msol
+# hmin ... minimum smoothing length in pc
+# returns the density above which the smoothing length is limited by hmin
+############################################################################
+def density_critical_hmin(mB, hmin):
+ mBsol = mB * msun
+ hminpc = hmin * pc
+ nH = XH / mH * np.power(eta_res,3) * mBsol / np.power(hminpc, 3.)
+ nH.convert_to_units('cm**-3')
+ return nH
+
+
+############################################################################
+# mB ... particle mass resolution in Msol
+# nH ... gas density in cm-3
+# hmin ... minimum smoothing length in pc
+# returns the smoothing length in pc
+############################################################################
+def h_smooth(mB, nH, hmin):
+ h_pc = 82. * np.power(mB / 1.e6, 1./3.) * np.power(nH / 100., -1./3.)
+ h_pc = np.maximum(h_pc, hmin)
+ return h_pc
+
+
+
+
+
+