program cloudsetup
      implicit none

      double precision :: masstab(6),mass
      real,    dimension(3) :: pos1
      real,    dimension(:,:), allocatable :: pos,vel
      real,    dimension(:),   allocatable :: ein,rad
      integer, dimension(:),   allocatable :: id
      integer :: npart,i
      real :: G,Rstar,Mstar,rgrid,r

      G     = 1.0
      Rstar = 1.0
      Mstar = 1.5815
! Increasing this will use more particles in the simulation:
      rgrid = 0.2 ! must be less than or equal to 0.5

! Our base grid has 1e6 particles with coordinates in the range 0...1,
! and we count how many there are within a sphere of radius rgrid,
      npart = 0
      open(1,file='glass_100x100x100.txt',status='old')
    1 read(1,*,end=2) pos1
      if(dot_product(pos1-0.5,pos1-0.5).lt.rgrid**2) npart=npart+1
      goto 1
    2 write(0,'("using ",i0," particles")') npart
! then we allocate our arrays and populate them with these particles.
      allocate(pos(3,npart),vel(3,npart),ein(npart),id(npart),rad(npart))
      i = 0
      rewind(1)
    3 read(1,*,end=4) pos1
      if(dot_product(pos1-0.5,pos1-0.5).lt.rgrid**2) then
        i = i+1
        id(i) = i
        pos(:,i) = (pos1-0.5)*Rstar/rgrid
      endif
      goto 3
    4 close(1)
! Get a 1/r distribution with maximum at Rstar by suitably scaling r.
! (We had previously defined the radii so that the maximum is at 1.)
      do i=1,npart
        r        = sqrt(dot_product(pos(:,i),pos(:,i)))
        pos(:,i) = pos(:,i)*sqrt(r)*Rstar
        rad(i)   = sqrt(dot_product(pos(:,i),pos(:,i)))
      enddo
      vel = 0.
      ein = 0.05*G*Mstar/Rstar
      mass = Mstar/npart
      masstab = 0.
      masstab(1) = mass

! dump coordinates to text file for debugging
      open(1,file='grid.txt')
      write(1,'(i6,4f8.4)') (i,rad(i),pos(:,i),i=1,npart)
      close(1)

! write the initial conditions file
      open(1,file='cloud.ic',form='unformatted')
      write(1) 'HEAD',8+256
      write(1) [npart,0,0,0,0,0],masstab, &
        (0,i=1,(256-6*4-6*8)/4) ! pad to 256 bytes
      write(1) 'ID  ',8+4*npart
      write(1) id
      write(1) 'POS ',8+3*4*npart
      write(1) pos
      write(1) 'VEL ',8+3*4*npart
      write(1) vel
      write(1) 'U   ',8+4*npart
      write(1) ein
      close(1)

      end