C ______________________________ Convert PM format to ASCII C November 1997 A. Klypin (aklypin@nmsu.edu) C C Scales internal PM coordinates and velocities C to coordinates in Mpc/h and km/s C In PM the coordinates are in the range 1 - (NGRID+1) C velocities are P = a_expansion*V_pec/(x_0H_0) C where x_0 = comoving cell_size=Box/Ngrid C H_0 = Hubble at z=0 C C NROW = number of particles in 1D C NGRID= number of cells in 1D INCLUDE 'PMparameters.h' INCLUDE 'PMrotation.h' REAL INPUT Character FileASCII*50 C................................................................... C Read data and open files WRITE (*,'(A,$)') 'Enter Name for output ASCII file = ' READ (*,'(A)') FileASCII OPEN(17,FILE=FileASCII,STATUS='UNKNOWN') WRITE (*,'(A,$)') 'Enter Name for input Centers = ' READ (*,'(A)') FileASCII OPEN(30,FILE=FileASCII,STATUS='UNKNOWN') CALL RDTAPE 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) Box =25. c 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 Nlines =25 CALL ReadMax(Ncent,Nlines) Radius =0. Do i=1,Ncent Radius =MAX(Rmc(i),Radius) EndDo CALL SetWeights CALL ReadPnt(N,ScaleC,ScaleV) ! particles CALL Points(Box,N,Ncp,Radius) ! Make more points periodically CALL Rotation(Ncent,Ncp) Stop END C-------------------------------------------------------------- C SUBROUTINE Rotation(Ncent,Ncp) C Lambda = J*sqrt(/2)/[V_circ_vir*R_vir] C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMrotation.h' real*8 wm,V2,ax,ay,az,Vvx,Vvy,Vvz, & ddx,ddy,ddz Do i=1,Ncent ! loop over centers nn =0 wm=0. V2 =0. ax =0. ay =0. az =0. x = Xm(i) y = Ym(i) z = Zm(i) wx =Vmx(i) wy =Vmy(i) wz =Vmz(i) R =Rmc(i) R2 =R**2 Vvx =0. Vvy =0. Vvz =0. ddx =0. ddy =0. ddz =0. Do j=1,Ncp ! loop over particles dx = Xp(j)-x dy = Yp(j)-y dz = Zp(j)-z dd =dx**2+dy**2+dz**2 If(dd.lt.R2)Then ! inside halo nn = nn+1 ww=iWeight(i) ux =VXp(j) uy =VYp(j) uz =VZp(j) Vvx =ux*ww +Vvx Vvy =uy*ww +Vvy Vvz =uz*ww +Vvz ddx =ddx + dx*ww ddy =ddy + dy*ww ddz =ddz + dz*ww wm=wm +ww V2 =V2 + (ux**2+uy**2+uz**2)*ww ax = ax + (dy*uz-dz*uy) ay = ay + (dz*ux-dx*uz) az = az + (dx*uy-dy*ux) Endif EndDo If(nn.eq.0)Then write(*,*) ' Error: no particles' stop Else ax = (ax - (ddy*Vvz-ddz*Vvy)/wm)/wm*AEXPN ay = (ay - (ddz*Vvx-ddx*Vvz)/wm)/wm*AEXPN az = (az - (ddx*Vvy-ddy*Vvx)/wm)/wm*AEXPN aJm = sqrt(ax**2+ay**2+az**2) V2 = (V2 - (Vvx**2+Vvy**2+Vvz**2)/wm)/wm V =sqrt(max(V2,1.e-10)) Vcirc =6.581e-5*sqrt(Amc(i)/Rmc(i)/AEXPN) aLambda1(i) =aJm*V/sqrt(2.)/Vcirc**2/Rmc(i)/AEXPN aLambda2(i) =aJm/Vcirc/Rmc(i)/AEXPN aJ(i) =aJm Vvx =Vvx/wm Vvy =Vvy/wm Vvz =Vvz/wm write (*,100) i,nn,wm,Vcirc,V,m1(i),m2(i),m3(i), & aJ(i),aLambda1(i) ,aLambda2(i), & Vvx,Vvy,Vvz,Vmx(i),Vmy(i),Vmz(i) write (17,120) AEXPN,x,y,z, Vvx,Vvy,Vvz,aMc(i), & Rmc(i)*1000.,Vrm(i),Vcmax(i), & m1(i),m2(i),m3(i),m4(i), & aLambda1(i),aJ(i),ax,ay,az EndIf EndDo 100 format(i3,i8,g11.3,2f8.2,2i8,i5,3g11.3,/ & 3x,3F8.2,3x,3F8.2,g11.3) 120 format(F7.4,F8.3,2F8.3,3F7.1,g10.3,f8.2,2F7.2,2i8,i7,i5, & f7.4,4g12.4 ) Return End C-------------------------------------------------------------- C Add points around the Box to take into account periodical conditions C It replicates points periodically and keeps those which are at distance C less than Radius from the Box. C N = number of points inside the Box C Ncp = total number of points SUBROUTINE Points(Box,N,Ncp,Radius) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMrotation.h' Logical Inside B1 =-Radius B2 =Box+Radius B3 =Box-Radius Ncp =N Do i=1,N ! Add dark matter particles x0 =xp(i) y0 =yp(i) z0 =zp(i) Inside = (x0.gt.Radius).and.(x0.lt.B3) Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3)) Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3)) If(.not.Inside)THEN ux =VXp(i) uy =VYp(i) uz =VZp(i) ww=iWeight(i) Do k1 = -1,1 zr =z0+k1*Box If(zr.GT.B1.AND.zr.LT.B2)Then kk =k1**2 Do j1 = -1,1 yr =y0+j1*Box If(yr.GT.B1.AND.yr.LT.B2)Then jj = kk +j1**2 Do i1 = -1,1 xr =x0+i1*Box ii =jj +i1**2 If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0) . Then Ncp =Ncp +1 If(Ncp.le.Np)Then xp(Ncp) =xr yp(Ncp) =yr zp(Ncp) =zr VXp(Ncp) =ux VYp(Ncp) =uy VZp(Ncp) =uz iWeight(Ncp) =ww EndIf EndIf EndDo EndIf EndDo EndIf EndDo EndIf EndDo If(Ncp.Gt.Np)Then write(*,*)' Too many points:',Ncp,' Maximum is =', Np STOP EndIf Return End C-------------------------------------------------- C Set Weights of particles for C fast access SUBROUTINE SetWeights C-------------------------------------------------- INCLUDE 'PMparameters.h' Wpr =NGRID**3/FLOAT(NROW**3) If(Nspecies.eq.0)Then ! old constant weights N_particles =lspecies(1) Ww = wspecies(1)/Wpr Do i=1,N_particles iWeight(i) =Ww EndDo Else N_particles =lspecies(Nspecies) jstart =1 write (*,*) ' N_Particles=',N_particles If(N_particles.le.0)Then write (*,*) ' Wrong number of particles:',N_particles write (*,*) ' Nspecies=',Nspecies,' N=',lspecies STop endif Do j=1,Nspecies jend =lspecies(j) Do k=jstart ,jend iWeight(k) =wspecies(j)/Wpr EndDo jstart =jend EndDo EndIf write (*,*) ' Set Weights for ',N_particles,' particles:', & ' w_first=' ,iWeight(1),' w_last=', iWeight(N_particles) RETURN END C--------------------------------------------------------- C Read particles SUBROUTINE ReadPnt(N,ScaleC,ScaleV) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMrotation.h' N =0 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 c write (*,*) ' N_in_last=',N_in_last DO IROW=1, Npages ! Loop over particle pages In_page =NPAGE If(IROW.eq.Npages)In_page =N_in_last iL = NPAGE*(IROW-1) C Loop over particles READ (21,REC=IROW) RECDAT DO IN=1,In_page N = N +1 Xp(N) =ScaleC*(XPAR(IN)-1.) Yp(N) =ScaleC*(YPAR(IN)-1.) Zp(N) =ScaleC*(ZPAR(IN)-1.) VXp(N)=ScaleV*VX(IN) VYp(N)=ScaleV*VY(IN) VZp(N)=ScaleV*VZ(IN) Ipart =IN+iL ! current particle number iWeight(N) =iWeight(Ipart) ! particles weight ENDDO ENDDO Return End C--------------------------------------------------------- C Read list of maxima c Xnew =X-1: coord. now are 0-NGRID SUBROUTINE ReadMax(Ncent,Nlines) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMrotation.h' Character*100 Line Do i=1,Nlines read(30,'(A)') Line write(*,'(A)') Line EndDo Ncent =0 xxm =1000. xxn =-1000. 10 read (30,*,err=5,end=5) Aa,X0, Y0, Z0, VvX,VvY,VvZ,aM, & Rv,V1,V2,nn1,nn2,nn3,nn4 Ncent=Ncent+1 If(Ncent.le.Nc)Then Xm(Ncent) =X0 Ym(Ncent) =Y0 Zm(Ncent) =Z0 Vmx(Ncent)=VvX Vmy(Ncent)=VvY Vmz(Ncent)=VvZ Amc(Ncent)=aM Rmc(Ncent)=Rv / 1000. ! scale to Mpc/h Vrm(Ncent)=V1 Vcmax(Ncent)=V2 m1(Ncent)=nn1 m2(Ncent)=nn2 m3(Ncent)=nn3 m4(Ncent)=nn4 xxm =min(xxm,Xm(Ncent)) xxn =max(xxn,Xm(Ncent)) EndIf GoTo 10 5 write(*,*) ' Read Centers: X{Min,Max}=',xxm,xxn,' Nc=',Ncent If(Ncent.GT.Nc)Then write (*,*) ' Too many centers: Max=',Nc Stop EndIf Return End