# Super simple script to measure total area under an input curve and calculate
# how many particles of fixed mass should be in each bin to reach a specified
# total.
# This version corrects for the fact we're dealing with rings.
# "Idealised density profile" comes from the observational M33 data. The values
# have been normalised by the peak value in topcat. The inner values have been
# manually set to 1.0 and the radius truncated at 25 kpc.
# This script outputs how many particles we need in each ring. The next stage is
# to run "IsoDiscProfile" which creates a disc of particles with the specified
# density pofile. In principle this procedure should work with ANY radial density
# profile.

import math
from math import pi as pi

totalp = 20000

infile = open('IdealisedDensityProfile.txt','r')
outfile = open('ParticleProfile.txt','w')
outfile.write('#r/kpc nparticles'+'\n')

totalarea = 0.0	#

h = 1.0


# First find the total area
line = infile.readline()
#line = infile.readline()
#rf, rho = line.split()
#rf = float(rf)
#rho = float(rho)


for line in infile:
	rf, rho = line.split()
	rf = float(rf)
	rho = float(rho)
	ri = rf - 1.0	# Known fixed ring step size
	
	area = ((pi*rf*rf)-(pi*ri*ri))*rho	
	ri = rf
	
	totalarea = totalarea + area
	
	
# Now find the normalised area under each bin so that we can set the number
# of particles in each bin.	
infile.seek(0)
infile.readline()

#line = infile.readline()
#ri, rho = line.split()
#ri = float(ri)
#rho = float(rho)

nparticles = 0

ncount = 0

for line in infile:
	rf, rho = line.split()
	rf = float(rf)
	ri  = rf - 1.0
	rho = float(rho)
	
	area = ((pi*rf*rf)-(pi*ri*ri))*rho	
	
	#if ri < 25.0:
	#	nparticles =  int(totalp*(area/totalarea))*2
	#if ri >= 25.0:
	nparticles = int(totalp*(area/totalarea))
	
	ncount = ncount + nparticles
	print ncount

	outfile.write(str(rf)+' '+str(nparticles)+'\n')
	
	ri = rf
	
	
infile.close()
outfile.close()