implicit none
c     Hydra simulation header file, should not touch this.
      include 'pinfo.inc'
c     Parameter file, this is what you want to modify!
      include 'pars.inc'
c     Particle data arrays, don't touch (cannot include this in pars.inc).
      include 'particledata.inc'
c
      real sphgrid(n1d,n1d,n1d)
      real tmass, CurrentHubbleCt, box4kms
      integer iselect
      
      CALL readdata_tt(simdatafile,.true.,nobj,rm,r
     &     ,v,hi,hei,heiii,h,dn,e,zmetal,itype)
      CALL readdata_tt(simdatafile,.false.,nobj,rm,r
     &     ,v,hi,hei,heiii,h,dn,e,zmetal,itype)

      write (*,*) ' data read-in '
      call flush(6)
c     
      CurrentHubbleCt = 100. * h100 * 
     &     sqrt(1. + omega0*(1./atime-1.) + xlambda0*
     &     (atime**2-1.)) /atime
      boxkms = box100 / h100 * atime * CurrentHubbleCt
      box4kms = boxkms

c     compute density on grid using TSC interpolation
      iselect = 1
      call tsc_density(iselect,tmass,sphgrid)
      write (*,*) ' tsc density: total mass = ',tmass
      open(unit=1,file='tscgrid_bar.d',status='unknown',form='unformatted')
      write (1) n1d
      write (1) tmass,box100,h100,atime,box4kms
      write (1) sphgrid
      close (1)
      write (*,*) ' step 1 '
      call flush(6)
c
      iselect = 2
      call tsc_density(iselect,tmass,sphgrid)
      write (*,*) ' tsc density: total mass = ',tmass
      open(unit=1,file='tscgrid_gas.d',status='unknown',form='unformatted')
      write (1) n1d
      write (1) tmass,box100,h100,atime,box4kms
      write (1) sphgrid
      close (1)
      write (*,*) ' step 2 '
      call flush(6)

c     
      iselect = 3
      call tsc_density(iselect,tmass,sphgrid)
      write (*,*) ' tsc density: total mass = ',tmass
      open(unit=1,file='tscgrid_dark.d',status='unknown',form='unformatted')
      write (1) n1d
      write (1) tmass,box100,h100,atime,box4kms
      write (1) sphgrid
      close (1)
      write (*,*) ' step 3 '
      call flush(6)

c      stop 'tsc_density written'

c     compute density on grid using SPH interpolation
      call sph_density(tmass,sphgrid)
      write (*,*) ' sph density: total mass = ',tmass
      open(unit=1,file='sphgrid.d',status='unknown',form='unformatted')
      write (1) n1d
      write (1) tmass,box100,h100,atime,box4kms
      write (1) sphgrid
      close (1)


      end

      subroutine tsc_density(iselect,tmass,sphgrid)
      implicit none
      include 'pinfo.inc'
      include 'pars.inc'
      include 'particledata.inc'
c
      real    sphgrid(n1d,n1d,n1d),rmp
c     compute density on rectangular grid using SPH interpolation
      real tmass,dx,rdx,box,box2,hh,h2,h4,hinv2,hinv3,hinv4,xx,yy,zz,xgrid,ygrid,zgrid,dx,dy,dz,dr2,q,kernel
      real*8 tmass8
      real dxgrid,rdxgrid,norm,DensCon
      integer i,j,k,ii,jj,kk,ipart,ioff,ix,iy,iz
      integer imass
      integer nper
      integer iselect

c     initialize grid
      do k=1,n1d
         do j=1,n1d
            do i=1,n1d
               sphgrid(i,j,k)=0.0
            enddo
         enddo
      enddo
      tmass = 0.0
c     
      dxgrid   = 1./float(n1d)
      rdxgrid  = 1./dxgrid
      box      = 1.0
      box2     = 0.5*box
      tmass8   = tmass
      nper     = 0
c
      do ipart=1,nobj
         if(iselect .eq. 1)then
            if(itype(ipart) .eq. itdark) goto 200 ! stars and gas
         else if (iselect .eq. 2) then
            if(itype(ipart) .ne. itgas) goto 200 ! gas only
         else
            if(itype(ipart) .ne. itdark) goto 200 ! dark matter only
         endif
