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