program galaxysetup
      implicit none

! number of particles in our base galaxy
      integer, parameter :: n1=15238, n2=141380, n3=60952, n4=32610
      integer, parameter :: npgal=n1+n2+n3+n4, npart=2*npgal
      integer, dimension(4), parameter :: idx=[n1,n2,n3,n4]
      real,    dimension(3,npgal) :: pos0,vel0
      real,    dimension(  npgal) :: ein0,mass0
      integer, dimension(  npgal) :: id0
      real,    dimension(3,npart) :: pos,vel
      real,    dimension(  npart) :: ein,mass
      integer, dimension(  npart) :: id
      double precision :: masstab(6)
      integer :: i,j,k,l,m,n
      real, dimension(3,3) :: rotmat

! read the base galaxy
      open(1,file='galaxy.txt')
      read(1,*)
      read(1,*) (id0(i),mass0(i),ein0(i),pos0(:,i),vel0(:,i),i=1,npgal)
      close(1)

      rotmat = reshape([ 1., 0., 0., &
                         0., 0., 1., &
                         0.,-1., 0.],[3,3])

! set up two copies of the galaxy
      j = 0
      m = 0
      do l=1,4
        n = idx(l)
        k = j+n
        do i=1,n
          id(i+j) = id0(i+m) + m
          id(i+k) = id0(i+m) + m + n
          pos(:,i+j) =               pos0(:,i+m)  + [-100.,0.,0.]
          pos(:,i+k) = matmul(rotmat,pos0(:,i+m)) + [ 100.,0.,0.]
          vel(:,i+j) =               vel0(:,i+m)  + [0.,0.,-20.]
          vel(:,i+k) = matmul(rotmat,vel0(:,i+m)) + [0.,0., 20.]
          mass(i+j) = mass0(i+m)
          mass(i+k) = mass0(i+m)
          ein(i+j) = ein0(i+m)
          ein(i+k) = ein0(i+m)
        enddo
        j = j+2*n
        m = m+n
      enddo

! dump coordinates to text file for debugging
      open(1,file='grid.txt')
      write(1,'(i8,3es12.4)') (id(i),pos(:,i),i=1,npart)
      close(1)

! write the initial conditions file
      masstab = 0.
      masstab(1) = 0.
      masstab(2) = mass0(n1+1)
      masstab(3) = mass0(n1+n2+1)
      masstab(4) = mass0(n1+n2+n3+1)
      open(1,file='galaxies.ic',form='unformatted')
      write(1) 'HEAD',8+256
      write(1) [2*n1,2*n2,2*n3,2*n4,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+4*npart*3
      write(1) pos
      write(1) 'VEL ',8+4*npart*3
      write(1) vel
      write(1) 'MASS',8+4*n1*2
      write(1) mass(1:2*n1)
      write(1) 'U   ',8+4*n1*2
      write(1) ein(1:2*n1)
      close(1)

      end