makeIC.py 8.32 KB
 Tom Theuns committed Aug 08, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 ############################################################################### # This file is part of SWIFT. # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk) # Tom Theuns (tom.theuns@durham.ac.uk) # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published # by the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with this program. If not, see . # ############################################################################## import h5py import sys import numpy import math import random import matplotlib.pyplot as plt # Generates a disk-patch in hydrostatic equilibrium # see Creasey, Theuns & Bower, 2013, for the equations: # disk parameters are: surface density sigma # scale height b # density: rho(z) = (sigma/2b) sech^2(z/b) # isothermal velocity dispersion = 0): # set seed for random number numpy.random.seed(1234) grp1 = file.create_group("/PartType1") radius = Radius * (numpy.random.rand(N))**(1./3.) ctheta = -1. + 2 * numpy.random.rand(N) stheta = numpy.sqrt(1.-ctheta**2) phi = 2 * math.pi * numpy.random.rand(N) r = numpy.zeros((numPart, 3))  Matthieu Schaller committed Aug 19, 2016 226   Tom Theuns committed Aug 08, 2016 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246  speed = vrot v = numpy.zeros((numPart, 3)) omega = speed / radius period = 2.*math.pi/omega print 'period = minimum = ',min(period), ' maximum = ',max(period) v[:,0] = -omega * r[:,1] v[:,1] = omega * r[:,0] ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd') ds[()] = r ds = grp1.create_dataset('Velocities', (numPart, 3), 'f') ds[()] = v v = numpy.zeros(1) m = numpy.full((numPart, ),10) ds = grp1.create_dataset('Masses', (numPart,), 'f') ds[()] = m m = numpy.zeros(1)  Matthieu Schaller committed Aug 19, 2016 247   Tom Theuns committed Aug 08, 2016 248 249 250 251 252 253 254 255  ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L') ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L') ds[()] = ids file.close() sys.exit()