c
         if(float(ipart) .ge. float(nper)*float(nobj)/100.) then
            write (*,*) ' percent done = ',nper
            nper = nper + 5
            call flush(6)
         endif
c     
         rmp    = rm(ipart)
         tmass8 = tmass8 + rmp
         imass  = imass + 1
c     
         xx = r(1,ipart)
         ix = int(xx*rdxgrid)
         dx = xx*rdxgrid-float(ix)
         ix = ix+1
         ii = mod(ix,n1d)+1
c     
         yy = r(2,ipart)
         iy = int(yy*rdxgrid)
         dy = yy*rdxgrid-float(iy)
         iy = iy+1
         jj = mod(iy,n1d)+1
c     
         zz = r(3,ipart)
         iz = int(zz*rdxgrid)
         dz = zz*rdxgrid-float(iz)
         iz = iz+1
         kk = mod(iz,n1d)+1
c     
         
         sphgrid(ix,iy,iz) = sphgrid(ix,iy,iz)+(1.-dx)*(1.-dy)*(1.-dz)*rmp
c     
         sphgrid(ii,iy,iz) = sphgrid(ii,iy,iz)+    dx *(1.-dy)*(1.-dz)*rmp
         sphgrid(ix,jj,iz) = sphgrid(ix,jj,iz)+(1.-dx)*    dy *(1.-dz)*rmp
         sphgrid(ix,iy,kk) = sphgrid(ix,iy,kk)+(1.-dx)*(1.-dy)    *dz *rmp
c     
         sphgrid(ii,jj,iz) = sphgrid(ii,jj,iz)+    dx *    dy *(1.-dz)*rmp
         sphgrid(ii,iy,kk) = sphgrid(ii,iy,kk)+    dx *(1.-dy)    *dz *rmp
         sphgrid(ix,jj,kk) = sphgrid(ix,jj,kk)+(1.-dx)    *dy     *dz *rmp
c     
         sphgrid(ii,jj,kk) = sphgrid(ii,jj,kk)+    dx     *dy     *dz *rmp
c     
 200     continue               ! particle of right type
      enddo                     ! i=1,nobj
c     normalise density to mean of 1.
      tmass  = tmass8
c$$$  tmass8 = 0.0
c$$$  do k=1,n1d
c$$$         do j=1,n1d
c$$$            do i=1,n1d
c$$$               tmass8 = tmass8 + sphgrid(i,j,k)
c$$$            enddo
c$$$         enddo
c$$$      enddo

c     convert to equivalent Hydrogen atoms/cm^3.
      DensCon = 10.**(alog10(munit) -3.*alog10(lunit)-alog10(m_H))*float(n1d)**3

      do k=1,n1d
         do j=1,n1d
            do i=1,n1d
               sphgrid(i,j,k)=sphgrid(i,j,k)*DensCon
            enddo
         enddo
      enddo

      return
      end

      subroutine sph_density(tmass,sphgrid)
      implicit none
      include 'pinfo.inc'
      include 'pars.inc'
      include 'particledata.inc'
c
      real    sphgrid(n1d,n1d,n1d)
      real*8 tmass8
c     compute density on rectangular grid using SPH interpolation
      real tmass,dx,rdx,box,box2,hh,h2,h4,hinv2,hinv3,hinv4,xx,yy,zz,xgrid,ygrid,zgrid,dx,dy,dz,dr2,q,kernel
      real dxgrid,rdxgrid,norm,DensCon
      integer i,j,k,ii,jj,kk,ipart,ioff,ix,iy,iz
      integer nper

c     initialize grid
      do k=1,n1d
         do j=1,n1d
            do i=1,n1d
               sphgrid(i,j,k)=0.0
            enddo
         enddo
      enddo
      tmass8 = 0.0
c     
      dxgrid   = 1./float(n1d)
      rdxgrid  = 1./dxgrid
      box      = 1.0
      box2     = 0.5*box
      nper     = 0
      do ipart=1,nobj
         if(itype(ipart) .eq. itgas)then
            if(float(ipart) .ge. float(nper)*float(nobj)/100.) then
               write (*,*) ' percent done = ',nper
               nper = nper + 5
               call flush(6)
            endif
            tmass8 = tmass8 + rm(ipart)
