program meandensity
      implicit none

! data from snapshot file
      character(len=4) :: label
      integer, dimension(6) :: nparttab,nparttottab,nparttothwtab
      integer :: blocksize,npart,num_files,lpt_scalingfactor
      integer :: flag_sfr,flag_feedback,flag_cooling,flag_stellarage, &
                 flag_metals,flag_entropy_instead_u, &
                 flag_doubleprecision,flag_ic_info
      double precision :: time,redshift,boxsize, &
                          omega0,omegalambda,hubbleparam
      double precision, dimension(6) :: masstab
      real, dimension(:), allocatable :: mass,rho
      real, dimension(:,:), allocatable :: pos

! other variables
      real :: rhoini,rhosim,rhovar
      integer :: i
      logical :: havemass

! read data from snapshot file
      write(0,'("Reading snapshot file...")')
      open(1,file='output_000',status='old',form='unformatted')
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      read(1) nparttab,masstab,time,redshift,flag_sfr,flag_feedback, &
              nparttottab,flag_cooling,num_files,boxsize, &
              omega0,omegalambda,hubbleparam,flag_stellarage, &
              flag_metals,nparttothwtab,flag_entropy_instead_u, &
              flag_doubleprecision,flag_ic_info,lpt_scalingfactor
      write(0,'("  box size: ",g10.3)') boxsize
      write(0,'("  total:  ",i8," particles")') sum(nparttab)
      write(0,'("  type ",i1,": ",i8," particles, mass",es12.4)') &
                           (i,nparttab(i),masstab(i),i=1,6)
      npart = nparttab(1) ! we use particle type 1 here (gas particles)
      allocate(mass(npart),rho(npart),pos(3,npart))
      mass = 0.
      rho  = 0.
      havemass = .false.
    1 continue
      read(1,end=2) label,blocksize
      write(0,7) label,blocksize-8
      if(label.eq.'MASS') then
        read(1) mass
        havemass = .true.
      elseif(label.eq.'RHO ') then
        read(1) rho
      elseif(label.eq.'POS ') then
        read(1) pos
      else
        read(1) ! just skip the blocks we're not interested in
      endif
      goto 1
    2 continue
      close(1)
    7 format(a4,": ",i8," bytes")

      if(.not.havemass) then
        write(0,'("no individual masses present")')
        mass = masstab(1)
      endif

! compute mean density and variance
      write(0,'("Computing mean density...")')
      rhoini = sum(mass)/boxsize**3
      rhosim = sum(rho)/npart
      rhovar = sum((rho-rhosim)**2)/npart
      write(*,'(4es12.4)') rhoini,rhosim,rhovar,rhosim/rhoini

! dump masses and densities to text file for plotting
      write(0,'("Writing positions and densities to file...")')
      open(1,file='density.txt')
      write(1,'(i5,7es12.4)') &
        (i,mass(i),rho(i),rho(i)/rhosim,rho(i)/rhoini,pos(:,i),i=1,npart)
      close(1)

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