C ______________________________ INCLUDE 'PMbox.h' REAL INPUT Character*90 Line PI =3.1415926535 C................................................................... C Read data and open files Box =160. BoxV =Box*100. ! Box size in km/s write(*,*) ' Enter: Center(xyz) Radius N_part_limit V_circ_limit ' read (*,*) Xc,Yc,Zc,Rsphere,N_plim,V_clim c Open(1,file='Catshort.999.small.DAT') Open(1,file='Catshort.DAT') open(18,file='Box.dat') CALL ReadHalos(N_plim,V_clim) c CALL CONE(ux,uy,uz,ACone,ScaleC) CALL SPHERE(Rsphere,Xc,Yc,Zc) c CALL DumpParticles(z0,dz0,ScaleC) STOP END C--------------------------------------- C Read Particles SUBROUTINE ReadHalos(N_plim,V_clim) C--------------------------------------- INCLUDE 'PMbox.h' Character*90 Line Lines_header =20 Do i=1,Lines_header Read(1,'(A)') Line If(i.le.4)write(18,'(A)') Line write(*,'(A)') Line ENDDO Nh =0 10 read(1,*,err=20,end=20) a,x,y,z,vvx,vvy,vvz,am,ar,v1,v2,n1,n2 If(n2.ge.N_plim.and.v2.ge.V_clim) Then ! take the halo Nh =Nh + 1 Xm(Nh) =x Ym(Nh) =y Zm(Nh) =z Vx(Nh) =vvx Vy(Nh) =vvy Vz(Nh) =vvz Amc(Nh) =am Rmc(Nh) =ar Vrm(Nh) =v2 Endif goto 10 20 write (*,*)'Number of halos read =',Nh,' Np>',N_plim,' Vc>',V_clim RETURN END C--------------------------------------- C Particles in a sphere SUBROUTINE Sphere(Rsphere,Xc,Yc,Zc) C--------------------------------------- C Find points close to the center INCLUDE 'PMbox.h' write (*,*) ' Center=',Xc,Yc,Zc,' R=',Rsphere V =0. V2 =0. Vm =0. Vm2=0. Nn = 0 R2 = Rsphere**2 Ic =0 Vvx =0. Vvy =0. Vvz =0. Vcx =-215.1 Vcy =72.3 Vcz =-5.6 DO IN=1,Nh ! find central halo dX=(Xm(IN)-Xc) dY=(Ym(IN)-Yc) dZ=(Zm(IN)-Zc) dd =dX**2 +dY**2 +dZ**2 c write (*,50)IN,Xm(IN),Ym(IN),Zm(IN),dX,dY,dZ c 50 format(i5,3f8.2,3g11.3) If(dd.lt.R2)Then Nn =Nn +1 Vvx =Vvx +Vx(IN) Vvy =Vvy +Vy(IN) Vvz =Vvz +Vz(IN) c If(dd.lt.1.e-4)Then c Ic =IN c Vcx =Vx(Ic) c Vcy =Vy(Ic) c Vcz =Vz(Ic) c Endif EndIf ENDDO If(Nn.le.1)Then write(*,*) ' Too few halos inside the radius: ',Nn stop Endif Vvx =Vvx/Nn Vvy =Vvy/Nn Vvz =Vvz/Nn If(Ic.eq.0)Then ! no center halo was found Vcx =Vvx Vcy =Vvy Vcz =Vvz endif Nc = 0 DO IN=1,Nh ! Main statistics dX=(Xm(IN)-Xc) dY=(Ym(IN)-Yc) dZ=(Zm(IN)-Zc) dd =dX**2 +dY**2 +dZ**2 If(dd.lt.R2.and.dd.gt.1.e-4)Then d = max(sqrt(dd+1.e-10),1.e-10) Nc =Nc +1 wx =Vx(IN)-Vcx wy =Vy(IN)-Vcy wz =Vz(IN)-Vcz vr =(dX*wx+dY*wy+dZ*wz)/d ww =sqrt(wx**2+wy**2+wz**2) c Xn =dX*(1.+vr/100./d)+Xc c Yn =dY*(1.+vr/100./d)+Yc c Zn =dZ*(1.+vr/100./d)+Zc write (*,90)d,ww,vr,100.*d+vr,Amc(IN),Vrm(IN) 90 format(f8.3,f9.2,2f8.2, & g11.3,f7.2) V2 = V2 + ww**2 Vm2= Vm2 + vr**2 endif Enddo write (*,*) ' There are ',Nc,' particles inside R=',Rsphere write (*,*) & ' Center(Mpc/h) =',Xc,Yc,Zc write (*,*) & ' V(km/s)=', Vcx, Vcy, Vcz write (*,*) ' Rms Radial velocity=',sqrt(Vm2/Nc), & ' Rms total velocity=',sqrt(V2/Nc) Return End C--------------------------------------- C write particles in a slice SUBROUTINE DumpParticles(z0,dz0,ScaleC) C--------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMfields.h' C Find points close to the center Xc =NGRID/2.+1. Yc =NGRID/2.+1. Zc =NGRID/2.+1. Rc =1. 10 R2 =Rc**2 C Loop over particles Nc =0 Vxc =0. Vyc =0. Vzc =0. DO IN=1,Np dX=(Xp(IN)-Xc) !coords dY=(Yp(IN)-Yc) dZ=(Zp(IN)-Zc) dd =dX**2 +dY**2 +dZ**2 If(dd.lt.R2)Then Nc =Nc+1 Vxc =Vxc +Vxp(IN) Vyc =Vyc +Vyp(IN) Vzc =Vzc +Vzp(IN) endif Enddo If(Nc.eq.0)Then Rcold =Rc Rc =Rc*1.2 write (*,*) ' Too small radius=',Rcold,' Increase it to ',Rc goto 10 endif Vxc =Vxc/Nc Vyc =Vyc/Nc Vzc =Vzc/Nc Xc =ScaleC*(Xc-1.) Yc =ScaleC*(Yc-1.) Zc =ScaleC*(Zc-1.) write (*,*) ' There are ',Nc,' particles inside Rc=',Rc*ScaleC, & ' Center(Mpc/h) =',Xc, & Yc,Zc, & ' V(km/s)=', Vxc, Vyc, Vzc C Loop over particles DO IN=1,Np X=ScaleC*(Xp(IN)-1.) !coords Y=ScaleC*(Yp(IN)-1.) Z=ScaleC*(Zp(IN)-1.) If(Z.gt.z0.and.Z.lt.z0+dz0)Then dX=(X-Xc) !coords dY=(Y-Yc) dZ=(Z-Zc) dd =dX**2 +dY**2 +dZ**2 d =max(sqrt(dd+1.e-10),1.e-10) ux =Vxp(IN)-Vxc uy =Vyp(IN)-Vyc uz =Vzp(IN)-Vxc vr =(dX*ux+dY*uy+dZ*uz)/d Xn =dX*(1.+vr/100./d)+Xc Yn =dY*(1.+vr/100./d)+Yc Zn =dZ*(1.+vr/100./d)+Zc write (18,100) X,Y,Z,ux,uy,uz,Wp(IN),vr, & Xn,Yn,Zn 100 format(3f8.3,3f8.1,g11.3,4f8.3) endif Enddo Return End C--------------------------------------- C Rotate around Y-axis by angle alpha SUBROUTINE RotateYParticles(alpha) C--------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMfields.h' XN =NGRID+1 c XC =(NGRID+1.)/2. ! center of rotation c ZC =(NGRID+1.)/2. cosA =cos(alpha*3.1415926/180.) sinA =sin(alpha*3.1415926/180.) C Loop over particles DO IN=1,Np X=Xp(IN) -XC !coords Y=Yp(IN) Z=Zp(IN) -ZC Xnew =X*cosA -Z*sinA + XC Znew =Z*cosA +X*sinA + ZC Vxp(IN) =Vxp(IN)*cosA-Vzp(IN)*sinA Vyp(IN) =Vyp(IN) Vzp(IN) =Vzp(IN)*cosA+Vxp(IN)*sinA X =Xnew If(X.ge.XN)X=X-NGRID If(X.lt.1.)X=X+NGRID Z =Znew If(Z.ge.XN)Z=Z-NGRID If(Z.lt.1.)Z=Z+NGRID Xp(IN) =X Yp(IN) =Y Zp(IN) =Z Enddo Return End C--------------------------------------- C Shift particles with periodical conditions SUBROUTINE ShiftParticles(xsh,ysh,zsh) C--------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMfields.h' XN =NGRID+1 C Loop over particles DO IN=1,Np X=Xp(IN) !coords Y=Yp(IN) Z=Zp(IN) X =X+xsh If(X.ge.XN)X=X-NGRID If(X.lt.1.)X=X+NGRID Y =Y+ysh If(Y.ge.XN)Y=Y-NGRID If(Y.lt.1.)Y=Y+NGRID Z =Z+zsh If(Z.ge.XN)Z=Z-NGRID If(Z.lt.1.)Z=Z+NGRID Xp(IN) =X Yp(IN) =Y Zp(IN) =Z c Vxpb(IN) =Vxp(IN) c Vypb(IN) =Vyp(IN) c Vzpb(IN) =Vzp(IN) Enddo Return End C---------------------------------- Read in variables REAL FUNCTION INPUT(text) C------------------------------------------------ Character text*(*) write (*,'(A,$)')text read (*,*) x INPUT =x Return End