C ______________________________ Select particles, find density C Simulations of a galaxy 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),AMlast(Np),Daverg(0:500) Real*8 Xc,Yc,Zc DIMENSION Iaverg(0:500) REAL INPUT Character FileASCII*50 EQUIVALENCE (Ndisk,extras(90)) 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 = 7.e-4 ! Cell size in Mpc/h dRaver = 0.00035 ! step in radius for the averge number of neighbors 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*1000.,Rsearch*1000. 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 (km/s) ScaleC = Box/NGRID*AEXPN/hubble ! scale factor for Coordinates (real Mpc) 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) Call ChangePnts(N,xmin,ymin,zmin,Xc,Yc,Zc,Box) CALL List(Box,N,xmin,ymin,zmin) c ScaleInc = 4. 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 c --------------------------- write (*,*) ' Search neighbors inside large radius=',Radius Do i=1,Nsm c if(i.eq.100000)STOP if(i/1000*1000.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 Enddo Else write (*,*) ' Center=',Xc,Yc,Zc Do i=1,Nsm AMprevious(i)=0. EndDo Endif Radius =Rsearch c --------------------------- write (*,*) ' Search neighbors inside small radius=',Radius,Box Do i=0,500 Daverg(i) =0. Iaverg(i) =0. EndDo Do i=1,Nsm if(i/1000*1000.eq.i)write (*,*) ' particle=',i,' of ',N X =Xp(i) Y =Yp(i) Z =Zp(i) dist = sqrt((X-Box/2.)**2+(Y-Box/2.)**2+(Z-Box/2.)**2) c Radius =Rsearch*sqrt(dist*1000./0.3) Call Neib(Amnew,xmin,ymin,zmin,X,Y,Z,Radius) AMlast(i)=AMnew dist = sqrt((X-Box/2.)**2+(Y-Box/2.)**2+(Z-Box/2.)**2) irad = MIN(INT(dist/dRaver),500) Daverg(irad) = Daverg(irad) +AMnew Iaverg(irad) = Iaverg(irad) + 1 EndDo c --------------------------- Do i=0,500 Daverg(i) =Daverg(i)/max(Iaverg(i),1) if(Iaverg(i).gt.10)write (*,*) i,i*dRaver,Daverg(i),Iaverg(i) EndDo Daverg(1) =MAX(Daverg(2),Daverg(2)) Do i=1,Nsm X =Xp(i) Y =Yp(i) Z =Zp(i) dist = sqrt((X-Box/2.)**2+(Y-Box/2.)**2+(Z-Box/2.)**2) irad = MIN(INT(dist/dRaver),500) irad2 = MIN(irad+1,500) Dav = (Daverg(irad)*((irad+1)-dist/Draver)+ & Daverg(irad+1)*(dist/dRaver-irad)) Amextr =Amlast(i) -AMprevious(i)/ScaleInc**3 If(AMprevious(i).lt.1.05*Amlast(i))Amextr =0.001 If(Amextr.lt.0.01)Amextr=0.001 NpMax = INT(MAX(Amlast(i),AMprevious(i))) c If(iWeight(i).lt.2.) c & write (17) X-Xc,Y-Box/2.,Z-Box/2.,Int(Amnew), c & INT(AMprevious(i)),INT(AMprevious2(i)) !,AMextr, If(NpMax.gt.1.and.iWeight(i).lt.200.)write (17,200) & (X-Box/2.)*1000.,(Y-Box/2.)*1000.,(Z-Box/2.)*1000., & Int(Amlast(i)), & INT(AMprevious(i)),AMextr,Dav,Daverg(irad), & INT(iWeight(i)) 200 Format(3F9.3,i6,i7,3g11.3,i3) Enddo write (*,*) ' Scaled Coordinates were written to:', & FileASCII write (*,*) ' Number of particles =',N END C------------------------------------------ C Change positions: center, inclination SUBROUTINE ChangePnts(N,xmin,ymin,zmin,Xc,Yc,Zc,Box) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' Real*8 xshift,yshift,zshift,Lx,Ly,Lz,aL,Xc,Yc,Zc,V2 Real*8 Xcc,Ycc,Zcc,Vcx,Vcy,Vcz,dR,dR2,aCos,Ld Ncenter = Min(15000,N) ! find center using central particles only If(Ncenter.eq.0)write (*,*) ' Error!! No particles in ChangePnts' Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. Xc = 0. Yc = 0. Zc = 0. Lx = 0. Ly = 0. Lz = 0. dR = 0. dR2 = 0. Vcx = 0. Vcy = 0. Vcz = 0. V2 =0. xmax =-1.d+9 xmin = 1.d+9 ymax =-1.d+9 ymin = 1.d+9 zmax =-1.d+9 zmin = 1.d+9 c Find center of mass and mean V Do jp=1,Ncenter Xc = Xc +Xp(jp) Yc = Yc +Yp(jp) Zc = Zc +Zp(jp) Vcx = Vcx + Vxp(jp) Vcy = Vcy + Vyp(jp) Vcz = Vcz + Vzp(jp) EndDo Xc = Xc/Ncenter Yc = Yc/Ncenter Zc = Zc/Ncenter Vcx = Vcx/Ncenter Vcy = Vcy/Ncenter Vcz = Vcz/Ncenter xshift = Xcc - Xc ! displacements to bring yshift = Ycc - Yc ! the center to {Xcc,Ycc,Zcc} zshift = Zcc - Zc write (*,100) Ncenter,Xc,Yc,Zc, & 1000.*xshift,1000.*yshift,1000.*zshift, & Vcx,Vcy,Vcz 100 format('--- Parameters of the center: n=',i8,' X(Mpc)=',3g11.4,/ & 30x,'dX(kpc)=',3g11.4,/30x,' V(km/s)=',3g11.4) Do jp=1,N Xp(jp) =Xp(jp) +xshift Yp(jp) =Yp(jp) +yshift Zp(jp) = Zp(jp) +zshift IF(Xp(jp).lt.0.)Xp(jp)=Xp(jp) +Box IF(Yp(jp).lt.0.)Yp(jp)=Yp(jp) +Box IF(Zp(jp).lt.0.)Zp(jp)=Zp(jp) +Box IF(Xp(jp).GE.Box)Xp(jp)=Xp(jp) -Box IF(Yp(jp).GE.Box)Yp(jp)=Yp(jp) -Box IF(Zp(jp).GE.Box)Zp(jp)=Zp(jp) -Box xmax =MAX(xmax,Xp(jp)) xmin =MIN(xmin,Xp(jp)) ymax =MAX(ymax,Yp(jp)) ymin =MIN(ymin,Yp(jp)) zmax =MAX(zmax,Zp(jp)) zmin =MIN(zmin,Zp(jp)) Lx =Lx +(Yp(jp)-Ycc)*(VZp(jp)-Vcz) & -(Zp(jp)-Zcc)*(VYp(jp)-Vcy) Ly =Ly +(Zp(jp)-Zcc)*(VXp(jp)-Vcx) & -(Xp(jp)-Xcc)*(VZp(jp)-Vcz) Lz =Lz +(Xp(jp)-Xcc)*(VYp(jp)-Vcy) & -(Yp(jp)-Ycc)*(VXp(jp)-Vcx) dd = (Xp(jp)-Xcc)**2 + & (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd V2 = V2 + (VXp(jp)-Vcx)**2 +(VYp(jp)-Vcy)**2+(VZp(jp)-Vcz)**2 EndDo Lx =Lx/N Ly =Ly/N Lz =Lz/N V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/N write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/N), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR 110 format(' --- All particles: Xmin/max=',2g11.4,/ & 30x,2g11.4,/30x,2g11.4/, & 25x,' Aver Distance(kpc)=',g11.4, & ' rms distance(kpc)=',g11.4,/ & 25x,' rms Velocity=',g11.4,' km/s',/ & 25x,' Ang.Momentum: direction=',3g11.3,/ & 25x,' Ang.Momentum: L/ =',g11.4,' km/s') c Rotate the system of coordinates If(Lz.le.0.5*aL)Then write (*,*) ' z-component of angular momentum', & ' is too small. Something is very wrong.' Stop EndIf aCos = 0.99999d-0 If(Lz.gt.aCos*aL)Then write (*,*) ' Tilt is too small (Cos>',aCos,'). No change' Return EndIf c Find maxtrix for rotation Lx =Lx/aL Ly =Ly/aL Lz =Lz/aL c Ld =sqrt(Lx**2 +Ly**2) c If(Ly.gt.0.)Then c ex1 = Ly/Ld c ey1 = -Lx/Ld c ex2 = Lx*Lz/Ld c ey2 = Ly*Lz/Ld c ez2 = -Ld c Else c ex1 = -Ly/Ld c ey1 = Lx/Ld c ex2 = -Lx*Lz/Ld c ey2 = -Ly*Lz/Ld c ez2 = Ld c EndIf c ez1 =0. Ld =sqrt(Lx**2 +Lz**2) ex1 = Lz/Ld ey1 = 0. ez1 = -Lx/Ld ex2 = -Lx*abs(Ly)/Ld ey2 = Ld ez2 = -Lz*abs(Ly)/Ld c If(Ly.lt.0.)ey2 = -Ld ex3 =Lx ey3 =Ly ez3 =Lz c Check if the vectors are (1) unit (2) orthogonal prod1 = ex1*Lx +ey1*Ly +ez1*Lz prod2 = ex2*Lx +ey2*Ly +ez2*Lz prod3 = ex1*ex2 +ey1*ey2 +ez1*ez2 a1 = Lx**2 +Ly**2 +Lz**2 a2 = ex1**2 +ey1**2 +ez1**2 a3 = ex2**2 +ey2**2 +ez2**2 write (*,*) ' Vector 1=',ex1,ey1,ez1 write (*,*) ' Vector 2=',ex2,ey2,ez2 write (*,*) ' Vector 3=',ex3,ey3,ez3 write (*,*) ' Vector lengths=', a1,a2,a3,' must be 1' write (*,*) ' Vector products=', prod1,prod2,prod3,' must be 0' if(abs(a1-1.).gt.1.e-2)stop if(abs(a2-1.).gt.1.e-2)stop if(abs(a2-1.).gt.1.e-2)stop if(abs(prod1).gt.1.e-2)stop if(abs(prod2).gt.1.e-2)stop if(abs(prod3).gt.1.e-2)stop xmax =-1.d+9 xmin = 1.d+9 ymax =-1.d+9 ymin = 1.d+9 zmax =-1.d+9 zmin = 1.d+9 dR = 0. dR2 = 0. V2 = 0. Lx = 0. Ly = 0. Lz = 0. Do jp=1,N Xnew = ex1*(Xp(jp)-Xcc) +ey1*(Yp(jp)-Ycc) +ez1*(Zp(jp)-Zcc)+Xcc Ynew = ex2*(Xp(jp)-Xcc) +ey2*(Yp(jp)-Ycc) +ez2*(Zp(jp)-Zcc)+Ycc Znew = ex3*(Xp(jp)-Xcc) +ey3*(Yp(jp)-Ycc) +ez3*(Zp(jp)-Zcc)+Zcc VXnew = ex1*(Vxp(jp)-Vcx) +ey1*(Vyp(jp)-Vcy) +ez1*(Vzp(jp)-Vcz) VYnew = ex2*(Vxp(jp)-Vcx) +ey2*(Vyp(jp)-Vcy) +ez2*(Vzp(jp)-Vcz) VZnew = ex3*(Vxp(jp)-Vcx) +ey3*(Vyp(jp)-Vcy) +ez3*(Vzp(jp)-Vcz) Xp(jp) = Xnew Yp(jp) = Ynew Zp(jp) = Znew IF(Xp(jp).lt.0.)Xp(jp)=Xp(jp) +Box IF(Yp(jp).lt.0.)Yp(jp)=Yp(jp) +Box IF(Zp(jp).lt.0.)Zp(jp)=Zp(jp) +Box IF(Xp(jp).GE.Box)Xp(jp)=Xp(jp) -Box IF(Yp(jp).GE.Box)Yp(jp)=Yp(jp) -Box IF(Zp(jp).GE.Box)Zp(jp)=Zp(jp) -Box Vxp(jp) = VXnew +Vcx Vyp(jp) = VYnew +Vcy Vzp(jp) = VZnew +Vcz xmax =MAX(xmax,Xp(jp)) xmin =MIN(xmin,Xp(jp)) ymax =MAX(ymax,Yp(jp)) ymin =MIN(ymin,Yp(jp)) zmax =MAX(zmax,Zp(jp)) zmin =MIN(zmin,Zp(jp)) Lx =Lx +(Yp(jp)-Ycc)*(VZp(jp)-Vcz) & -(Zp(jp)-Zcc)*(VYp(jp)-Vcy) Ly =Ly +(Zp(jp)-Zcc)*(VXp(jp)-Vcx) & -(Xp(jp)-Xcc)*(VZp(jp)-Vcz) Lz =Lz +(Xp(jp)-Xcc)*(VYp(jp)-Vcy) & -(Yp(jp)-Ycc)*(VXp(jp)-Vcx) dd = (Xp(jp)-Xcc)**2 + & (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd V2 = V2 + (VXp(jp)-Vcx)**2 +(VYp(jp)-Vcy)**2+(VZp(jp)-Vcz)**2 EndDo Lx =Lx/N Ly =Ly/N Lz =Lz/N V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/N write (*,*) '--- After rotation of the coordinates' write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/N),V2,Lx/aL,Ly/aL,Lz/aL,aL/dR Return End C------------------------------------------ C Read particles SUBROUTINE ReadPnts(N,Nsm,ScaleC,ScaleV,xmin,ymin,zmin) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMlists.h' EQUIVALENCE (Ndisk,extras(90)) Integer Xdist(0:25),Ydist(0:25),Zdist(0:25) Box =ScaleC*NGRID Radius2 =Radius**2 write (*,*) ' Number of disk particles=',Ndisk 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 N_particles =Ndisk 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) 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) 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 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 Vxp(N) =Vxs Vyp(N) =Vys Vzp(N) =Vzs 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 ENDDO ENDDO write (*,*) ' Number of particles selected:=',N,' Small=',Nsm write (*,*) ' Limits =',xmin,xmax write (*,*) ' =',ymin,ymax write (*,*) ' =',zmin,zmax 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 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) c write (*,'(5x,6f8.2,g11.3)') c & 1000.*Xc,1000.*Xp(jp),1000.*Yc,1000.*Yp(jp), c & 1000.*Zc,1000.*Zp(jp), c & sqrt(dd)*1000. EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo c write (*,100) nn,i1,i2,j1,j2,k1,k2,Xc,Yc,Zc,Amnew c 100 format(' Nn=',i6,3(2i4,2x),4f8.3) Return End