From 684b9aa1c756a5e2ca38b79274ba91bdcd9143a6 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 8 Feb 2017 11:59:44 +0000
Subject: [PATCH] More simplifications of script to generate ICs
---
examples/ExternalPointMass/makeIC.py | 19 ++++++++++---------
1 file changed, 10 insertions(+), 9 deletions(-)
diff --git a/examples/ExternalPointMass/makeIC.py b/examples/ExternalPointMass/makeIC.py
index 8bda95681a..ba415daf9e 100644
--- a/examples/ExternalPointMass/makeIC.py
+++ b/examples/ExternalPointMass/makeIC.py
@@ -49,10 +49,10 @@ print 'G in internal units: ', const_G
# Parameters
-periodic= 1 # 1 For periodic box
-boxSize = 100. #
-Radius = boxSize / 4. # maximum radius of particles
-Mass = 1e10
+periodic = 1 # 1 For periodic box
+boxSize = 100. #
+max_radius = boxSize / 4. # maximum radius of particles
+Mass = 1e10
print "Mass at the centre: ", Mass
numPart = int(sys.argv[1]) # Number of particles
@@ -95,22 +95,23 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
grp1 = file.create_group("/PartType1")
#generate particle positions
-radius = Radius * (numpy.random.rand(numPart))**(1./3.)
-ctheta = -1. + 2 * numpy.random.rand(numPart)
-stheta = numpy.sqrt(1.-ctheta**2)
-phi = 2 * math.pi * numpy.random.rand(numPart)
+radius = max_radius * (numpy.random.rand(numPart))**(1./3.)
+print '---------------------'
+print 'Radius: minimum = ',min(radius)
+print 'Radius: maximum = ',max(radius)
+radius = numpy.sort(radius)
r = numpy.zeros((numPart, 3))
r[:,0] = radius
#generate particle velocities
speed = numpy.sqrt(const_G * Mass / radius)
-v = numpy.zeros((numPart, 3))
omega = speed / radius
period = 2.*math.pi/omega
print '---------------------'
print 'Period: minimum = ',min(period)
print 'Period: maximum = ',max(period)
+v = numpy.zeros((numPart, 3))
v[:,0] = -omega * r[:,1]
v[:,1] = omega * r[:,0]
--
GitLab