C ______________________________ Select particles, find density C C C Scales internal PM coordinates and velocities C to coordinates in Mpc/h and km/s INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' REAL AMprevious(Np),AMprevious2(Np) REAL INPUT Character FileASCII*50 C................................................................... C Read data and open files WRITE (*,'(A,$)') 'Enter File Name for processed data = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') c OPEN(17,FILE=FileASCII,STATUS='UNKNOWN',FORM='UNFORMATTED') CALL RDTAPE Cell = 0.15 ! Cell size in Mpc/h write (*,*) ' RDTAPE is done' Fraction = 1. c Fraction =INPUT(' Enter fraction of particles =') If(Fraction.le.0.)Stop Ifraction = INT(1./Fraction+0.1) write (*,*) ' Every ',Ifraction,' particle will be selected' write (*,*) + ' Enter Center and Radius XYZ R (Mpc/h)=' read (*,*) Xc,Yc,Zc,Radius c Rsearch =0.001 Rsearch =INPUT(' Enter radius to look for neighbors (Mpc/h)=') WRITE (17,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv,Xc,Yc,Zc,Radius,Rsearch 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, + ' Center=',3f8.3,' Rad=',f7.2,' Rsearch=',f7.3) Box =INPUT(' Enter box size in comoving Mpc/h =') BoxV =Box*100. ! Box size in km/s ScaleV = BoxV/AEXPN/NGRID ! scale factor for Velocities ScaleC = Box/NGRID ! scale factor for Coordinates xshift =0. yshift =0. zshift =0. c CALL GetCentre(ScaleC,ScaleV,Xc,Yc,Zc, c & xshift, yshift , zshift ) c write (*,*) ' Center=',Xc,Yc,Zc c Radius =0.5 CALL ReadPnts(N,Nsm,ScaleC,ScaleV,xmin,ymin,zmin, & xshift,yshift,zshift,Xc,Yc,Zc,Radius,Ifraction) CALL List(Box,N,xmin,ymin,zmin) ScaleInc = 4. c ScaleInc =INPUT(' Enter factor for larger radius or 1 for none=') c WRITE (17) AEXPN,ISTEP,Box,Rsearch,ScaleInc write (*,*) ' Center=',Xc,Yc,Zc If(ScaleInc.lt.1.1)ScaleInc=1. ! background scale Radius =Rsearch *ScaleInc If(ScaleInc.gt.1.1)Then write (*,*) ' Search neighbors inside large radius=',Radius Do i=1,Nsm c if(i.eq.100000)STOP if(i/100000*100000.eq.i)write (*,*) ' particle=',i,' of ',N X =Xp(i) Y =Yp(i) Z =Zp(i) Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Radius) AMprevious(i)=AMnew c Rr =Radius *ScaleInc c Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Rr) AMprevious2(i)=Amnew Enddo Else write (*,*) ' Center=',Xc,Yc,Zc Do i=1,Nsm AMprevious(i)=0. AMprevious2(i)=0. EndDo Endif Radius =Rsearch write (*,*) ' Search neighbors inside small radius=',Radius Do i=1,Nsm if(i/100000*100000.eq.i)write (*,*) ' particle=',i,' of ',N X =Xp(i) Y =Yp(i) Z =Zp(i) Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Radius) If(ScaleInc.lt.1.1)AMprevious(i)=0. Amextr =Amnew -AMprevious(i)/ScaleInc**3 If(AMprevious(i).lt.1.05*Amnew)Amextr =0.001 If(Amextr.lt.0.01)Amextr=0.001 NpMax = INT(MAX(Amnew,AMprevious(i))) c If(iWeight(i).lt.2.) c & write (17) X-Xc,Y-Yc,Z-Zc,Int(Amnew), c & INT(AMprevious(i)),INT(AMprevious2(i)) !,AMextr, If(NpMax.gt.1.and.iWeight(i).lt.2.) & write (17,200) X-Xc,Y-Yc,Z-Zc,Int(Amnew), & INT(AMprevious(i)),AMextr, & INT(iWeight(i)) 200 Format(3F9.4,i6,i7,g11.3,i3) Enddo write (*,*) ' Scaled Coordinates were written to:', & FileASCII write (*,*) ' Number of particles =',N END C------------------------------------------ C Read particles C SUBROUTINE GetCentre(ScaleC,ScaleV,Xc,Yc,Zc, & xshift, yshift , zshift ) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Dimension ListPart(10000) Character*90 Line Logical Tcen Open(18,file='CenterParticles.dat') Tcen =.false. Box =ScaleC*NGRID ifile =1 If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies C -------------------------------------- If(Tcen)Then ! assign the center Xc = 0.988 + xshift Yc = 0.923 + yshift Zc = 11.638 + zshift Vxcc = -98.643 Vycc = 110.17 Vzcc = 35.71 Radius = 0.002 c find particles close to the center Radius2 = Radius**2 ! square of radius write (*,*) ' Center =',Xc,Yc,Zc,' shift=',xshift, & yshift, zshift,' Radius=', Radius write (18,*) ' Center =',Xc,Yc,Zc,' shift=',xshift, & yshift, zshift,' Radius=', Radius Vxcen =0. Vycen =0. Vzcen =0. Ncen =0 C Loop over particles DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last write (*,*)' Read page=',IROW,' file=',ifile,' N=',N iL = NPAGE*(IROW-1) CALL GETROW(IROW,ifile) DO IN=1,In_page X =ScaleC* (XPAR(IN)-1.) Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) X =X+xshift Y =Y+yshift Z =Z+zshift IF(X.lt.0.)X=X +Box IF(Y.lt.0.)Y=Y +Box IF(Z.lt.0.)Z=Z +Box IF(X.GE.Box)X=X -Box IF(Y.GE.Box)Y=Y -Box IF(Z.GE.Box)Z=Z -Box Vxs=ScaleV* VX(IN) Vys=ScaleV* VY(IN) Vzs=ScaleV* VZ(IN) If(Nspecies.eq.0)Then ! find weight of particle W =PARTW Else Ipart =IN+iL ! current particle number Do i =1,Nspecies if(Ipart.le.lspecies(i))Then W =wspecies(i) goto 50 endif EndDo EndIf 50 Ipart =IN+iL ! current particle number iWeigh = INT(W/wspecies(1)+0.1) ! W ! ! W ! particles weight iWeight(Ipart) If(iWeigh.lt.1.001)Then DD =(X-Xc)**2+(Y-(Yc))**2 & +(Z-(Zc))**2 If(DD.lt.Radius2)Then ! inside the limits -- take it Ncen =Ncen +1 Vxcen = Vxcen +Vxs Vycen = Vycen +Vys Vzcen = Vzcen +Vzs Vv = SQRT((Vxs-Vxcc)**2+ & (Vys-Vycc)**2+(Vzs-Vzcc)**2) If(Vv.lt.30.) write (18,'(i8,3f9.4,4f8.2)') Ipart,X,Y,Z, & Vxs-Vxcc,Vys-Vycc,Vzs-Vzcc,Vv EndIf Endif ENDDO ENDDO 300 write (*,*) ' Number of particle at the center=',Ncen write (*,*) ' Radius =',Radius,' Centre=',Xc,Yc,Zc write (*,*) ' Vmean =', & Vxcen/Ncen,Vycen/Ncen,Vzcen/Ncen Close (18) Return c ELSE ! read the list of particles and find the center C -------------------------------------- c read(18,'(A)') Line read(18,'(A)') Line Ncen =0 Xc = 0. Yc = 0. Zc = 0. 400 read(18,*,err=410,end=410) In,X,Y,Z Ncen =Ncen +1 ListPart(Ncen) =In goto 400 c find coords of the particles 410 Ncurrent =1 IpCurrent =ListPart(Ncurrent) C Loop over particles DO IROW=1,Npages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last write (*,*)' Read page=',IROW,' file=',ifile,' Nc=',Ncurrent iL = NPAGE*(IROW-1) CALL GETROW(IROW,ifile) DO IN=1,In_page Ipart =IN+iL ! current particle number If(Ipart.eq.IpCurrent)Then X =ScaleC* (XPAR(IN)-1.) Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) X =X+xshift Y =Y+yshift Z =Z+zshift IF(X.lt.0.)X=X +Box IF(Y.lt.0.)Y=Y +Box IF(Z.lt.0.)Z=Z +Box IF(X.GE.Box)X=X -Box IF(Y.GE.Box)Y=Y -Box IF(Z.GE.Box)Z=Z -Box Xc =Xc +X Yc =Yc +Y Zc =Zc +Z If(Ncurrent.ge.Ncen) Goto 500 Ncurrent =Ncurrent +1 IpCurrent =ListPart(Ncurrent) EndIf EndDo EndDo 500 write (*,*) ' Number of centers found=',Ncurrent write(*,*) ' it should be equal to =',Ncen Xc =Xc/Ncurrent Yc =Yc/Ncurrent Zc =Zc/Ncurrent write (*,*) ' Center=',Xc,Yc,Zc ENDIF Return End C------------------------------------------ C Read particles C SUBROUTINE ReadPnts(N,Nsm,ScaleC,ScaleV,xmin,ymin,zmin, & xshift,yshift ,zshift ,Xc,Yc,Zc,Radius,Ifraction) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Integer Xdist(0:25),Ydist(0:25),Zdist(0:25) Box =ScaleC*NGRID Radius2 =Radius**2 N =0 ! Number of particles Nsm = 0 xmax =-1.e+9 xmin = 1.e+9 ymax =-1.e+9 ymin = 1.e+9 zmax =-1.e+9 zmin = 1.e+9 vmax =-1.e+9 vmin = 1.e+9 Do i=0,25 Xdist(i) =0 Ydist(i) =0 Zdist(i) =0 Enddo Icount =0 Icount_p =Icount-1 ifile =1 If(Nspecies.eq.0)Then ! old case: 1 complete set Npages =NROW ! Number of pages N_in_last=NPAGE ! Number of particles in last page Else ! multiple masses N_particles =lspecies(Nspecies) ! Total number of particles Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf 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,' file=',ifile,' N=',N iL = NPAGE*(IROW-1) if(Icount_p.eq.Icount.and.Icount_p.gt.100)goto300 ! stop if number does not change CALL GETROW(IROW,ifile) Icount_p=Icount DO IN=1,In_page X =ScaleC* (XPAR(IN)-1.) Y =ScaleC* (YPAR(IN)-1.) Z =ScaleC* (ZPAR(IN)-1.) C write (*,'(i8,10g11.4)') IN+iL, c & X,Y,Z,XPAR(IN),YPAR(IN),ZPAR(IN) X =X+xshift Y =Y+yshift Z =Z+zshift IF(X.lt.0.)X=X +Box IF(Y.lt.0.)Y=Y +Box IF(Z.lt.0.)Z=Z +Box IF(X.GE.Box)X=X -Box IF(Y.GE.Box)Y=Y -Box IF(Z.GE.Box)Z=Z -Box DD =(X-Xc)**2+(Y-(Yc))**2 & +(Z-(Zc))**2 If(DD.lt.Radius2)Then ! inside the limits -- take it Vxs=ScaleV* VX(IN) Vys=ScaleV* VY(IN) Vzs=ScaleV* VZ(IN) Icount =Icount +1 vmax =MAX(vmax,Vxs,Vys,Vzs) vmin =MIN(vmin,Vxs,Vys,Vzs) If(Nspecies.eq.0)Then ! find weight of particle W =PARTW Else Ipart =IN+iL ! current particle number Do i =1,Nspecies if(Ipart.le.lspecies(i))Then W =wspecies(i) goto 50 endif EndDo EndIf 50 If(Icount/Ifraction*Ifraction.eq.Icount)Then N = N +1 If(Ipart.le.lspecies(1))Nsm=Nsm +1 IF(N.gt.Np)Then ! too many particles write (*,*) ' Too many particles. Max=',Np STOP Endif Xp(N) =X Yp(N) =Y Zp(N) =Z Ipart =IN+iL ! current particle number iWeight(N) = INT(W/wspecies(1)+0.1) ! W ! ! W ! particles weight iWeight(Ipart) If(iWeight(N).lt.1.001)Then xmax =MAX(xmax,Xp(N)) xmin =MIN(xmin,Xp(N)) ymax =MAX(ymax,Yp(N)) ymin =MIN(ymin,Yp(N)) zmax =MAX(zmax,Zp(N)) zmin =MIN(zmin,Zp(N)) ii =Min(INT(X/Box*25.),25) jj =Min(INT(Y/Box*25.),25) kk =Min(INT(Z/Box*25.),25) Xdist(ii) =Xdist(ii) +1 Ydist(jj) =Ydist(jj) +1 Zdist(kk) =Zdist(kk) +1 EndIf Endif Endif ! end test for center ENDDO ENDDO 300 write (*,*) ' Number of particles selected:=',N,' Small=',Nsm write (*,*) ' Limits =',xmin,xmax write (*,*) ' =',ymin,ymax write (*,*) ' =',zmin,zmax write (*,*) ' Radius =',Radius,' Centre=',Xc,Yc,Zc write (*,'(i3,2x,i9,2x,i9,2x,i9)') & (i,Xdist(i),Ydist(i),Zdist(i),i=0,25) Return End C-------------------------------------------------------------- C Make linker lists of particles in each cell SUBROUTINE List(Box,N,xmin,ymin,zmin) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Do i=1,N Lst(i)=-1 EndDo Do k=Nm,Nb Do j=Nm,Nb Do i=Nm,Nb Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,N i=INT((Xp(jp)-xmin)/Cell) j=INT((Yp(jp)-ymin)/Cell) k=INT((Zp(jp)-zmin)/Cell) i=MIN(MAX(Nm,i),Nb) j=MIN(MAX(Nm,j),Nb) k=MIN(MAX(Nm,k),Nb) Lst(jp) =Label(i,j,k) Label(i,j,k) =jp EndDo Return End C-------------------------------------------------------------- C Find all neibhours for a center C Xc,Yc,Zc - center; ww-its weight C SUBROUTINE Neib(Amnew,xmin,ymin,zmin,Xc,Yc,Zc,Radius) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' C Initiate counters Amnew =0. Radius2=Radius**2 c limits for Label i1=INT((Xc-xmin-Radius)/Cell) j1=INT((Yc-ymin-Radius)/Cell) k1=INT((Zc-zmin-Radius)/Cell) i1=MIN(MAX(Nm,i1),Nb) j1=MIN(MAX(Nm,j1),Nb) k1=MIN(MAX(Nm,k1),Nb) i2=INT((Xc-xmin+Radius)/Cell) j2=INT((Yc-ymin+Radius)/Cell) k2=INT((Zc-zmin+Radius)/Cell) i2=MIN(MAX(Nm,i2),Nb) j2=MIN(MAX(Nm,j2),Nb) k2=MIN(MAX(Nm,k2),Nb) C Look for neibhours c nn =0 Do k=k1,k2 Do j=j1,j2 Do i=i1,i2 jp =Label(i,j,k) 10 If(jp.ne.0)Then dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2 If(dd.lt.Radius2)Then c nn =nn +1 Amnew=Amnew + iWeight(jp) EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo c write (*,100) nn,i1,i2,j1,j2,k1,k2,Xc,Yc,Zc c 100 format(' Nn=',i6,3(2i4,2x),3f8.3) Return End