c
            hh    = h(ipart)
            h2    = hh*hh
            h4    = 4.*h2
            hinv2 = 1. / h2
            hinv3 = hinv2 / hh
c
            xx = r(1,ipart)
            ix = int(xx*rdxgrid)+1
            yy = r(2,ipart)
            iy = int(yy*rdxgrid)+1
            zz = r(3,ipart)
            iz = int(zz*rdxgrid)+1
            ioff = int(hh*2.*rdxgrid)+2
c     loop over all contributing vertices
            do k=iz-ioff,iz+ioff
               kk = k
               if (kk .gt. n1d) kk=kk-n1d
               if (kk .le.   0) kk=kk+n1d
               zgrid = float(kk-1)*dxgrid
               do j=iy-ioff,iy+ioff
                  jj = j
                  if (jj .gt. n1d) jj=jj-n1d
                  if (jj .le.   0) jj=jj+n1d
                  ygrid = float(jj-1)*dxgrid

                  do i=ix-ioff,ix+ioff
                     ii = i
                     if (ii .gt. n1d) ii=ii-n1d
                     if (ii .le.   0) ii=ii+n1d
                     xgrid = float(ii-1)*dxgrid
c     test whether vertex is close enough
                     dx = abs(xx-xgrid)
                     if(dx .gt. box2) dx=box-dx
                     dr2 = dx*dx
                     if(dr2 .lt. h4)then
                        dy = abs(yy-ygrid)
                        if(dy .gt. box2) dy=box-dy
                        dr2 = dr2+dy*dy
                        if(dr2 .le. h4)then
                           dz = abs(zz-zgrid)
                           if(dz .gt. box2) dz=box-dz
                           dr2 = dr2+dz*dz
                           if(dr2 .le. h4)then                           
                              q = sqrt(dr2 * hinv2)
                              if (q .le. 1.) then
                                 kernel = (1.+q**2*(-1.5+0.75*q))/pi
                              else
                                 kernel = 0.25*(2.-q)**3/pi
                              endif
                              kernel = kernel * hinv3
                              kernel = kernel * rm(ipart)
                              sphgrid(ii,jj,kk) = sphgrid(ii,jj,kk) + kernel
                           endif ! dr2< h4
                        endif   ! dr2 < h4
                     endif      ! dx^2<h4
                  enddo         ! ixoff
               enddo            ! iyoff
            enddo               ! izoff
         endif                  ! itype = itgas
      enddo                     ! i=1,nobj

c     convert to equivalent Hydrogen atoms/cm^3. Currently in code units
      DensCon = 10.**(alog10(munit) -3.*alog10(lunit)-alog10(m_H)+alog10(1.-ymass))
      
c     normalise density to mean of 1.
      tmass = tmass8
      norm  = 1./tmass
      do k=1,n1d
         do j=1,n1d
            do i=1,n1d
               sphgrid(i,j,k)=sphgrid(i,j,k)*DensCon
            enddo
         enddo
      enddo

      return
      end

      SUBROUTINE readdata_tt(filename,honly,N,rm,r,v,hi,hei,heiii,h,dn,e
     &     ,zmetal,itype)

c     Read in data of Tom Theuns' simulation (modified Hydra).

      implicit none

      include 'pinfo.inc'
#include "pars.inc"

      integer N

      real*4 rm(N),r(3,N),v(3,N),dn(N),e(N),h(N)
      real*4 hi(N),hei(N),heiii(N),zmetal(N)
      integer*4 itype(N)


      real*8 zcurrent
      real*8 rmtot
      real*8 meanmol
      character*(*) filename
      character*72 longfile
      character dumbo
      logical honly             ! Read header only?

      integer lnblnk, ier, i, ic, iczfill
      real*8 fdark, fgas, fstar, ftot, rhoc, rscale, vscale, tempscale
      real*8 dscale, tmin, tmax, dnmin, dnmax, zmax, tmin_ib, tmax_ib
      real*8 dmin_ib, dmax_ib, OverDen, heliumiii, hydrogeni, hydrogenii
      real*8 heliumi, heliumii, enumber
      real*8 rlow, rhigh, ran3
      integer idum
c
      logical testout           ! write data to ascii file
      parameter (testout=.false.)
