program blastsetup
      implicit none

      real, parameter :: pi=4.*atan(1.)
! we use 10x10x10 cubes of our base grid
      integer, parameter :: nbase=1000, npart=10*10*10*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,k,l
      real :: gamma,cs,rho0,eps0,p0
      double precision :: mass
      real, dimension(3) :: pos0

! static properties of medium
      gamma = 5./3.
      mass  = 1.
      rho0  = npart*mass/2.**3 ! 1000 particles of mass 1 per cube of (0.2)^3
      eps0  = 1.e-20
      p0    = (gamma-1.)*eps0*rho0
      cs    = sqrt(gamma*p0/rho0)

! the base grid has coordinate values 0...1 in all three dimensions
      open(1,file='glass.txt')
      read(1,*) grid
      close(1)

! we scale the coordinates of our initial setup to 2x2x2
      do i=1,10
        do j=1,10
          do k=1,10
            l = (100*(i-1)+10*(j-1)+k-1)*nbase
            pos(1,l+1:l+nbase) = 0.2*(grid(1,:) + (i-1))
            pos(2,l+1:l+nbase) = 0.2*(grid(2,:) + (j-1))
            pos(3,l+1:l+nbase) = 0.2*(grid(3,:) + (k-1))
          enddo
        enddo
      enddo
      forall(i=1:npart) id(i)=i
      vel = 0.
      ein = eps0

! set up the blast
      do i=1,npart
        pos0 = pos(:,i)-[1.,1.,1.]
        if(dot_product(pos0,pos0).lt.0.04**2) ein(i) = 1.
      enddo

! dump coordinates to text file for debugging
      open(1,file='grid2.txt')
      write(1,'(i8,f6.0,3f8.4)') (id(i),ein(i),pos(:,i),i=1,npart)
      close(1)

! write the initial conditions file
      masstab = 0.
      masstab(1) = mass
      open(1,file='blast.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