# Sets up a uniform disc of particles in solid body rotation

import numpy
from numpy import *
import random
from random import random as rand
import math
from math import pi as pi

# Disc radius kpc
r = 8.65

# Disc thickness kpc
t = 1.0

# No. particles
Npart = 10000

# Sound speed km/s
cs = 6.5

# Total mass Msolar
mtot = 3.0E7

# Maximum speed km/s
v = 5.0

# Gamma
gamma = 1.666

# Initial position
xoffset = 0.0
yoffset = 0.0
zoffset = 0.0

# Initial overall velocity
vxoffset = 0.0
vyoffset = 0.0
yzoffset = 0.0


# Derived properties
mpart = mtot / Npart					# Particle mass
utherm = cs*cs*1e6/(gamma-1.0) 	# m^2/s^2
theta = (2.0*pi) / Npart			# Angular increment
linespace = (2.0*r) / Npart		# Linear separation
rhop = Npart / (pi*r*r)				# Particle density


outfile = open('stream.dat','w')

vecfile = open('vectors.txt','w')

#iseed = -1234

# Initial sweep positions
angle = 0.0
z     = r
ri = math.sqrt(1.0/(pi*rhop))

print ri

for i in range(0,Npart):

	# METHOD ONE - sweep out a rotating line, populating each line with a single particle. Doesn't work, gives
	# a cored distribution. Requires a higher change of angle to give particles near the centre an equal 
	# separation.
	#dist  = r*rand()
	#depth = t*rand() - (0.5*t)	# Random over the range -0.5t -> +0.5t
	
	#phi = (pi/2.0) - angle
	
	# Calculate velocity
	#vx = -(v*(dist/r))*cos(phi)
	#vy = 0.0
	#vz = (v*(dist/r))*sin(phi)
	
	# Find cartesian coordinates
	#x = dist*cos(angle)
	#z = dist*sin(angle)
	#y = depth
	
	
	# METHOD TWO - sweep out a line moving linearly across the disc, populating each one with a single particle.
	#ztemp = abs(z)
	#if ztemp > r:
	#	ztemp = r
	#xdist = math.sqrt((r*r) - (ztemp*ztemp))
	#x = 2.0*xdist*rand() - xdist
	
	#depth = t*rand() - (0.5*t)	# Random over the range -0.5t -> +0.5t
	#y = depth
		
	# METHOD THREE - sweep out a line which increments a different amount each time, so that the area sampled
	# is always equivalent to containing one particle with a uniform density
	# Try... except prevents divide by zero errors at edges
	#try:
	#	linespace = (1.0/(2.0*rhop * math.sqrt((r*r) - (ztemp*ztemp))))
	#except:
	#	linespace = (2.0*r) / Npart
		
		
		
	# METHOD FOUR - increase the radial value by an amount equal to sweeping out equal area rings each time,
	# populating each ring with a particle positioned at a random angle.
	angle = 2.0*pi*rand()
	phi = (pi/2.0) - angle
	x = ri*cos(phi)	
	z = ri*sin(phi)
	y = t*rand() - (0.5*t)
	
	rout = math.sqrt((1.0+rhop*pi*ri*ri)/(rhop*pi))
	deltar = rout - ri
	
	# Calculate velocity
	vx = -(v*(ri/r))*sin(phi)
	vy = 0.0
	vz = (v*(ri/r))*cos(phi)
	
	# Write the parameters to the file
	thing  = str(x+xoffset)+' '+str(y+yoffset)+' '+str(z+zoffset)+' '+str(vx)+' '+str(vy)+' '+str(vz)+' '+str(mpart)+' '+str(utherm)
	#outfile.write(str(i)+' '+str(dist)+' '+str(angle)+'\n')
	outfile.write(thing+'\n')	
	
	#angle = angle + theta
	z = z - linespace
	ri = ri + deltar