Skip to content
Snippets Groups Projects
Commit b03440e4 authored by Sylvia Plöckinger's avatar Sylvia Plöckinger
Browse files

Added Fig 2

parent b6d0dc0e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Fig2_SPH_density_sketch.py 0 → 100644
+ 184
− 0
View file @ b03440e4
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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment