diff --git a/Fig2_SPH_density_sketch.py b/Fig2_SPH_density_sketch.py
new file mode 100644
index 0000000000000000000000000000000000000000..a14176f22c506213d262c44f2107277b5f1adc3e
--- /dev/null
+++ b/Fig2_SPH_density_sketch.py
@@ -0,0 +1,184 @@
+from unyt import cm, Msun, pc
+from unyt import proton_mass_cgs as mH
+import numpy as np
+import random
+import matplotlib.pyplot as plt
+import matplotlib.patches as patches
+
+############################################################################
+# Produce plots for various gas densities, represented by nparticles
+############################################################################
+nH_highdens_arr = np.asarray([1.e2, 1.e3, 1.e4, 1.e5, 1.e6])* cm**(-3)
+nparticles_arr = np.asarray([100, 100000])
+kerneltype = "WendlandSWIFT" # options: Wendland, WendlandSWIFT, QuarticSWIFT
+
+############################################################################
+# Simulation settings
+############################################################################
+mB = 1.e5 * Msun # mass per gas particle
+XH = 0.74 # hydrogen mass fraction
+hmin_in_pc = 1.55 * 10. # minimum smoothing length in pc
+
+############################################################################
+# Creates random positions within radius for nparticles
+############################################################################
+def create_particle_position(nparticles, radius):
+ i = 0
+ x = np.zeros(nparticles)
+ y = np.zeros(nparticles)
+ z = np.zeros(nparticles)
+ r = np.zeros(nparticles)
+ while i < nparticles:
+ xx = -radius + 2. * radius * random.random()
+ yy = -radius + 2. * radius * random.random()
+ zz = -radius + 2. * radius * random.random()
+
+ rr = np.sqrt(xx*xx + yy*yy + zz*zz)
+ if rr <= radius:
+ x[i] = xx
+ y[i] = yy
+ z[i] = zz
+ r[i] = rr
+ i = i+1
+ return x, y, z, r
+
+############################################################################
+# Kernel function to calculate the density: Wendland C2
+# source: https://pysph.readthedocs.io/en/latest/reference/kernels.html
+# class pysph.base.kernels.WendlandQuintic
+############################################################################
+def kernel(q,h):
+ alpha = 21./ (16. * np.pi * h**3)
+ k = alpha * np.power(1. - q/2, 4) * (2.*q + 1.)
+ k[q>2] = 0.
+ return k
+
+############################################################################
+# Kernel function to calculate the density: Wendland C2
+# source: SWIFT paper
+############################################################################
+def wendland_kernel(u,h):
+ kernel_gamma = 1.936492
+ q = u / kernel_gamma
+ alpha = 21. / (2. * np.pi * h**3) / (kernel_gamma * kernel_gamma * kernel_gamma)
+ k = alpha * (1. + q*q * (-10. + q * (20. + q * (-15. + q*4.))))
+ k[q>1] = 0.
+ return k
+
+############################################################################
+# Kernel function to calculate the density: quartic spline
+############################################################################
+def quartic_spline_kernel(u,h):
+ kernel_gamma = 2.018932
+ q = u / kernel_gamma
+ alpha = 15625. / (np.pi * 512. * h**3) / (kernel_gamma * kernel_gamma * kernel_gamma)
+ k = alpha * ( 0.368 + q*q * (-2.4 + q*q* 6.))
+ k[(q>=0.2)&(q<0.6)] = alpha * ( 0.352 + q[(q>=0.2)&(q<0.6)] * (0.32 + \
+ q[(q>=0.2)&(q<0.6)] * (-4.8 + \
+ q[(q>=0.2)&(q<0.6)] * (8.00 - \
+ q[(q>=0.2)&(q<0.6)] * 4.))))
+ k[(q>=0.6)&(q<1.0)] = alpha * ( 1.000 + q[(q>=0.6)&(q<1.0)] * (-4.0 + \
+ q[(q>=0.6)&(q<1.0)] * ( 6.0 + \
+ q[(q>=0.6)&(q<1.0)] * (-4. + \
+ q[(q>=0.6)&(q<1.0)] ))))
+ k[q>=1.] = 0.
+ return k
+
+############################################################################
+# Calculate SPH density at the center (x = y = z = 0) of the particle distribution
+############################################################################
+def get_SPH_density(r, mB, h):
+ q = r/h
+ if kerneltype == "Wendland":
+ W = kernel(q,h)
+ elif kerneltype == "WendlandSWIFT":
+ W = wendland_kernel(q,h)
+ elif kerneltype == "QuarticSWIFT":
+ W = quartic_spline_kernel(q,h)
+ else:
+ print ("UNKNOWN kerneltype; default to WendlandSWIFT")
+ W = wendland_kernel(q,h)
+ return mB * np.sum(W)
+
+############################################################################
+# Make individual plot for one density and one particle number
+############################################################################
+def make_plot(x, y, hsmooth, dens_input, dens_output, hcorrect, hmin, nparticles):
+ xsize = 150.
+ csize = xsize/20.
+ fig = plt.figure(figsize = (2.5,2.5))
+ fig.subplots_adjust(bottom = 0.02, top =0.98, right = 0.98,left = 0.02)
+ ax = fig.add_subplot(111)
+ ax.set_aspect('equal')
+ ax.set_xlabel('x [pc]')
+ ax.set_ylabel('y [pc]')
+ ax.set_xlim(-xsize/2., xsize/2.)
+ ax.set_ylim(-xsize/2., xsize/2.)
+ ax.xaxis.set_ticks(np.arange(-xsize/2., xsize/2. + 50., 50.))
+ ax.yaxis.set_ticks(np.arange(-xsize/2., xsize/2. + 50., 50.))
+
+ ax.axes.xaxis.set_visible(False)
+ ax.axes.yaxis.set_visible(False)
+
+ ax.scatter(x,y, color = 'black', s = 2, alpha = 0.3)
+
+ rect1 = patches.Rectangle((-75., 45), 60, 30, linewidth=1, edgecolor='grey', facecolor='white', alpha = 1.0)
+ rect2 = patches.Rectangle((-75., -75), 60, 30, linewidth=1, edgecolor='grey', facecolor='white', alpha = 1.0)
+ rect3 = patches.Rectangle((15., 45), 60, 30, linewidth=1, edgecolor='grey', facecolor='white', alpha = 1.0)
+ rect4 = patches.Rectangle((15., -75), 60, 30, linewidth=1, edgecolor='grey', facecolor='white', alpha = 1.0)
+
+ ax.add_patch(rect1)
+ ax.add_patch(rect2)
+ ax.add_patch(rect3)
+ ax.add_patch(rect4)
+
+ circle2 = plt.Circle((0.0, 0.0), hcorrect, color='#CC0066', fill = False, lw = 3)
+ circle3 = plt.Circle((0.0, 0.0), hmin, color='#0066CC', lw = 3, ls = 'dashed', fill = False)
+ circle4 = plt.Circle((0.0, 0.0), hmin, color='white', lw = 3, ls = 'solid', fill = False)
+ ax.add_artist(circle4)
+ ax.add_artist(circle2)
+ ax.add_artist(circle3)
+
+ ax.plot([30,60],[-65,-65], color = 'black', linestyle = 'solid')
+ ax.text(45., -63., '30 pc', va = 'bottom', ha = 'center', fontsize = 14)
+
+ circle2leg = plt.Circle((-xsize/2. + csize, -xsize/2. + 2.5 * csize), csize /2., color = '#CC0066', fill = False)
+ circle3leg = plt.Circle((-xsize/2. + csize, -xsize/2. + 1.0 * csize), csize /2., color = '#0066CC', ls = 'dashed', fill = False)
+ ax.add_artist(circle2leg)
+ ax.add_artist(circle3leg)
+
+ ax.text(-xsize/2. + 2. * csize, -xsize/2. + 2.0 * csize, 'h$_{\mathrm{correct}}$', fontsize = 14)
+ ax.text(-xsize/2. + 2. * csize, -xsize/2. + 0.5 * csize, 'h$_{\mathrm{min}}$', fontsize = 14)
+
+ ax.text(0.01, 0.99, 'log n$_{\mathrm{H}}$' + '\n' + '%.2f'%(np.log10(dens_input)), ha = 'left', va = 'top', transform=ax.transAxes, fontsize = 14)
+ ax.text(0.99, 0.99, 'log n$_{\mathrm{H,SPH}}$' + '\n' + '%.2f'%(np.log10(dens_output)), ha = 'right', va = 'top', transform=ax.transAxes, fontsize = 14)
+
+ outputname = "Fig2_particles%i_hsmooth_logn_input%.2f_log_hsmooth%.2f"%(nparticles, np.log10(dens_input), np.log10(hsmooth))
+ outputname = outputname.replace(".", "p") + ".png"
+ fig.savefig(outputname, dpi = 150)
+ plt.close()
+ print ("Figure saved: "+outputname)
+
+ return
+
+############################################################################
+# Make plots for all densities and all particle numbers
+############################################################################
+def make_all_plots():
+ for nH_highdens in nH_highdens_arr:
+ for nparticles in nparticles_arr:
+ radius_highdens = np.power( 3. * nparticles * mB * XH / (4. * np.pi * mH * nH_highdens) , 1./3.)
+ radius_highdens.convert_to_units('pc')
+
+ # this is the correct smoothing length
+ h = np.power(1./ nH_highdens * XH / mH * mB, 1./3.) * 1.2348
+ h.convert_to_units('pc')
+
+ h_min = np.maximum(hmin_in_pc, h.value)
+ x, y, z, r = create_particle_position(nparticles, radius_highdens)
+ dens = get_SPH_density(r * pc, mB, h_min * pc) / mH
+ dens.convert_to_units('cm**(-3)')
+ make_plot(x, y, h_min, nH_highdens, dens.value, h.value, hmin_in_pc, nparticles)
+
+
+make_all_plots()