C ______________________________	 
C		             November 2003         A. Klypin (aklypin@nmsu.edu) 
C		                     
C          Compare two snapshots
C          
C
C
      INCLUDE 'PMparameters.h'
      REAL     INPUT
      DIMENSION         RECDAT2(NRECL)
      COMMON / ROW2/XPAR2(NPAGE),YPAR2(NPAGE),ZPAR2(NPAGE),
     +	VX2(NPAGE),VY2(NPAGE),VZ2(NPAGE)
      REAL*8   disAvr,disMax
      INTEGER  iDispl(0:200)
      logical inside
      Character  FileASCII*50,snap1*4,snap2*4
      EQUIVALENCE    (RECDAT2(1),XPAR2(1))
      DATA disAvr,disMax,xmin,xmax,vmin,vmax,Icount,iDispl
     &    /0.d0,0.d0,1.e+9,-1.e+9,1.e+9,-1.e+9,0,201*0/
C...................................................................
C			Read data and open files
      ifile1      =1
      ifile2      =2
      snap1 ='1657'
      snap2='1677'
      hDispl   =0.5/200
      FileASCII ='SnapCompare.'//snap1//'.dat'
      open(17,file=FileASCII )
      CALL RDTAPE(snap1,9,ifile1)
      CALL RDTAPE(snap2,10,ifile2)

      write (*,*) ' RDTAPE is done'

      WRITE (17,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
100   FORMAT(1X,'Header=>',A45,/
     +          1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +           1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/
     +           1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/
     +           1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +           1x,' Omega_curvature=',F7.3)

      If(extras(100).ne.0.)Then
         Box = extras(100)
      else
        Box =INPUT(' Enter box size in Mpc/h   =')
      endif 
      write (*,*) ' Box size is =  ',Box,'Mpc/h'
      BoxV =Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
c      ScaleC = Box/NGRID         ! scale factor for Coordinates
      ScaleC = 1.

         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
      write (*,*) ' N_in_last=',N_in_last
C				       Loop over particles
       DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            write (*,*)' Read page=',IROW,' Max=',DisMax,' N=',In_page
            iL = NPAGE*(IROW-1)
            CALL GETROW(IROW,ifile2)
              Do i =1,NRECL
                 RECDAT2(i) =RECDAT(i)
              EndDo 
            CALL GETROW(IROW,ifile1)
                DO IN=1,In_page
                  Icount =Icount +1
                  X  =XPAR(IN)
                  Y  =YPAR(IN)
                  Z  =ZPAR(IN)
                  X1 =XPAR2(IN)
                  Y1 =YPAR2(IN)
                  Z1 =ZPAR2(IN)
                  dx =min(abs(X-X1),NGRID-abs(X-X1))
                  dy =min(abs(Y-Y1),NGRID-abs(Y-Y1))
                  dz =min(abs(Z-Z1),NGRID-abs(Z-Z1))
                  if(dx.gt.1.)write(*,300)dx,X,X1,Y,Y1,Z,Z1,Icount
 300              format('X:',g12.3,2x,3(2g12.4,2x),i9)
                  displ     = sqrt(dx**2+dy**2+dz**2)
                  disAvr  =disAvr +displ**2
                  disMax =max(displ,disMax)
                  indx = MIN(INT(displ/hDispl),200)
                  iDispl(indx) = iDispl(indx)+1
                  Vxs=ScaleV* VX(IN)
                  Vys=ScaleV* VY(IN)
                  Vzs=ScaleV* VZ(IN)
                  xmax =MAX(xmax,X,Y,Z)
                  xmin   =MIN(xmin,X,Y,Z)
                  vmax =MAX(vmax,Vxs,Vys,Vzs)
                  vmin  =MIN(vmin,Vxs,Vys,Vzs)
               ENDDO
        ENDDO 
      write (*,*) ' Compare two snapshots        =',snap1,'  ',snap2 
      write (*,*) ' Number of particles          =',Icount
      write (*,*) ' Min/Max of coordinates(Mpc/h)=',xmin,xmax
      write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax
      write (*,*) ' Max displacement (0-cell)    =',disMax
      write (*,*) ' Avr displacement             =',sqrt(disAvr/Icount)
      write (*,*) '    Number of particles with given displacement'
      imax = MIN(INT(disMax/hDispl)+2,200)
      write (*,*) ' dx             N'
      write(*,'(g11.4,i8)')((i-1)*hDispl,iDispl(i),i=0,imax,10)

      write (17,*) ' Compare two snapshots        =',snap1,'  ',snap2 
      write (17,*) ' Number of particles          =',Icount
      write (17,*) ' Min/Max of coordinates(Mpc/h)=',xmin,xmax
      write (17,*) ' Min/Max of velocities (km/s) =',vmin,vmax
      write (17,*) ' Max displacement (0-cell)    =',disMax
      write (17,*) ' Avr displacement             =',sqrt(disAvr/Icount)
      write (17,*) '    Number of particles with given displacement'
      
      write (17,*) ' dx             N'
      write(17,'(g11.4,i8)')((i-1)*hDispl,iDispl(i),i=0,imax)
      END



C---------------------------------------------------
C                                  Read  current data from disk/tape,
C                                  Open files
C                                  Nrecl is the number of values in a record
C                                  Npage is the number of particles in a page
      SUBROUTINE RDTAPE(snapshot,ifile1,ifile2)
C---------------------------------------------------
      INCLUDE 'PMparameters.h'
      CHARACTER snapshot*(*),file1Name*50,file2Name*50
C                                     Open file on a tape
      file1Name = 'FILES/PMcrda0.'//snapshot//'.DAT'
      file2Name = 'FILES/PMcrs0a0.'//snapshot//'.DAT'
      write(*,'(A)')' File=',file1Name,file2Name
      OPEN(UNIT=ifile1,FILE=file1Name,
     +                FORM='UNFORMATTED', STATUS='UNKNOWN')
C                                 Read control information
C                                 and see whether it has proper
C                                 structure
      READ  (ifile1,err=20,end=20) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,
     +                  NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5
     +                   ,Ocurv,extras
      WRITE (*,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,EKIN1,EKIN2,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
      WRITE (16,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,EKIN1,EKIN2,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
100   FORMAT(1X,'Header=>',A45,/
     +           1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',3E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I8,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3)
      IF(NROW.NE.NROWC) THEN
         WRITE (*,*)
     +            ' NROW in PARAMETER and in TAPE-FILE are different'
      ENDIF
      IF(NGRID.NE.NGRIDC) THEN
         WRITE (*,*)
     +           ' NGRID in PARAMETER and in TAPE-FILE are different:'
         write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC
      ENDIF
C                                         Open work file on disk
 10   NBYTE = NRECL*4
      NACCES= NBYTE / nbyteword

      OPEN(UNIT=20+ifile2,FILE=file2Name,ACCESS='DIRECT',
     +	               STATUS='UNKNOWN',RECL=NACCES)

 
      REWIND ifile1
      RETURN
 20   write (*,*) ' Error reading the header file: Abort',ifile1
      stop
      END
C--------------------------------------------
C                            Read current PAGE of particles  from disk
C                            NRECL - length of ROW block in words
      SUBROUTINE GETROW(IROW,Ifile)
C--------------------------------------------
      INCLUDE 'PMparameters.h'
         READ  (20+Ifile,REC=IROW) RECDAT
      RETURN
      END
C---------------------------------- Read in variables      
      REAL FUNCTION INPUT(text)
C------------------------------------------------
      Character text*(*)
          write (*,'(A,$)')text
          read (*,*) x
          INPUT =x
      Return
      End
