Added appendix figure

FigA1_Jeans_length_root_finding.py

+from scipy import optimize
+import numpy as np
+from scipy.optimize import curve_fit
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MultipleLocator
+
+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)
+
+colorP = 'black'
+colorC = 'black'
+colorM = '#CC0066'
+dP = 10 # show only every dP point
+dC = 10 # show only every dC point 
+
+############################################################################
+# Set output filename
+############################################################################
+outputname = "FigA1_Jeans_length_mass.png"
+
+############################################################################
+# x ... Plummer softened Jeans length over Newtonian Jeans length
+# f ... Plummer softening scale over Newtonian Jeans length 
+############################################################################
+def Plummer_func_to_find_root(x,f):
+    return np.power(x,10./3.) - np.power(x,2) -  np.power(2., 2./3.)*np.power(f,2)
+
+############################################################################
+# x ... Wendland C2 softened Jeans length over Newtonian Jeans length
+# f ... H (3 times Plummer equivalent softening scale) over Newtonian Jeans length
+############################################################################
+def Wendland_func_to_find_root(x,f):
+    return np.power(x,5) - np.power(x,3) - 1./7. * np.power(f,3)
+
+############################################################################
+# fit to the Plummer softened Jeans length over Newtonian Jeans length
+# f ... Plummer softening scale over Newtonian Jeans length 
+# a has been determined in a different fitting routine
+############################################################################
+def Plummer_fit(f):
+    a = 1.42
+    #best fit:
+    #a = 1.42313008
+    return np.power(1. + a * np.power(f, 3./2.), 2./5.)
+
+############################################################################
+# fit to the Wendland softened Jeans length over Newtonian Jeans length
+# f ... Wendland softening scale over Newtonian Jeans length 
+# a has been determined in a different fitting routine
+############################################################################
+def Wendland_fit(f):
+    a = 0.27
+    #best fit:
+    #a = 0.27381175
+    return np.power(1. + a * np.power(f, 2), 0.3)
+
+############################################################################
+# Find the value for softened over Newtonian Jeans length for an array of values
+# for the softening scale over the Newtonian Jeans length (farr) for 
+# Plummer softening (Parr) and Wendland C2 softening (Carr)
+############################################################################
+farr = np.logspace(-3,3,num = 200)
+Parr = np.zeros_like(farr)
+Carr = np.zeros_like(farr)
+
+for i in range(len(farr)):
+    Parr[i] = optimize.bisect(Plummer_func_to_find_root, 1.e-8 , 1.e8, args = ( farr[i] ))
+    Carr[i] = optimize.bisect(Wendland_func_to_find_root, 1.e-8 , 1.e8, args = ( farr[i] ))
+    
+############################################################################
+# make plot
+############################################################################
+
+def make_plot(outputname):
+    fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True, squeeze=True, figsize = (5*1.1,4*1.1), gridspec_kw=dict(height_ratios=[1.,0.5]))
+    plt.subplots_adjust(top = 0.85, bottom = 0.15, right = 0.9)
+    fig.suptitle('Jeans length and mass for softened potentials', fontsize = 14)
+
+    ax[0,0].set_title('Plummer potential')
+    ax[0,0].yaxis.set_major_locator(MultipleLocator(1.0))
+    ax[0,0].set_ylabel('log $\lambda_{\mathrm{J,soft}}$ / $\lambda_{\mathrm{J,0}}$')
+    ax[0,0].set_ylim(-0.5, 4)
+    ax[0,0].tick_params(labelleft=True, right = True)
+    ax[0,0].plot(np.log10(farr), np.log10(Plummer_fit(farr)), color = colorP, label = r'Fit: (1 + 1.42 $\left ( \frac{\epsilon}{\lambda_{\mathrm{J,0}}}\right )^{3/2}$)$^{2/5}$')
+    ax[0,0].plot(np.log10(farr), 3.*np.log10(Plummer_fit(farr)), color = colorM)
+    ax[0,0].scatter(np.log10(farr[::dP]), np.log10(Parr[::dP]), color = colorP)
+    ax[0,0].scatter(np.log10(farr[::dP]), 3. * np.log10(Parr[::dP]), color = colorM)
+
+    ax[1,0].set_xlabel('log $\epsilon$ / $\lambda_{\mathrm{J,0}}$')
+    ax[1,0].set_ylabel('exact / fit')
+    ax[1,0].set_ylim(0.9, 1.5)
+    ax[1,0].yaxis.set_minor_locator(MultipleLocator(0.1))
+    ax[1,0].plot(np.log10(farr), Parr / Plummer_fit(farr), color = colorP)
+    ax[1,0].plot(np.log10(farr), np.power(Parr / Plummer_fit(farr),3), color = colorM, label = 'Jeans mass')
+    ax[1,0].legend(loc = 'upper left')
+
+    ax[0,1].set_title('Wendland C2 potential')
+    ax[0,1].yaxis.set_label_position("right")
+    ax[0,1].yaxis.tick_right()
+    ax[0,1].set_ylabel('log $M_{\mathrm{J,soft}}$ / $M_{\mathrm{J,0}}$', color = colorM)
+    ax[0,1].yaxis.set_major_locator(MultipleLocator(1.0))
+    ax[0,1].set_ylim(-0.5, 4)
+    ax[0,1].tick_params(labelright=True, left = True)
+    ax[0,1].plot(np.log10(farr), np.log10(Wendland_fit(farr)), color = colorC, label = r'Fit: (1 + 0.27 $\left ( \frac{H}{\lambda_{\mathrm{J,0}}}\right )^{2}$)$^{0.3}$')
+    ax[0,1].plot(np.log10(farr), 3. * np.log10(Wendland_fit(farr)), color = colorM)
+    ax[0,1].scatter(np.log10(farr[::dC]), np.log10(Carr[::dC]), color = colorC)
+    ax[0,1].scatter(np.log10(farr[::dC]), 3. * np.log10(Carr[::dC]), color = colorM)
+
+    ax[1,1].set_xlabel('log H / $\lambda_{\mathrm{J,0}}$')
+    ax[1,1].plot(np.log10(farr), Carr / Wendland_fit(farr), color = colorC, label = 'Jeans length')
+    ax[1,1].plot(np.log10(farr), np.power(Carr / Wendland_fit(farr),3), color = colorM)
+    ax[1,1].set_ylim(0.9, 1.5)
+    ax[1,1].yaxis.set_minor_locator(MultipleLocator(0.1))
+    ax[1,1].legend(loc = 'upper left')
+
+    fig.savefig(outputname, dpi = 150)
+    plt.close()
+    print ("Figure saved: "+outputname)
+
+make_plot(outputname)