diff --git a/README.md b/README.md
index 35e4fe5321c8368377ea7cc552ff1d05ebba0cd4..58d88f1000f09c457fdcacc5a11a917fef6d2c1d 100644
--- a/README.md
+++ b/README.md
@@ -4,9 +4,11 @@ Git for summerstudent project on evolution of open cluster (see project descript
# TODO:
- [ ] Get nbody6 working and get to know the code
+ - [ ] Make appointment with max and go through things
- [ ] Make it fairly easy for the student to work with nbody6
- [ ] make output scripts
- [ ] Create environment to run things on. i.e. go to IT
- [ ] Ask someone to send us a couple of review papers on clusters
- [x] Write up a more solid form of the project description
-- [ ] stay in contact with Mark and Fabio, keep them in the loop of any further advice. Also discuss with them semi regularly about the topic.
+- [ ] Stay in contact with Mark and Fabio, keep them in the loop of any further advice. Also discuss with them semi regularly about the topic.
+
diff --git a/python_scripts_for_nbody6/readout_script.py b/python_scripts_for_nbody6/readout_script.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ba12d534861b30958ac36d07d9b93d46adac6f2
--- /dev/null
+++ b/python_scripts_for_nbody6/readout_script.py
@@ -0,0 +1,385 @@
+#!/usr/bin/python
+
+#17/06/2014 FC: change line308,col73 (rh=compute_rh) // before: np.array(alist[6:9]) - after: np.array(alist[6:9])*alist[2]
+
+# Initially written by Mark Gieles.
+# Modified by David Hendriks
+
+"""
+I've attached mark's script. It looks through the "out3" file (the primary output file from an Nbody6 simulation run) and writes a load of interesting properties to the HDF5 file format which can be read in Python.
+In order to read this file format, you'll need to import module h5py. You open the file with:
+f = h5py.File("path_to_item", 'r')
+The properties can then be accessed like:
+f['00001/stars/id'].value
+... which would retrieve the star particle IDs for the first snapshot. Some of the units are a bit weird and you'll need to subtract the density centre and core velocity (these are properties you can find in the HDF5 file) if you want to centre on the cluster itself.
+Best, -Matt
+"""
+
+
+################################
+
+# Code to convert NBODY6 OUT3 file to HDF5 file format
+import numpy as np
+import struct
+import h5py
+import sys
+import argparse
+
+import pycuda.autoinit
+import pycuda.driver as drv
+from pycuda.compiler import SourceModule
+
+import matplotlib.pyplot as plt
+import os
+import argparse
+
+#Use GPU with cuda support
+dev=drv.Device(1)
+dev.make_context()
+
+def comp_pot_gpu(m, pos):
+ #Function to compute potential energy on the GPU using PyCUDA, given m,x,y,z
+ mod = SourceModule("""
+ __global__ void compute_potential_gpu(float *m, float *x, float *y, float *z, float *pot,int N) {
+ int i,j;
+ float potential, G = 1.0;
+
+ i = threadIdx.x + blockIdx.x*blockDim.x;
+ if (i<N) {
+ potential = 0.0f;
+ for (j=0; j<N; j++) {
+ if (j!=i){
+ potential -= G*m[i]*m[j]/sqrt(powf(x[i]-x[j],2.0) + powf(y[i]-y[j],2.0) + powf(z[i]-z[j],2.0));
+ }
+ }
+ pot[i] = potential;
+ }
+ }
+ """)
+
+ compute_potential_gpu = mod.get_function("compute_potential_gpu")
+ x = pos[:,0].astype(np.float32)
+ y = pos[:,1].astype(np.float32)
+ z = pos[:,2].astype(np.float32)
+ m = m.astype(np.float32)
+
+ pot = np.zeros_like(m)
+ N = len(m)
+
+ compute_potential_gpu(drv.In(m), drv.In(x), drv.In(y), drv.In(z), drv.Out(pot), np.int32(N), block=(1024,1,1), grid=(5000,1))
+ return pot
+
+def get_kin_pot(m,pos,vel):
+ # Compute kinetic and potential energy of all stars, ignoring tides
+ vcom = np.sum(np.multiply(m,vel.T).T,axis=0)/np.sum(m)
+ kin = 0.5*m*np.sum(np.square(vel-vcom), axis=1)
+ pot = comp_pot_gpu(m,pos)
+ return kin,pot,kin+pot,vcom
+
+def compute_energies(m,pos,vel,rdens,nbin):
+ # Compute bound stars, in 3 iterations
+
+ # Distances to density centres
+ r = np.sqrt(np.sum(np.square(pos - rdens), axis=1))
+
+ # 1st iteration: include only stars not too far from density centre and exclude binary members
+ c1 = (m>0)&(r<100) # This 500 values is a bit arbitrary and potentially tricky
+ c1[0:2*nbin]=False
+
+ # Compute energy and centre of mass velocity (vcom)
+ kin,pot,etot,vcom = get_kin_pot(m[c1],pos[c1,:],vel[c1,:])
+
+ # 2nd iteration: keep bound stars
+ c2 = (etot < 0)
+ kin,pot,etot,vcom = get_kin_pot(m[c1][c2],pos[c1,:][c2],vel[c1,:][c2])
+
+ # 3d iteration: keep bound stars
+ c3 = (etot<0)
+ bound=np.zeros(len(m),dtype='bool')
+
+ # Create array with 0 [unbound] and 1 [bound]
+ ind=np.flatnonzero(c1)
+ ind=ind[c2]
+ ind=ind[c3]
+ bound.flat[ind]=1
+
+ nbound = sum(c3)
+ mbound = sum(m[c1][c2][c3])
+
+ kinbound = sum(kin[c3])
+ potbound = sum(pot[c3])/2
+ return nbound,mbound,vcom,bound,kinbound,potbound
+
+def compute_rh(mass, pos, rdens, kstar, bound, sev):
+ # Compute half-mass radius of bound stars
+ r = np.sqrt(np.sum(np.square(pos - rdens), axis=1))
+
+ rb = r[bound]
+ mb = mass[bound]
+ args = np.argsort(rb)
+
+ ms = mb[args]
+ rs = rb[args]
+ cumm = np.cumsum(ms)
+
+
+ c300 = (r<300)
+ cn = (bound)
+ if (sev):
+ cn = (bound)&(kstar<=9)&(mass>=0.5)
+ c300 = (r<300)&(kstar<=9)&(mass>=0.5)
+
+ rhn = np.interp(0.5*sum(cn),np.arange(sum(cn)),np.sort(r[cn]))
+ rhn300 = np.interp(0.5*sum(c300),np.arange(sum(c300)),np.sort(r[c300]))
+ rhm = np.interp(0.5*cumm[-1],cumm,rs)
+
+# print " TEST ",rhn, rhm, cumm[-1], sum(cn), min(ms), max(ms)
+ return rhn,rhn300,rhm
+
+def read_nbody6(f5, file, sev):
+ """Read in the nbody6 output file(s).
+
+ Keyword arguments:
+ filename -- the OUT3 file generated by nbody6
+
+ """
+ # open NBODY6 file
+ f = open(file,'r')
+
+ # start looping over time and reading the nbody6 data.
+ # itime tracks the number of timesteps outputted.
+ itime = -1
+
+ # runs on some computers seem to generate extra bytes after the record length.
+ # 'extra' and its associated reads below take care of that.
+ extra = 0
+
+ byte = f.read(4) # begin header block
+ model = 0
+
+ while byte != "":
+ blocksize1 = struct.unpack('i',byte)[0] # header block size
+ if extra == 1:
+ f.read(4)
+ itime = itime + 1
+ # read the header for this time
+ ntot = struct.unpack('i',f.read(4))[0] #ntot
+ if ntot < 10 and itime == 0:
+ # check for a nonsensical ntot value,
+ # which is symptomatic of bizarre extra record padding that occurs on some
+ # systems. if needed, include extra reads to deal with these.
+ extra = 1
+ ntot = struct.unpack('i',f.read(4))[0]
+ print 'extra! '+str(extra)
+
+ model = struct.unpack('i',f.read(4))[0] #model
+ nrun = struct.unpack('i',f.read(4))[0] #nrun
+ nk = struct.unpack('i',f.read(4))[0] #nk
+
+ blocksize2 = struct.unpack('i',f.read(4))[0] #end header block size
+ if extra == 1:
+ f.read(4)
+
+ # check for consistency
+ if blocksize1 != blocksize2:
+ print 'header trouble! t = '+str(itime)+' '+str(blocksize1)+' '+str(blocksize2)
+ sys.exit(1)
+
+ # now read the star data
+ blocksize1 = struct.unpack('i',f.read(4))[0] #begin data block size
+ if extra == 1:
+ f.read(4)
+
+ alist = []
+ for i in range(nk):
+ alist.append(struct.unpack('f',f.read(4))[0]) #Sverre's 'as'
+ alist.append(ntot)
+
+ stardata = np.zeros((ntot, 11))
+ datalist=[] # masses
+ for i in range(ntot):
+ datalist.append(struct.unpack('f',f.read(4))[0])
+
+ stardata[:,0] = datalist
+
+ datalistx=[] # positions
+ datalisty=[]
+ datalistz=[]
+ for i in range(ntot):
+ datalistx.append(struct.unpack('f',f.read(4))[0])
+ datalisty.append(struct.unpack('f',f.read(4))[0])
+ datalistz.append(struct.unpack('f',f.read(4))[0])
+
+ stardata[:,1] = datalistx
+ stardata[:,2] = datalisty
+ stardata[:,3] = datalistz
+
+ datalistx=[] # velocities
+ datalisty=[]
+ datalistz=[]
+
+ for i in range(ntot):
+ datalistx.append(struct.unpack('f',f.read(4))[0])
+ datalisty.append(struct.unpack('f',f.read(4))[0])
+ datalistz.append(struct.unpack('f',f.read(4))[0])
+
+ stardata[:,4] = datalistx
+ stardata[:,5] = datalisty
+ stardata[:,6] = datalistz
+
+ datalist=[] # names
+ for i in range(ntot):
+ datalist.append(struct.unpack('i',f.read(4))[0])
+ stardata[:,7] = datalist
+
+ if (sev):
+ datalist=[] # kstar
+ for i in range(ntot):
+ datalist.append(struct.unpack('i',f.read(4))[0])
+ stardata[:,8] = datalist
+
+ datalist=[] # lum
+ for i in range(ntot):
+ datalist.append(struct.unpack('f',f.read(4))[0])
+ stardata[:,9] = datalist
+
+ datalist=[] # teff
+ for i in range(ntot):
+ datalist.append(struct.unpack('f',f.read(4))[0])
+ stardata[:,10] = datalist
+
+ blocksize2 = struct.unpack('i',f.read(4))[0] #end data block size
+ if extra == 1:
+ f.read(4)
+
+ # check for consistency
+ if blocksize1 != blocksize2:
+ print 'star data trouble! '+str(itime)
+ sys.exit(3)
+
+ splitarray = np.hsplit(stardata, [1,4,7,8,9,10])
+
+ # Output snapshot in hdf5 file
+ h5output(f5, alist, splitarray, model-1, sev)
+
+ # next header blocksize for loop control
+ byte = f.read(4) # begin header block
+
+ # close OUT3 and other files if present
+ f.close()
+
+
+def h5output(f5, alist, splitarray, model, sev):
+ """output the nbody6 snapshot data in hdf5 format. Plots the first snapshot
+ just to make sure it's sensible looking.
+
+ Keyword arguments:
+ outputfilename -- the name of the hdf5 file to create
+ alist -- header data from the nbody6 data. this is sverre's AS list.
+ splitarray -- the data that has been parsed from the nbody6 files
+ """
+
+ # Create model group
+ modgrp = f5.create_group("%05d"%model)
+ f5.attrs["nmod"]=model+1 # Number of models, updated everytime a snapshot is written
+
+ # create groups for the "stars" and "cluster" data
+ stargrp = f5.create_group("%05d/stars"%model)
+ clustergrp = f5.create_group("%05d/cluster"%model)
+
+ # Compute bound particles
+ nbound,mbound,vcom,bound,kinbound,potbound = compute_energies(splitarray[0][:,0],splitarray[1],splitarray[2], np.array(alist[6:9]),alist[1])
+
+
+ # Assume the following for AS list:
+ # AS(1) = TTOT
+ # AS(2) = NBIN
+ # AS(3) = RBAR
+ # AS(4) = ZMBAR
+ # AS(5) = TSTAR
+ # AS(6) = VSTAR
+ # AS(7:9) = RDENS
+ # AS(10:12) = RG
+ # AS(13:15) = VG
+ # AS(16) = RC
+ # AS(17) = NC
+
+
+ # Populate the groups with data
+ massdset = stargrp.create_dataset('mass', data = alist[3]*splitarray[0][:,0],dtype='float32')
+ posdset = stargrp.create_dataset('pos', data = alist[2]*splitarray[1],dtype='float32')
+ veldset = stargrp.create_dataset('vel', data = alist[5]*splitarray[2],dtype='float32')
+ namedset = stargrp.create_dataset('id', data = splitarray[3][:,0],dtype='i')
+ if (sev):
+ kstardset = stargrp.create_dataset('kstar', data = splitarray[4][:,0],dtype='i')
+ lumdset = stargrp.create_dataset('lum', data = splitarray[5][:,0],dtype='float32')
+ teffdset = stargrp.create_dataset('teff', data = splitarray[6][:,0],dtype='float32')
+ bounddset = stargrp.create_dataset('bound', data = bound,dtype='bool')
+
+ timedset = clustergrp.create_dataset('age', data = alist[0]*alist[4])
+ nbindset = clustergrp.create_dataset('nbin', data = alist[1],dtype='i')
+ rdensdset = clustergrp.create_dataset('rdens',data = np.array(alist[6:9])*alist[2])
+ vcomdset = clustergrp.create_dataset('vcom' ,data = np.array(vcom)*alist[5])
+ rgdset = clustergrp.create_dataset('rg', data = np.array(alist[9:12]))
+ vgdset = clustergrp.create_dataset('vg', data = np.array(alist[12:15]))
+ rcdset = clustergrp.create_dataset('rcore',data = alist[2]*alist[15])
+
+ # Compute half-mass radii
+ rh = compute_rh(alist[3]*splitarray[0][:,0],alist[2]*splitarray[1], np.array(alist[6:9])*alist[2],splitarray[4][:,0],bound, sev)
+ rhndset = clustergrp.create_dataset('rhn', data = rh[0],dtype='float32')
+ rhn300dset = clustergrp.create_dataset('rhn300',data = rh[1],dtype='float32')
+ rhmdset = clustergrp.create_dataset('rhm', data = rh[2],dtype='float32')
+
+ ncdset = clustergrp.create_dataset('ncore', data = alist[16],dtype='i')
+
+ nbdset = clustergrp.create_dataset('nbound',data = nbound,dtype='i')
+ mbdset = clustergrp.create_dataset('mbound',data = mbound*alist[3],dtype='float32')
+ kinbdset = clustergrp.create_dataset('kinbound',data = kinbound,dtype='float32')
+ potbdset = clustergrp.create_dataset('potbound',data = potbound,dtype='float32')
+
+ print " T = %10.3f [NBODY] %10.3f [MYR]; M(bound) = %10.3f, rhn = %10.3f; rhm = %10.3f"%(alist[0],alist[0]*alist[4],mbound, rh[0], rh[2])
+
+
+
+ # Units information as attributes
+ massdset.attrs["unit"] = "[Msun]"
+ posdset.attrs["unit"] = "[pc] in guide centre reference frame"
+ veldset.attrs["unit"] = "[km/s] in guide centre reference frame"
+
+ if (sev):
+ lumdset.attrs["unit"] = "[log10 Lsun]"
+ teffdset.attrs["unit"] = "[log10 Teff]"
+
+ timedset.attrs["unit"] = "[Myr]"
+ rdensdset.attrs["unit"] = "[pc] in guide centre reference frame"
+ vcomdset.attrs["unit"] = "[km/s] in guide centre reference frame"
+
+ rgdset.attrs["unit"] = "[kpc] in Galactic reference frame"
+ vgdset.attrs["unit"] = "[km/s] in Galactic reference frame"
+
+ rcdset.attrs["unit"] = "[pc]"
+ mbdset.attrs["unit"] = "[Msun]"
+
+ return()
+
+
+def main():
+ parser = argparse.ArgumentParser(description='Converts OUT3 data from NBODY6 to HDF5 format')
+ parser.add_argument('-f',dest='files',action="store",default='OUT3',help='NBODY6 files to convert [default="OUT3"]')
+ parser.add_argument('-s', dest='sev',action="store_true", default=False,help='activate reading stellar evolution data [default=False]')
+
+ args = parser.parse_args()
+
+ # open up the hdf5 file!
+ f5 = h5py.File('out.hdf5','w')
+
+ # Loop over all NBODY6 OUT3 files
+ for file in args.files.split():
+ print " File = ",file
+ read_nbody6(f5, file, args.sev)
+
+ f5.close()
+
+
+if __name__ == '__main__':
+ main()
\ No newline at end of file