program jetsetup
      implicit none

! maximum number of particles for container, pipe, water, air, syrup, and plunger
      integer, parameter :: ncx=10, ncy=1,  ncz=1
      integer, parameter :: npx=2,  npy=17, npz=1
      integer, parameter :: nwx=10, nwy=9,  nwz=1
      integer, parameter :: nax=5,  nay=10, naz=1 ! (not all are used)
      integer, parameter :: nsx=2,  nsy=32, nsz=2
      integer, parameter :: nux=1,  nuy=1,  nuz=1
      integer, parameter :: nbase=1000
      integer, parameter :: npmax=nbase*(ncx*ncy*ncz + npx*npy*npz + &
        nwx*nwy*nwz + nax*nay*naz + nsx*nsy*nsz + nux*nuy*nuz)
      real,    dimension(3,nbase) :: grid
      real,    dimension(3,npmax) :: pos,vel
      real,    dimension(  npmax) :: ein,mass
      integer, dimension(  npmax) :: id
      double precision, dimension(6) :: masstab
      integer :: i,j,k,l,ntot,nfluid
      real :: gamma,cs,rho0,eps0,m0,v0,x,y,z

! static properties of medium
      gamma = 5./3.
      mass  = 1.
      rho0  = 1.
      eps0  = 1.
      m0    = rho0/nbase
      v0    = -0.1

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

      ntot   = 0
      nfluid = 0
! set up the container (bottom wall)
      do i=1,ncx
        do j=1,ncy
          do k=1,ncz
            pos(1,ntot+1:ntot+nbase) = grid(1,:) + (i-1)
            pos(2,ntot+1:ntot+nbase) = grid(2,:) + (j-1)
            pos(3,ntot+1:ntot+nbase) = grid(3,:) + (k-1)
            vel(:,ntot+1:ntot+nbase) = 0.
            ein(  ntot+1:ntot+nbase) = eps0
            mass( ntot+1:ntot+nbase) = m0*2.
            id (  ntot+1:ntot+nbase) = 0
            ntot = ntot+nbase
          enddo
        enddo
      enddo

! set up the pipe
      do i=1,npx
        do j=1,npy
          do k=1,npz
            pos(1,ntot+1:ntot+nbase) = grid(1,:) + (i-1)*2. + 3.5
            pos(2,ntot+1:ntot+nbase) = grid(2,:) + (j-1)    + 12.
            pos(3,ntot+1:ntot+nbase) = grid(3,:) + (k-1)
            vel(:,ntot+1:ntot+nbase) = 0.
            ein(  ntot+1:ntot+nbase) = eps0
            mass( ntot+1:ntot+nbase) = m0*2.
            id (  ntot+1:ntot+nbase) = 0
            ntot = ntot+nbase
          enddo
        enddo
      enddo

! set up the plunger
      do i=1,nux
        do j=1,nuy
          do k=1,nuz
            pos(1,ntot+1:ntot+nbase) = grid(1,:) + (i-1) + 4.5
            pos(2,ntot+1:ntot+nbase) = grid(2,:) + (j-1) + 28.
            pos(3,ntot+1:ntot+nbase) = grid(3,:) + (k-1)
            vel(:,ntot+1:ntot+nbase) = 0.
            vel(2,ntot+1:ntot+nbase) = v0
            ein(  ntot+1:ntot+nbase) = eps0
            mass( ntot+1:ntot+nbase) = m0*2.
            id (  ntot+1:ntot+nbase) = 0
            ntot = ntot+nbase
          enddo
        enddo
      enddo

! set up the water
      do i=1,nwx
        do j=1,nwy
          do k=1,nwz
            pos(1,ntot+1:ntot+nbase) = grid(1,:) + (i-1)
            pos(2,ntot+1:ntot+nbase) = grid(2,:) + (j-1) + 1.
            pos(3,ntot+1:ntot+nbase) = grid(3,:) + (k-1)
            vel(:,ntot+1:ntot+nbase) = 0.
            ein(  ntot+1:ntot+nbase) = eps0
            mass( ntot+1:ntot+nbase) = m0
            id (  ntot+1:ntot+nbase) = [(l,l=nfluid+1,nfluid+nbase)]
            ntot   = ntot+nbase
            nfluid = nfluid+nbase
          enddo
        enddo
      enddo

! set up the syrup
      do i=1,nsx
        do j=1,nsy
          do k=1,nsz
            pos(1,ntot+1:ntot+nbase) = grid(1,:)*0.5 + (i-1)*0.5 + 4.5
            pos(2,ntot+1:ntot+nbase) = grid(2,:)*0.5 + (j-1)*0.5 + 12.
            pos(3,ntot+1:ntot+nbase) = grid(3,:)*0.5 + (k-1)*0.5
            vel(:,ntot+1:ntot+nbase) = 0.
            vel(2,ntot+1:ntot+nbase) = v0
            ein(  ntot+1:ntot+nbase) = eps0/8.
            mass( ntot+1:ntot+nbase) = m0
            id (  ntot+1:ntot+nbase) = [(l,l=nfluid+1,nfluid+nbase)]
            ntot   = ntot+nbase
            nfluid = nfluid+nbase
          enddo
        enddo
      enddo

      write(0,'("Liquid ends at ",i0)') nfluid

! set up the air
      do i=1,nax
        do j=1,nay
          do k=1,naz
            do l=1,nbase
              x = grid(1,l)*2. + (i-1)*2.
              y = grid(2,l)*2. + (j-1)*2. + 10.
              z = grid(3,l)*2. + (k-1)*2.
              if(z.lt.1. .and. &
                (y.lt.12. .or. y.gt.29. .or. x.lt.3.5 .or. x.gt.6.5)) then
                ntot   = ntot+1
                nfluid = nfluid+1
                pos(1,ntot) = grid(1,l)*2. + (i-1)*2.
                pos(2,ntot) = grid(2,l)*2. + (j-1)*2. + 10.
                pos(3,ntot) = grid(3,l)*2. + (k-1)*2.
                vel(:,ntot) = 0.
                ein(  ntot) = eps0*8.
                mass( ntot) = m0
                id (  ntot) = nfluid
              endif
            enddo
          enddo
        enddo
      enddo

      write(0,'("Using ",i0," particles total, ",i0," fluid particles")') ntot,nfluid

! dump coordinates to text file for debugging
      open(1,file='grid2.txt')
      write(1,'(i6,6f8.4)') (id(i),pos(:,i),vel(:,i),i=1,ntot)
      close(1)

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

      end