# 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()