pro setup_evcoll
; Setup unit system
  M_star=1.5815
  R_star=1.0
  G=1.0

  fglass='glass_100x100x100'
  readnew,fglass,h,'HEAD'
  readnew,fglass,x,'POS'

  N=20000L
  fout='evcoll_ini_20k'
  v=fltarr(3,N)
  m=fltarr(N)
  id=LINDGEN(N)+1
  u=fltarr(N)

  ; Check if a sphere with N particle fits into the grid
  IF N GT 4 * !Pi / 3 * TOTAL(h.npart,/int) / 8 THEN stop

  ; Ensure that all is within [-1,1]
  x=x/h.BoxSize
  x=x-0.5
  x=x*2

  ; Compute radii for each particle
  rr=sqrt(x(0,*)^2+x(1,*)^2+x(2,*)^2)

  ; To get a 1/r density profile, we need to multiply by sqrt(r)
  FOR i=0,2 DO x(i,*)=x(i,*)*sqrt(rr)

  ; Now we want to have N particles
  rr=sqrt(x(0,*)^2+x(1,*)^2+x(2,*)^2)
  ii=SORT(rr)

  x=x(*,ii(0:N-1))
  rr=sqrt(x(0,*)^2+x(1,*)^2+x(2,*)^2)

  ; Scale to the units of this test
  x=FLOAT(x/max(rr)*R_star)
  m(*)=M_star/N
  u(*)=G*M_star/R_star*0.05

  h.npart(*)=0
  h.parttotal(*)=0
  h.npart(0)=N
  h.parttotal(0)=N

  write_head,fout,h
  add_block,fout,x,'POS'
  add_block,fout,v,'VEL'
  add_block,fout,id,'ID'
  add_block,fout,m,'MASS'
  add_block,fout,u,'U'

end