program slabsetup
      implicit none

      real, parameter :: pi=4.*atan(1.)
! we use 20 cubes of our base grid for the slab
      integer, parameter :: nbase=1000, npart=20*nbase
      real,    dimension(3,nbase) :: grid
      real,    dimension(3,npart) :: pos,vel
      real,    dimension(  npart) :: ein
      integer, dimension(  npart) :: id
      double precision, dimension(6) :: masstab
      integer :: i,j
      real :: gamma,cs,rho0,eps0,p0,dv,drho,deps

! static properties of medium
      gamma = 5./3.
      rho0  = 1000.
      eps0  = 1.
      p0    = (gamma-1.)*eps0*rho0
      cs    = sqrt(gamma*p0/rho0)

! amplitudes of perturbation
      dv    = 0.01
      drho  = dv*rho0/cs
      deps  = (gamma-1.)*eps0*drho/rho0

! the base grid has coordinates 0...1 in all 3 dimensions
      open(1,file='glass.txt',status='old')
      read(1,*) grid
      close(1)

! set up the unperturbed slab with 20 copies of the base grid
! we have 1000 particles of mass=1 per cube of length=1 => rho=1000
      do i=1,20
        j = (i-1)*nbase
        pos(1,j+1:j+nbase) = grid(1,:) + i-1
        pos(2,j+1:j+nbase) = grid(2,:)
        pos(3,j+1:j+nbase) = grid(3,:)
      enddo
      forall(i=1:npart) id(i)=i
      vel = 0.
      ein = eps0

! set up perturbation
! velocity, density, and internal energy are all in phase
      do i=1,npart
        vel(1,i) = vel(1,i) +   dv*sin(pi*pos(1,i))
        pos(1,i) = pos(1,i) + drho*cos(pi*pos(1,i))/(pi*rho0)
        ein(  i) = eps0     + deps*sin(pi*pos(1,i))
      enddo

! dump to file for debugging
      open(1,file='grid.txt')
      write(1,'(i6,7es12.4)') &
        (id(i),ein(i),pos(:,i),vel(:,i),i=1,npart)
      close(1)

! write the initial conditions file
      masstab = [1.,0.,0.,0.,0.,0.]
      open(1,file='slab.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