program streamlines
      implicit none

! our grid
      integer, parameter :: gridlen=20 ! number of bins in x,y,z
      real, dimension(  gridlen,gridlen,gridlen) :: gridmass
      real, dimension(3,gridlen,gridlen,gridlen) :: gridvel

! data from snapshot file
      character(len=4) :: label
      integer, dimension(6) :: nparttab
      integer :: blocksize,npart1
      double precision :: redshift,time
      double precision, dimension(6) :: masstab
      real, dimension(:,:), allocatable :: pos,vel
      real, dimension(:),   allocatable :: mass

! other variables
      integer :: i,ix,iy,iz,it,jx,jy,jz
      real    :: boxsize,x,y,z,dt,g,xn,yn,zn

! read positions and velocities from snapshot file
      write(0,'("Reading snapshot file...")')
      open(1,file='snap_144',status='old',form='unformatted')
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      read(1) nparttab,masstab,redshift,time
      write(0,'("  total:  ",i8," particles")') sum(nparttab)
      write(0,'("  type ",i1,": ",i8," particles")') (i,nparttab(i),i=1,6)
! we are only interested in particle type 1 (gas particles) here
      npart1  = nparttab(1)
      allocate(pos(3,npart1),vel(3,npart1),mass(npart1))
    1 continue
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      if(label.eq.'POS ') then
        read(1) pos
      elseif(label.eq.'VEL ') then
        read(1) vel
      elseif(label.eq.'MASS') then
        read(1) mass
      else
        read(1) ! just skip the blocks we're not interested in
      endif
      goto 1
    2 continue
      close(1)
    7 format(a4,": ",i8," bytes")

! scale particle positions and velocities to our grid length
      boxsize = 48000 ! maximum coordinate value in snapshot file
      pos = pos*gridlen/boxsize
      vel = vel*gridlen/boxsize

! output particle data scaled to our grid (i.e., with coordinates 0...gridlen)
      write(0,'("Writing particle data...")')
      open(1,file='particles.txt')
      write(1,'(i8,6es12.4)') (i,pos(:,i),vel(:,i),i=1,npart1)
      close(1)

      write(0,'("Computing grid data...")')
! initialize grid to zero
      gridmass = 0.
      gridvel  = 0.
! loop over all particles and update the grid bins they fall into
      do i=1,npart1
        ix = floor(pos(1,i))+1 ! bin number in x-direction
        iy = floor(pos(2,i))+1 ! bin number in y-direction
        iz = floor(pos(3,i))+1 ! bin number in z-direction
        gridmass(  ix,iy,iz) = gridmass(  ix,iy,iz) + mass(i)          ! mass in bin
        gridvel (:,ix,iy,iz) = gridvel (:,ix,iy,iz) + mass(i)*vel(:,i) ! weighted sum of velocities
      enddo
! for all bins: divide the summed-up velocities of the bin by the number
! of particles in the bin to obtain the average velocity for the bin
      where(gridmass.ne.0.) ! avoid division by zero
        gridvel(1,:,:,:) = gridvel(1,:,:,:)/gridmass(:,:,:) ! x components
        gridvel(2,:,:,:) = gridvel(2,:,:,:)/gridmass(:,:,:) ! y components
        gridvel(3,:,:,:) = gridvel(3,:,:,:)/gridmass(:,:,:) ! z components
      endwhere

! output velocity grid to file
      write(0,'("Writing grid data...")')
      open(1,file='velocitygrid.txt')
      do iz=1,gridlen
        do iy=1,gridlen
          do ix=1,gridlen
            write(1,'(3i4,4es12.4)') ix,iy,iz,&
              gridvel(:,ix,iy,iz), gridmass(ix,iy,iz)
          enddo
        enddo
        write(1,'(/)')
      enddo
      close(1)

! calculate some streamlines and write them to file
      write(0,'("Computing and writing streamlines...")')
      dt = 1.
      g = real(gridlen)
      open(1,file='streamlines.txt')
! use (multiples of) the grid bin positions as starting points for the streamlines
      do jz=1,gridlen,1
        do jy=1,gridlen,1
          do jx=1,gridlen,1
! take the cell center as starting position
            x = jx-0.5
            y = jy-0.5
            z = jz-0.5
! write starting position to file
            write(1,'(i6,3f8.4)') 0,x,y,z
! do steps to create the streamlines
            do it=1,100
! find grid cell we're in (to get the current velocity)
              ix = floor(x)+1
              iy = floor(y)+1
              iz = floor(z)+1
! update our position by following the current velocity for a small distance
              x = x + dt*gridvel(1,ix,iy,iz)
              y = y + dt*gridvel(2,ix,iy,iz)
              z = z + dt*gridvel(3,ix,iy,iz)
! if we have left the box, wrap around to the other side
              xn = x-g*floor(x/g)
              yn = y-g*floor(y/g)
              zn = z-g*floor(z/g)
! if we have wrapped, output an empty line (this is for gnuplot)
              if(xn.ne.x .or. yn.ne.y .or. zn.ne.z) write(1,'()')
              x = xn
              y = yn
              z = zn
! write current position to file
              write(1,'(i6,3f8.4)') it,x,y,z
            enddo
! output two empty lines to mark the end of the streamline
            write(1,'(/)')
          enddo
        enddo
      enddo
      close(1)

      write(0,'("Done.")')
      end