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=npgal
      integer, dimension(4), parameter :: idx=[n1,n2,n3,n4]
      real,    dimension(3,npart) :: pos,vel
      real,    dimension(  npart) :: ein,mass
      integer, dimension(  npart) :: id,type
      double precision :: masstab(6)
      integer :: i,j,k,l,m,n

! read the base galaxy
      open(1,file='galaxy.txt')
      read(1,*)
      read(1,*) (id(i),mass(i),ein(i),pos(:,i),vel(:,i),i=1,npgal)
      close(1)

! set up the type array for nicer plotting
      m = 1
      do l=1,4
        n = idx(l)
        type(m:m+n-1) = l
        m = m+n
      enddo

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

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

      end