c

      longfile = datadir 
      longfile = longfile(1:lnblnk(longfile))//filename
      open(3,file=longfile,status='old',form='unformatted',iostat=ier)
      if (ier .ne. 0) then
         write(*,*) 'ERROR in readdata_tt: Problem reading file ', longfile
         stop
      endif 

      read(3) ibuf,ibuf1,ibuf2

      if (nobj .gt. Nmax) then
         write (*,*) 'ERROR in readdata_tt: N>Nmax'
         write (*,'('' N= '',i12,'' nmax= '',i12)') N,Nmax
         stop
      endif

      atime_sim = atime
      ztime = 1. / atime - 1.

      if (honly) then           ! Read header only
         close(3)
         return
      endif 

      read(3) rm
      read(3) r
      read(3) v
      read(3) h
      read(3) e
      read(3) itype
      read(3) dn
      read(3) hi
      read(3) hei
      read(3) heiii

      if (Z_rel .lt. 0.) then
#ifndef SIGHT
         read(3)
         read(3)
#endif
         read(3) zmetal
      endif

      close(3)

c     determine baryon fraction
      if (omegab0 .eq. 0.0) then  ! if omegab0 not in header.
         if (verbose) write(*,*) 'Computing omegab0...'
         fdark = 0.0
         fgas  = 0.0
         fstar = 0.0
         do i = 1, nobj
            if (itype(i) .eq. itdark) then
               fdark = fdark + rm(i)
            elseif (itype(i) .eq. itgas) then
               fgas  = fgas + rm(i)
            elseif (itype(i) .eq. itstar) then
               fstar = fstar + rm(i)
            endif
         enddo
         ftot = fdark + fgas + fstar
         omegab0 = (fgas + fstar) / ftot * omega0
         if (verbose) write(*,*) 'done.'
      endif
      if (omegab0 .lt. 1e-20) then 
         omegab0 = 0.019 / h100**2
         write(*,'("WARNING: omegab0 in header must be wrong!")')
         write(*,'("WARNING: assuming omegab0 = 0.019 / h^2 = ",f5.3)') 
     &        omegab0
      endif 

      if (verbose .and. .not. sightlines) then 
         write(*,*)
         write(*,*) 'Simulation parameters (set in code):'
         write(*,'("omega_m      = ",f5.3)') omega0
         write(*,'("omega_lambda = ",f5.3)') xlambda0
         write(*,'("omega_b*h^2  = ",f5.3)') omegab0 * h100**2
         write(*,'("h            = ",f5.3)') h100
         write(*,'("box size     = ",f8.3," h^{-1} Mpc, comoving")')
     &        box100
         write(*,'("z            = ",f5.3)') ztime
         write(*,*)
      endif 


c     ---------------------------
c     Compute conversion factors:
c     ---------------------------

c     Code length unit is size of box.
      lunit = box100 * Mpc / h100 * atime ! cm, physical
c     Code time unit is present age of the universe.
      tunit = h0t0 / H0 / h100  ! s
      rmtot = 0.
      do i = 1, nobj
         if (itype(i) .gt. itnone) rmtot = rmtot + rm(i)
      enddo
c
      rhoc = 3 * (H0*h100)**2 / (8. * pi * G) ! g/cm^3
      munit = log10(rhoc / atime**3 * omega0) + 3.*log10(lunit) -
     &     log10(rmtot)         ! g
      munit = 10.**munit

      vunit = lunit / tunit     ! cm/s, physical
      dunit = log10(munit) - 3.*log10(lunit)  ! g/cm^3 
      dunit = 10.**dunit
c     Coordinates to physical Mpc
      rscale = lunit / Mpc      ! Mpc, physical
c     Peculiar velocity to physical km/s
      vscale = vunit / 1.e5     ! km/s
c     Temperature to K / mu (mu = mean particle mass in units of m_H)
      tempscale = m_H * vunit**2 * 2. / 3. / Boltz ! K/mu
c     Density to n_H (cm^-3)
      dscale = dunit / m_H * (1.-ymass) ! n_H/cm^3
c     Mean baryon density (n_H/cm^3)
      rhocb = rhoc / atime**3 / m_H * (1.-ymass) * omegab0
c
      return
      end