C ______________________________ Profiles for Galaxy Simulations C C INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' REAL AMprevious(Np),AMlast(Np),Daverg(0:500) Real*8 Xc,Yc,Zc DIMENSION Iaverg(0:500) REAL INPUT Character FileASCII*100 Character File1*100,File2*100,File3*100 EQUIVALENCE (Box,extras(100)) C................................................................... C Read data and open files File3 ='cooling_part.dat' write(*,*)' enter 3 digits of the snapshot' read(*,*) moment write(File1,'("FILES/PMcrda0.",i3.3,".DAT")')moment write(File2,'("FILES/PMcrs0a0.",i3.3,".DAT")')moment write(*,*) File1 write(*,*) File2 c File1 ='PMcrd.DAT' c File2 ='PMcrs0.DAT' CALL RDTAPE(File1,File2) iSnapshot =INT(AEXPN*10000.) write(FileASCII,'("RESULTS/Result.DM.",i5.5,".dat")') & iSnapshot OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result.SPH.",i5.5,".dat")') & iSnapshot OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') write(FileASCII,'("RESULTS/Result.CYL.",i5.5,".dat")') & iSnapshot OPEN(41,FILE=FileASCII,STATUS='UNKNOWN') OPEN(18,FILE='Stat.dat',STATUS='UNKNOWN',access='append') OPEN(20,FILE='particles.dat',STATUS='UNKNOWN') OPEN(30,FILE='particlesu.dat',STATUS='UNKNOWN',form='UNFORMATTED') OPEN(10,file=File3,form='UNFORMATTED',status='OLD') time = (age(AEXPN)-age(0.6))/1.e9 box = extras(100)! simulation box xfrac = extras(51) ! fraction of particles to be cooled cool_frac = extras(52) ! total cooling cfrac = extras(53) ! cooling per step aexpn_s = extras(54) ! start cooling aexpn_e = extras(55) ! and cooling i_cool_which = extras(57) ! cool every # step i_cool_all = extras(58) ! total # of cooling steps xm_part = extras(59) ! particle weight (h^-1 M_s) CALL ReadCool i_halo = 1 ! only one halo cooled n_cooled = iph_cool(i_halo,1) xm_part =1.e12/lspecies(Nspecies) aMassDisk = xm_part/hubble*n_cooled write (*,*) ' hubble =',hubble,' mass of disk=',aMassDisk write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN c Do i=37,41,2 c WRITE (i,100) HEADER,AEXPN,AEXP0,AMPLT, c + ASTEP,ISTEP,PARTW,EKIN, c + NROWC,NGRID,NRECL,Om0,Oml0,hubble, c + Ocurv,n_cooled,aMassDisk,time c EndDo 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=',I8,/ + 1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/ + 1x,' Omega_curvature=',F7.3,/ + 1x,' Ndisk =',i8,' Mass disk =',G11.3,' Age=',f7.3,'Gyr') c Box =INPUT(' Enter box size in comoving Mpc/h =') c Box =0.5 ! NG3109 c Box = 2 ! 'bulge' 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) ScaleH = 100.*hubble*sqrt(Om0/AEXPN**3+Oml0) ! H(a)=km/s/Mpc ScaleM = xm_part/hubble CALL ReadPnts(N,ScaleC,ScaleV,ScaleM) Call ChangePnts(N,Box,ScaleH) Call DumpPnts(N,Box,ScaleM) Call Statist(N,Box,ScaleM) c Call AngMom(N,30000,3.,Box,ScaleM) Call Profile(N,Box,ScaleM) c Call Projection(N,Box,0.) c Call Projection(N,Box,5.) END c-------------------------------------------------------------------- C Get basic statistics of the system SUBROUTINE Statist(N,Box,ScaleM) c-------------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Real*8 Lx,Ly,Lz,aL,V2,Lz1,Lz2,Lz3,Lz5,aM1,aM2,aM3,aM5 Real*8 Lxh,Lyh,Lzh,aLh,dRh,dR2h,Wh,Lhalo Real*8 Xcc,Ycc,Zcc,Vcx,Vcy,Vcz,dR,dR2,X,Y,Z Real*8 Lzh1,Lzh2,Lzh3,Lzh5,aMh1,aMh2,aMh3,aMh5 Real*8 R1,R2,R3,R5,Rh1,Rh2,Rh3,Rh5 Real*8 Xx2,Yy2,Zz2 DATA Vcx,Vcy,Vc,dR,dR2,V2/6*0.d+0/ DATA Lx,Ly,Lz,Lxh,Lyh,Lzh,Wh,dRh,dR2h/9*0.d+0/ DATA Lz1,Lz2,Lz3,Lz5,aM1,aM2,aM3,aM5/8*0.d+0/ DATA Lzh1,Lzh2,Lzh3,Lzh5,aMh1,aMh2,aMh3,aMh5/8*0.d+0/ DATA R1,R2,R3,R5,Rh1,Rh2,Rh3,Rh5/8*0.d+0/ Rad1 =1. Rad2 =2. Rad3 =3. Rad5 =5. Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. Xx2=0. ! find moments of inertia Yy2=0. Zz2=0. n_cooled = 0 Do jp=1,N X = (Xp(jp)- Xcc) *1000. ! scale to kpc inits Y = (Yp(jp)- Ycc) *1000. Z = (Zp(jp)- Zcc) *1000. aLx =Y*VZp(jp) -Z*VYp(jp) aLy =Z*VXp(jp) -X*VZp(jp) aLz =X*VYp(jp) -Y*VXp(jp) dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) If(iMask(jp).gt.0)Then n_cooled =n_cooled +1 Lx =Lx +aLx Ly =Ly +aLy Lz =Lz +aLz dR = dR + ds dR2= dR2+ dd If(ds.le.Rad5)Then Lz5 =Lz5 +aLz aM5 = aM5 + 1 R5 = R5 + ds If(ds.le.Rad3)Then Lz3 =Lz3 +aLz aM3 = aM3 + 1 R3 = R3 + ds If(ds.le.Rad2)Then Lz2 =Lz2 +aLz aM2 = aM2 + 1 R2 = R2 + ds If(ds.le.Rad1)Then Lz1 =Lz1 +aLz aM1 = aM1 + 1 R1 = R1 + ds EndIf EndIf EndIf EndIf Else W = iWeight(jp) Lxh =Lxh +aLx*W Lyh =Lyh +aLy*W Lzh =Lzh +aLz*W dRh = dRh + ds*W dR2h= dR2h+ dd*W Vcx = Vcx + VXp(jp) *W Vcy = Vcy + VYp(jp) *W Vcz = Vcz + VZp(jp) *W Wh = Wh + iWeight(jp) If(ds.le.Rad5)Then Lzh5 =Lzh5 +aLz*W aMh5 = aMh5 + W Rh5 = Rh5 + ds*W If(ds.le.Rad3)Then Lzh3 =Lzh3 +aLz*W aMh3 = aMh3 + W Rh3 = Rh3 + ds*W If(ds.le.Rad2)Then Lzh2 =Lzh2 +aLz*W aMh2 = aMh2 + W Rh2 = Rh2 + ds*W If(ds.le.Rad1)Then Lzh1 =Lzh1 +aLz*W aMh1 = aMh1 + W Rh1 = Rh1 + ds*W EndIf EndIf EndIf EndIf EndIf V2 = V2 + VXp(jp)**2 +VYp(jp)**2+VZp(jp)**2 EndDo Lx =Lx/n_cooled Ly =Ly/n_cooled Lz =Lz/n_cooled V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/n_cooled Lhalo = sqrt(Lxh**2 +Lyh**2 +Lzh**2) write (*,110) & dR,sqrt(dR2/n_cooled), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*n_cooled, & aL, & Lhalo,Lxh,Lyh,Lzh, & Lhalo/Wh,Lxh/Wh,Lyh/Wh,Lzh/Wh write (37,110) & dR,sqrt(dR2/n_cooled), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*n_cooled, & aL, & Lhalo,Lxh,Lyh,Lzh, & Lhalo/Wh,Lxh/Wh,Lyh/Wh,Lzh/Wh 110 format( & 5x,' Aver Distance Disk Particles(kpc)=',g11.4, & ' rms distance(kpc)=',g11.4,/ & 5x,' rms Velocity=',g11.4,' km/s',/ & 5x,' Disk Ang.Momentum: direction=',3g11.3,/ & 5x,' Disk Ang.Momentum: L/ =',g11.4,' km/s',/ & 5x,' Disk Ang.Momentum: L =',g11.4,' kpc*km/s',/ & 5x,' Disk Ang.Momentum: L/M =',g11.4,' kpc*km/s',/ & 5x,' Halo Ang.Momentum: L =',4g11.4,' kpc*km/s',/ & 5x,' Halo Ang.Momentum: L/M =',4g11.4,' kpc*km/s') write (*,120) Vcx/Wh, Vcy/Wh,Vcz/Wh,Wh write (37,120) Vcx/Wh, Vcy/Wh,Vcz/Wh,Wh 120 format( & 5x,' Aver Velocity of Dm Particles(km/s)=',3g11.3, & ' Total Weight=',g11.3) write (37,140) write (37,130) Rad1, & Lz1,Lz1/max(aM1,1.d+0),Lz1/max(R1,1.d-10), & R1/max(aM1,1.d+0),aM1, & Lzh1,Lzh1/max(aMh1,1.d+0),Lzh1/max(Rh1,1.d-10), & Rh1/max(aMh1,1.d+0),aMh1 write (37,130) Rad2, & Lz2,Lz2/max(aM2,1.d+0),Lz2/max(R2,1.d-10), & R2/max(aM2,1.d+0),aM2, & Lzh2,Lzh2/max(aMh2,1.d+0),Lzh2/max(Rh2,1.d-10), & Rh2/max(aMh2,1.d+0),aMh2 write (37,130) Rad3, & Lz3,Lz3/max(aM3,1.d+0),Lz3/max(R3,1.d-10), & R3/max(aM3,1.d+0),aM3, & Lzh3,Lzh3/max(aMh3,1.d+0),Lzh3/max(Rh3,1.d-10), & Rh3/max(aMh3,1.d+0),aMh3 write (37,130) Rad5, & Lz5,Lz5/max(aM5,1.d+0),Lz5/max(R5,1.d-10), & R5/max(aM5,1.d+0),aM5, & Lzh5,Lzh5/max(aMh5,1.d+0),Lzh5/max(Rh5,1.d-10), & Rh5/max(aMh5,1.d+0),aMh5 time = Age(AEXPN)-Age(0.60) write (18,150) AEXPN,time/1.e9,ISTEP,aL*n_cooled, & bar_angle*180./3.1415926, & bar_length,sqrt(1.-bar_ecc), & Rad1,Lz1/max(aM1,1.d+0),aM1, & Rad2,Lz2/max(aM2,1.d+0),aM2, & Rad5,Lz5/max(aM5,1.d+0),aM5 150 format(2f8.4,i4,g12.5,3f8.3,4(2x,f6.1,g11.3,g11.3)) 140 format(' Rad(kpc) Disk:Lz',4x,'Lz/M',2x,'km/s',1x, & 'kpc Npart',3x,' Halo: Lz',4x,'Lz/M',2x, & 'km/s',1x,'kpc Npart') 130 format(f7.2,2(2g11.3,2f7.2,g11.3,3x)) Return End C------------------------------------------ C write particles in a file SUBROUTINE DumpPnts(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Real*8 MassTot Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. MassTot =0. Do jp=1,N MassTot = MassTot + iWeight(jp) EndDo MassTot = MassTot * ScaleM i_halo = 1 ! only one halo cooled n_cooled = iph_cool(i_halo,1) write(20,100) write(20,110) N,n_cooled,age(AEXPN)-age(0.6),ISTEP, & ScaleM, ScaleM*n_cooled, MassTot write(20,115) c CALL GetDensity(Box) Do jp=1,N If(iMask(jp).gt.0)Then ! DM particle xx = (Xp(jp)-Xcc)*1000. yy = (Yp(jp)-Ycc)*1000. zz = (Zp(jp)-Zcc)*1000. write (20,130) xx,yy,zz, & VXp(jp) ,VYp(jp),VZp(jp) c xx*VYp(jp)-yy*VXp(jp) c & INT(iWeight(jp)+0.1) EndIf EndDo 100 format(' TotalPart N: Disk ', & T32,'Time',T42,'Step', & T48,'MassParticle',T62,'MassDisk', & T74,'MassTotal',/ & T32,'(yrs)',T48,'(Msun)',T62,'(Msun)', & T74,'(Msun)') 115 format(T3,'X(kpc) Y Z', & T38,'Vx(kms) Vy Vz') 110 format(i10,i7,T30,1P,g11.3,T41,i4, & T48,g11.4,T62,g11.4,T74,g11.4) 120 format(1P,3g13.5,2x,3g12.4,i3) 130 format(3g13.5,2x,4g12.4,i3) Return End C------------------------------------------ C Change positions: center, inclination SUBROUTINE ChangePnts(N,Box,ScaleH) C------------------------------------------ PARAMETER (pi =3.14159265) ! INCLUDE 'PMparameters.h' INCLUDE 'PMcool.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 Real*8 Xx,Yy,Zz,XY Ncenter = N/2 If(Ncenter.eq.0)write (*,*) ' Error!! No particles in ChangePnts' Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. 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 Rselect = 0.250 Xc0 = Xcc Yc0 = Ycc Zc0 = Zcc R2 = Rselect**2 Ncount = 0 dn =0. 300 nn = 0 dn0 =dn Xc = 0. Yc = 0. Zc = 0. Lx = 0. Ly = 0. Lz = 0. dR = 0. dR2 = 0. Vcx = 0. Vcy = 0. Vcz = 0. c Find center of mass and mean V Do jp=1,Ncenter dd = (Xp(jp) -Xc0)**2 +(Yp(jp) -Yc0)**2 +(Zp(jp) -Zc0)**2 If(dd.lt.R2)Then nn = nn + 1 Xc = Xc +Xp(jp) Yc = Yc +Yp(jp) Zc = Zc +Zp(jp) Vcx = Vcx + Vxp(jp) Vcy = Vcy + Vyp(jp) Vcz = Vcz + Vzp(jp) EndIf EndDo If(nn.eq.0)write (*,*) ' Error: no particles found. R=',Rselect Xc = Xc/nn Yc = Yc/nn Zc = Zc/nn Vcx = Vcx/nn Vcy = Vcy/nn Vcz = Vcz/nn xshift = Xcc - Xc ! displacements to bring yshift = Ycc - Yc ! the center to {Xcc,Ycc,Zcc} zshift = Zcc - Zc c write (*,100) nn,Rselect*1000.,nn/Rselect**3,Xc,Yc,Zc, c & 1000.*xshift,1000.*yshift,1000.*zshift, c & Vcx,Vcy,Vcz 100 format('--- Center: n=',i8,f7.2,g11.2, & ' X(Mpc)=',3g12.6,/ & 25x,'shift(kpc) =',3g11.5,/ & 25x,' average Velocity(km/s)=',3g11.4) Xc0 = Xc Yc0 = Yc Zc0 = Zc Ncount = Ncount +1 Rselect = Rselect/1.1 R2 = Rselect**2 dn = nn/Rselect**3 If(Ncount.lt.180.and.Rselect.gt.0.2e-3.and.nn.gt.20)goto300 write (37,100) nn,Rselect*1000.,nn/Rselect**3,Xc,Yc,Zc, & 1000.*xshift,1000.*yshift,1000.*zshift, & Vcx,Vcy,Vcz n_cooled =0 Do jp=1,N Xp(jp) =Xp(jp) +xshift Yp(jp) =Yp(jp) +yshift Zp(jp) = Zp(jp) +zshift VXp(jp) = VXp(jp) +ScaleH*(Xp(jp) -Xcc)-Vcx VYp(jp) = VYp(jp) +ScaleH*(Yp(jp) -Ycc)-Vcy VZp(jp) = VZp(jp) +ScaleH*(Zp(jp) -Zcc)-Vcz 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 If(iMask(jp).gt.0)Then ! cooled particle n_cooled = n_cooled +1 Lx =Lx +(Yp(jp)-Ycc)*VZp(jp) & -(Zp(jp)-Zcc)*VYp(jp) Ly =Ly +(Zp(jp)-Zcc)*VXp(jp) & -(Xp(jp)-Xcc)*VZp(jp) Lz =Lz +(Xp(jp)-Xcc)*VYp(jp) & -(Yp(jp)-Ycc)*VXp(jp) dd = (Xp(jp)-Xcc)**2 + & (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd V2 = V2 + VXp(jp)**2+VYp(jp)**2+VZp(jp)**2 EndIf xmin =min(xmin,Xp(jp)) ymin =min(ymin,Yp(jp)) zmin =min(zmin,Zp(jp)) xmax =max(xmax,Xp(jp)) ymax =max(ymax,Yp(jp)) zmax =max(zmax,Zp(jp)) EndDo Lx =Lx/n_cooled Ly =Ly/n_cooled Lz =Lz/n_cooled V2 =sqrt(V2/n_cooled) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/n_cooled write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/n_cooled), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*n_cooled 110 format(' --- All particles: Xmin/max=',2g11.4,/ & 30x,2g11.4,/30x,2g11.4/, & 25x,' Aver Distance Disk Particles(kpc)=',g11.4, & ' rms distance(kpc)=',g11.4,/ & 25x,' rms Velocity=',g11.4,' km/s',/ & 25x,' Disk Ang.Momentum: direction=',3g11.4,/ & 25x,' Disk Ang.Momentum: L/ =',g11.4,' km/s',/ & 25x,' Disk Ang.Momentum: L =',g12.5,' kpc*km/s') c Rotate the system of coordinates return c If(Lz.le.0.5*aL)Then write (*,*) ' z-component of angular momentum', & ' is too small. Something is very wrong.' c 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.5e-1)stop if(abs(prod2).gt.1.5e-1)stop if(abs(prod3).gt.1.5e-1)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) +ey1*Vyp(jp) +ez1*Vzp(jp) VYnew = ex2*Vxp(jp) +ey2*Vyp(jp) +ez2*Vzp(jp) VZnew = ex3*Vxp(jp) +ey3*Vyp(jp) +ez3*Vzp(jp) 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 Vyp(jp) = VYnew Vzp(jp) = VZnew 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)) If(iMask(jp).gt.0)Then Lx =Lx +(Yp(jp)-Ycc)*VZp(jp) & -(Zp(jp)-Zcc)*VYp(jp) Ly =Ly +(Zp(jp)-Zcc)*VXp(jp) & -(Xp(jp)-Xcc)*VZp(jp) Lz =Lz +(Xp(jp)-Xcc)*VYp(jp) & -(Yp(jp)-Ycc)*VXp(jp) dd = (Xp(jp)-Xcc)**2 + & (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd EndIf V2 = V2 + VXp(jp)**2 +VYp(jp)**2+VZp(jp)**2 EndDo Lx =Lx/n_cooled Ly =Ly/n_cooled Lz =Lz/n_cooled V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/n_cooled write (*,*) '--- After rotation of the coordinates' write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/n_cooled),V2, & Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*n_cooled Return End C------------------------------------------ C Distribution of the angular momentum C of disk and halo particles SUBROUTINE AngMom(N,N1,R1,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' PARAMETER ( Nang = 50, Ntang = 2*Nang+1 ) PARAMETER ( aMax = 4000. ) Real*8 Lzd(-Nang:Nang),Lzd1(-Nang:Nang),Lzd2(-Nang:Nang) Integer Mzd(-Nang:Nang),Mzd1(-Nang:Nang),Mzd2(-Nang:Nang) Real*8 Lzh(-Nang:Nang),Lzh1(-Nang:Nang),Lzh2(-Nang:Nang) Integer Mzh(-Nang:Nang),Mzh1(-Nang:Nang),Mzh2(-Nang:Nang) DATA Lzd/Ntang*0.d+0/,Lzd1/Ntang*0.d+0/,Lzd2/Ntang*0.d+0/ DATA Lzh/Ntang*0.d+0/,Lzh1/Ntang*0.d+0/,Lzh2/Ntang*0.d+0/ DATA Mzd/Ntang*0/,Mzd1/Ntang*0/,Mzd2/Ntang*0/ DATA Mzh/Ntang*0/,Mzh1/Ntang*0/,Mzh2/Ntang*0/ Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. hz = aMax/Nang Do jp=1,N X = (Xp(jp)- Xcc) *1000. ! scale to kpc inits Y = (Yp(jp)- Ycc) *1000. Z = (Zp(jp)- Zcc) *1000. aLx =Y*VZp(jp) -Z*VYp(jp) aLy =Z*VXp(jp) -X*VZp(jp) aLz =X*VYp(jp) -Y*VXp(jp) dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) ind = INT(aLz/hz+10000.) -10000 ind = Min(Max(ind,-Nang),Nang) If(iMask(jp).gt.0)Then Lzd(ind) =Lzd(ind) +aLz Mzd(ind) =Mzd(ind) + 1 If(ds.le.R1)Then Lzd2(ind) =Lzd2(ind) +aLz Mzd2(ind) =Mzd2(ind) + 1 EndIf If(jp.le.N1)Then Lzd1(ind) =Lzd1(ind) +aLz Mzd1(ind) =Mzd1(ind) + 1 EndIf Else W = iWeight(jp) Lzh(ind) =Lzh(ind) +aLz*W Mzh(ind) =Mzh(ind) + W If(ds.le.R1)Then Lzh2(ind) =Lzh2(ind) +aLz*W Mzh2(ind) =Mzh2(ind) + W EndIf EndIf EndDo write (*,100)N1,R1,R1 write (37,100)N1,R1,R1 100 format(5x,'Disk',12x,'n Npart', & ' Npart', & ' Npart', &' Npart', & ' Npart', & ' Npart') Do i=-Nang,Nang a = i*hz aL = Lzd(i)/Max(Mzd(i),1) aL2= Lzd2(i)/Max(Mzd2(i),1) aL1= Lzd1(i)/Max(Mzd1(i),1) bL = Lzh(i)/Max(Mzh(i),1) bL2= Lzh2(i)/Max(Mzh2(i),1) write(*,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i), & bL,Mzh(i),bL2,Mzh2(i) write(37,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i), & bL,Mzh(i),bL2,Mzh2(i) EndDo 110 format(f7.1,3(f7.1,i6),5x,f8.1,i8,2(f8.1,i6)) Return End C------------------------------------------ C SUBROUTINE Projection(N,Box,Alpha) c Alpha is angle between plain of the disk and c plain of projection c Zold Znew c ` / Disk (Xold) | c ` / | c ` / Alpha | c==>>> ------------------------------ Xnew ==> -------------------------------Xnew c / | c / | c c Znew is the distance in the plain of projects c Vlos = V_x,new c Average for two observers 90 degrees away C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Parameter (pi180 = 180./3.1415926535) Parameter (pi = 3.1415926535) EQUIVALENCE (Ndisk,extras(90)) drlog =0.05 d_lim = 0.100 cosa = COS(Alpha/pi180) sina = SIN(Alpha/pi180) write (*,*) ' Angle =',Alpha, ' cos/sin=',cosa,sina Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. Do i =-Nrad,Nrad Nradp(i) = 0 Radp(i) = 0. Vradp(i) = 0. Vrad2p(i) = 0. EndDo Do jp=1,Ndisk X = (Xp(jp)- Xcc) *1000. ! scale to kpc inits Y = (Yp(jp)- Ycc) *1000. Z = (Zp(jp)- Zcc) *1000. c x-z coordinates z_new = X*sina +Z*cosa ! z-coord in projection If(abs(z_new).le.d_lim)Then ! inside the strip ind = INT(LOG10(ABS(Y))/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Vradial = -VXp(jp)*cosa +VZp(jp)*sina If(Y.le.0.) Vradial =-Vradial Nradp(ind) = Nradp(ind) + 1 Radp(ind) = Radp(ind) + abs(Y) Vradp(ind) = Vradp(ind) + Vradial Vrad2p(ind) = Vrad2p(ind) + Vradial**2 EndIf c y-z coordinayes z_new = Y*sina +Z*cosa ! z-coord in projection If(abs(z_new).le.d_lim)Then ! inside the strip ind = INT(LOG10(ABS(X))/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Vradial = VYp(jp)*cosa -VZp(jp)*sina If(X.le.0.) Vradial =-Vradial Nradp(ind) = Nradp(ind) + 1 Radp(ind) = Radp(ind) + abs(X) Vradp(ind) = Vradp(ind) + Vradial Vrad2p(ind) = Vrad2p(ind) + Vradial**2 EndIf EndDo V1 =0 write(37,90) Alpha 90 format(1x,10('-'),'Projection along major axis',10('-'),/ & ' Rout(kpc) Rdisk ', & ' Vrot sig_v_rot Npart Alpha=',f7.2) Do i=-Nrad+1,Nrad-1 rad_out = 10.**((i+1)*drlog) nn = Nradp(i) Din = Din + nn n0 = max(nn,1) Raddisk = Radp(i)/n0 If(nn.eq.0)Raddisk=rad_out Vradial = Vradp(i)/n0 Vradial2=sqrt(max(Vrad2p(i)/n0-Vradial**2,1.e-10)) If(nn.gt.20)Then write(37,100)rad_out,Raddisk, & Vradial/cosa,Vradial2/cosa,nn EndIf EndDo 100 format(2g10.3,2g11.3,i7) Return End C------------------------------------------ C SUBROUTINE Profile(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Parameter (pi180 = 180./3.1415926535) Parameter (pi = 3.1415926535) drlog =0.075 Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. c DM profile --------------------- spherical shells Do jp=1,N If(iMask(jp).lt.0)Then ! not cooled particle W = iWeight(jp) X = (Xp(jp)- Xcc) *1000. ! scale to kpc inits Y = (Yp(jp)- Ycc) *1000. Z = (Zp(jp)- Zcc) *1000. dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) Vtot = (VXp(jp)**2+VYp(jp)**2+VZp(jp)**2) Vradial = (VXp(jp)*X+VYp(jp)*Y+VZp(jp)*Z) / & max(ds,1.e-5) Vtan2 = Vtot -Vradial**2 ind = INT(LOG10(ds)/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Nradh(ind) = Nradh(ind) + W Radh(ind) = Radh(ind) + ds*W Vfih2(ind) = Vfih2(ind) + Vtan2*W Vradh(ind) = Vradh(ind) + Vradial*W Vradh2(ind) = Vradh2(ind) + Vradial**2*W EndIf EndDo V1 = 0. Din = 0. write(37,90) 90 format(1x,20('-'),' DM profile',/ & ' Rout Rhalo Vchalo ', & ' Density ', & ' Vrad sigRadial V_tan Beta Nhalo') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 = 4.1888*rad_out**3 wn = Nradh(i) nn =wn /iWeight(1) Din = Din + wn w0 = max(wn,1.e-9) densD = wn/(V1-V0)*ScaleM*1.e-9 Radhalo = Radh(i)/w0 If(wn.eq.0)Radhalo=rad_out Vtan = sqrt(Vfih2(i)/w0) Vradial = Vradh(i)/w0 Vradial2 =sqrt(max(Vradh2(i)/w0-Vradial**2,1.e-10)) Vtot = sqrt((Vradh2(i)+Vfih2(i))/w0) Vcirc = 2.08e-3*sqrt(max(Din*ScaleM/rad_out,1.e-10)) beta = 1. -Vtan**2/2./(Vradial2**2+Vradial**2) If(nn.gt.4)Then write(37,100)rad_out,Radhalo,Vcirc,densD, & Vradial,Vradial2,Vtan,beta,nn,INT(Din/iWeight(1)) EndIf EndDo 100 format(2g10.3,f7.2, & g11.3,1x,4g11.4,2x,i7,i9) ! disk statistics: spherical and cylindrical c Do i=-Nrad,Nrad c Nradp(i) =0 c EndDo i_halo = 1 ! only one halo cooled n_cooled = iph_cool(i_halo,1) Do i=1,n_cooled jp = ip_cool(i) X = (Xp(jp)- Xcc) *1000. ! scale to kpc inits Y = (Yp(jp)- Ycc) *1000. Z = (Zp(jp)- Zcc) *1000. aLx =Y*VZp(jp) -Z*VYp(jp) aLy =Z*VXp(jp) -X*VZp(jp) aLz =X*VYp(jp) -Y*VXp(jp) dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) rd = X**2 + Y**2 rs = sqrt(rd) Vrotat = aLz/max(rs,1.e-5) Vradial = (VXp(jp)*X+VYp(jp)*Y+VZp(jp)*Z) / & max(ds,1.e-5) Vvert = VZp(jp) ind = INT(LOG10(ds)/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) ! sperical shells Nradd(ind) = Nradd(ind) + 1 Radd(ind) = Radd(ind) + ds Vfi(ind) = Vfi(ind) + Vrotat c Vfi2(ind) = Vfi2(ind) + Vrotat**2 c Vrad(ind) = Vrad(ind) + Vradial c Vrad2(ind) = Vrad2(ind) + Vradial**2 c Vzz(ind) = Vzz(ind) + Vvert c Vzz2(ind) = Vzz2(ind) + Vvert**2 ! projection on the plane ind = INT(LOG10(rs)/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Nradp(ind) = Nradp(ind) + 1 Radp(ind) = Radp(ind) + rs Vfip(ind) = Vfip(ind) + Vrotat Vfi2p(ind) = Vfi2p(ind) + Vrotat**2 Vradp(ind) = Vradp(ind) + Vradial Vrad2p(ind) = Vrad2p(ind) + Vradial**2 Vzzp(ind) = Vzzp(ind) + Vvert Vzz2p(ind) = Vzz2p(ind) + Vvert**2 ! projection along wedges inf = -1 If(abs(Y)/rs.lt.0.5) inf =1 ! along the bar If(abs(X)/rs.lt.0.5) inf =2 ! perpendic bar IF(inf.gt.0)Then Nradb(ind,inf) = Nradb(ind,inf) + 1 Radb(ind,inf) = Radb(ind,inf) + rs Vfib(ind,inf) = Vfib(ind,inf) + Vrotat Vfi2b(ind,inf) = Vfi2b(ind,inf) + Vrotat**2 Vradb(ind,inf) = Vradb(ind,inf) + Vradial Vrad2b(ind,inf) = Vrad2b(ind,inf) + Vradial**2 Vzzb(ind,inf) = Vzzb(ind,inf) + Vvert Vzz2b(ind,inf) = Vzz2b(ind,inf) + Vvert**2 EndIf EndDo ! print statistics ------------------------------ V1 = 0. Din = 0. write(39,92) 92 format(1x,20('-'),'All components: Spherical shells',/ & ' Rout Rdisk Vchalo Vcdisk Vcbulge Vctot', & ' Vrot_disk',6x,'Nhalo',11x,'Ndisk',14x,'Nbulge') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 = 4.1888*rad_out**3 wn = Nradh(i) nn =wn /iWeight(1) Din = Din + wn w0 = max(wn,1.e-9) Radhalo = Radh(i)/w0 If(wn.eq.0)Radhalo=rad_out nd = Nradd(i) Dind = Dind + nd n0 = max(nd,1) Raddisk = Radd(i)/n0 Vrotat =Vfi(i)/n0 If(nd.eq.0)Raddisk=rad_out nbl = Nradg(i) Dinbl = Dinbl + nbl n0 = max(nbl,1) Radbulge = Radg(i)/n0 If(nbl.eq.0)Radbulge=rad_out Vcirc = 2.08e-3*sqrt(max((Din+Dind+Dinbl)*ScaleM/rad_out, & 1.e-10)) VcD = 2.08e-3*sqrt(max(( Dind)*ScaleM/rad_out,1.e-10)) VcH = 2.08e-3*sqrt(max((Din )*ScaleM/rad_out,1.e-10)) VcBl = 2.08e-3*sqrt(max((Dinbl )*ScaleM/rad_out,1.e-10)) If(nn.gt.10)Then write(39,105)rad_out,Raddisk,VcH,VcD,VcBl,Vcirc, & Vrotat,INT(wn),INT(Din), & nd,INT(Dind),nbl,INT(Dinbl) EndIf EndDo 105 format(2g10.3,f8.2,4f7.2,i8,2x,i8,2x,i6,4(2x,i8)) ! cylindrical shells: print --------------------- V1 = 0. Din = 0. write(41,94) 94 format(1x,20('-'),'Cylindrical shells profile: Disk and Bulge',/ & ' Rout Rdisk Vrot sig_rotat sig_z sig_rad', & ' SurfDisk SurfBulge Bar: Vrot sig_rot Surf', & ' off bar: Vrot sig_rot Surf') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 = 3.141592*rad_out**2 nd = Nradp(i) Dind = Dind + nd n0 = max(nd,1) Raddisk = Radp(i)/n0 If(nd.eq.0)Raddisk=rad_out Vrotat =Vfip(i)/n0 ! all disk V2p =sqrt(max(Vfi2p(i)/n0-Vrotat**2,1.e-10)) Vwz =sqrt(Vzz2p(i)/n0) Vradial =sqrt(Vrad2p(i)/n0) Surf = nd/(V1-V0) SurfBl = Nradu(i)/(V1-V0) nb = Nradb(i,1) ! along bar +- 30 degrees n0b = max(nb,1) VrotatB =Vfib(i,1)/n0b V2B =sqrt(max(Vfi2b(i,1)/n0b-VrotatB**2,1.e-10)) SurfB = nb/(V1-V0)*3. ! correction for 30degrees nb = Nradb(i,2) ! off bar +- 30 degrees n0b = max(nb,1) VrotatN =Vfib(i,2)/n0b V2N =sqrt(max(Vfi2b(i,2)/n0b-VrotatN**2,1.e-10)) SurfN = nb/(V1-V0)*3. vvv =sqrt(Vfi2b(i,2)/n0b) c write (*,'(4i6,6g11.3)') nd,n0,Nradb(i,1),Nradb(i,2), c & Radp(i),rad_out,VrotatB ,V2B,VrotatN ,V2N If(nd.gt.10)Then write(41,107)rad_out,Raddisk, & Vrotat ,V2p,Vwz,Vradial,Surf,SurfBl, & VrotatB,V2B,SurfB, & VrotatN,V2N,SurfN EndIf EndDo 107 format(2g10.3,4f7.2,2g11.3,3x,2f7.2,g11.3,3x,2f7.2,g11.3) Return End C------------------------------------------ C Read cooled particles, set Mask SUBROUTINE ReadCool C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Real xmax(2),xmin(2),ymax(2),ymin(2), + zmax(2),zmin(2),vmax(2),vmin(2) DATA xmax,ymax,zmax,vmax/8*0.e-0/ DATA xmin,ymin,zmin,vmin/8*0.e+15/ read(10) n_halo IF(n_halo .gt. nh_cool) THEN write(*,*) 'n_halo .gt. nh_cool', n_halo, nh_cool write(*,*) 'increase nh_cool in a_cooling.h' STOP ENDIF kkL = 0 ! last index of the previous halo DO k = 1, n_halo read(10) (iph_cool(k,i),i =1,3),(halo_prop(k,kk),kk = 1,nh_prop) write(*,*) 'halo # :',k, ' iph_cool(k,i) = ', + (iph_cool(k,i),i =1,3) IF(iph_cool(k,1) .gt. n_cool) THEN write(*,*) 'iph_cool(k,1) .gt. nh_cool',iph_cool(k,1),n_cool write(*,*) 'increase n_cool in a_cooling.h' STOP ENDIF kkF = kkL + 1 ! read First particle if (kkF .NE. iph_cool(k,2)) write(*,*) 'ERROR in kkF' kkL = kkF + iph_cool(k,1) - 1 ! read Last particle if (kkL .NE. iph_cool(k,3)) write(*,*) 'ERROR in kkL' read(10) (ip_cool(kk), kk = iph_cool(k,2),iph_cool(k,3)) ENDDO close (10) i_halo = 1 ! only one halo cooled n_cooled = iph_cool(i_halo,1) write(*,*) 'n_cooled = ', n_cooled Do i=1,Np iMask(i) = -1 EndDo DO k = 1, n_cooled i = ip_cool(k) iMask(i) = k ENDDO Return End C------------------------------------------ C Read particles SUBROUTINE ReadPnts(N,ScaleC,ScaleV,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Integer Xdist(0:25),Ydist(0:25),Zdist(0:25) Real*8 TotWeight Box = ScaleC*NGRID Radius2 = Radius**2 N = 0 ! Number of particles Nsm = 0 TotWeight =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 c Do i=0,25 c Xdist(i) =0 c Ydist(i) =0 c Zdist(i) =0 c 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 cc N_particles = N_particles/3 !!!!! Npages = (N_particles -1)/NPAGE +1 N_in_last=N_particles -NPAGE*(Npages-1) EndIf write (*,*) ' Pages=',Npages,' Species=',Nspecies c 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 c 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) ! particles weight iWeight(Ipart) TotWeight = TotWeight + iWeight(N) 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)) c ii =Min(INT(X/Box*25.),25) c jj =Min(INT(Y/Box*25.),25) c kk =Min(INT(Z/Box*25.),25) c Xdist(ii) =Xdist(ii) +1 c Ydist(jj) =Ydist(jj) +1 c 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 (*,'(" Total weight =",1P,g11.3)')TotWeight write (*,'(" Total Mass =",1P,g11.3)')TotWeight*ScaleM c write (*,'(i3,2x,i9,2x,i9,2x,i9)') c & (i,Xdist(i),Ydist(i),Zdist(i),i=0,25) Return End c ----------------------------- Function Age(aexp) c ----------------------------- c age for LCDM model PARAMETER (Om0 =0.3) PARAMETER (tconst =9.31e+9) x = ((1.-Om0)/Om0)**0.333333*aexp Age = tconst /sqrt(1.-Om0)*(log(x**1.5+sqrt(1.+x**3))) Return End C ______________________________ Assign initial values C C BLOCK DATA INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' DATA Nradd/Nrad2*0/,Vfi/Nrad2*0./,Vrad/Nrad2*0./, & Vzz/Nrad2*0./,Vfi2/Nrad2*0./,Vrad2/Nrad2*0./, & Vzz2/Nrad2*0./,Radd/Nrad2*0./ DATA Nradb/Nbt*0/,Vfib/Nbt*0./,Vradb/Nbt*0./, & Vzzb/Nbt*0./,Vfi2b/Nbt*0./,Vrad2b/Nbt*0./, & Vzz2b/Nbt*0./,Radb/Nbt*0./ DATA Nradp/Nrad2*0/,Vfip/Nrad2*0./,Vradp/Nrad2*0./, & Vzzp/Nrad2*0./,Vfi2p/Nrad2*0./,Vrad2p/Nrad2*0./, & Vzz2p/Nrad2*0./,Radp/Nrad2*0./ DATA Nradu/Nrad2*0/ DATA Nradh/Nrad2*0/,Vfih/Nrad2*0./,Vradh/Nrad2*0./, & Vfih2/Nrad2*0./,Vradh2/Nrad2*0./,Radh/Nrad2*0./ DATA Dens/Np*0./ End C--------------------------------------------------- SUBROUTINE RDTAPE(File1,File2) C--------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' Character File1*100,File2*100 C Open file on a tape OPEN(UNIT=9,FILE=File1, + FORM='UNFORMATTED', STATUS='UNKNOWN') C Read control information C and see whether it has proper C structure READ (9,err=10,end=10) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, + NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5 + ,Ocurv,extras Ocurv =0. WRITE (*,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv WRITE (16,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv 100 FORMAT(1X,'Header=>',A45,/ + 1X,' A=',F8.5,' 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(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!/4 !for athena c write (*,*) ' Open direct access file' OPEN(UNIT=21,FILE=File2,ACCESS='DIRECT', + STATUS='UNKNOWN',RECL=NACCES) REWIND 9 RETURN END C-------------------------------------------- C Read current PAGE of particles from disk C NRECL - length of ROW block in words C-------------------------------------------- SUBROUTINE GETROW(IROW,Ifile) C-------------------------------------------- c INCLUDE 'Statparam.h' INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' READ (20+Ifile,REC=IROW) RECDAT RETURN END C-------------------------------------------------------------- C Make linker lists of particles in each cell SUBROUTINE GetDensity(Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' EQUIVALENCE (Ndisk,extras(90)) Cell = Box/Nblink dStep = 15.e-6 Do i=1,5 RadDens(i) =dStep*i**2 EndDo write (*,*) ' GetDens: Cell=',Cell,dStep CALL List(Box) C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (jp,Density) Do jp=1,Ndisk CALL Neib(Xp(jp),Yp(jp),Zp(jp),Density) Dens(jp) = Density EndDo Return End C-------------------------------------------------------------- C Make linker lists of particles in each cell SUBROUTINE List(Box) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Do i=1,Ndisk Lst(i)=-1 EndDo Do k=Nmlink,Nblink Do j=Nmlink,Nblink Do i=Nmlink,Nblink Label(i,j,k)=0 EndDo EndDo EndDo Do jp=1,Ndisk i=INT(Xp(jp)/Cell) j=INT(Yp(jp)/Cell) k=INT(Zp(jp)/Cell) i=MIN(MAX(Nmlink,i),Nblink) j=MIN(MAX(Nmlink,j),Nblink) k=MIN(MAX(Nmlink,k),Nblink) 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; a0 -its weight C SUBROUTINE Neib(Xc,Yc,Zc,Density) C-------------------------------------------------------------- INCLUDE 'PMparameters.h' INCLUDE 'PMcool.h' DIMENSION iDen(5) C Initiate counters Radius2= RadDens(5)**2 Radius = RadDens(5) Do i=1,5 iDen(i) =0 EndDo c write (*,*) '==== Neib:',Xc,Yc,Zc,dStep c limits for Label i1=INT((Xc-Radius)/Cell) j1=INT((Yc-Radius)/Cell) k1=INT((Zc-Radius)/Cell) i1=MIN(MAX(Nmlink,i1),Nblink) j1=MIN(MAX(Nmlink,j1),Nblink) k1=MIN(MAX(Nmlink,k1),Nblink) i2=INT((Xc+Radius)/Cell) j2=INT((Yc+Radius)/Cell) k2=INT((Zc+Radius)/Cell) i2=MIN(MAX(Nmlink,i2),Nblink) j2=MIN(MAX(Nmlink,j2),Nblink) k2=MIN(MAX(Nmlink,k2),Nblink) c write (*,'(" limits=",6i5)')i1,i2,j1,j2,k1,k2 C Look for neibhours 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 disd = sqrt(sqrt(dd)/dStep) ind = min(INT(disd) +1,5) iDen(ind) = iDen(ind) +1 c write (*,'(6x,i9,g11.4,3f9.5,2i5)') c & jp,sqrt(dd)*1.e6, c & (Xp(jp)-0.5)*1.e3,(Yp(jp)-0.5)*1.e3, c & (Zp(jp)-0.5)*1.e3, c & ind,iDen(ind) EndIf jp =Lst(jp) GoTo10 EndIf EndDo EndDo EndDo nn =0 Do i=1,5 nn =nn +iDen(i) If(nn.ge.5)Then Density = nn/RadDens(i)**3 goto 20 EndIf EndDo Density = nn/RadDens(5)**3 20 Density =Density/4.188e9 ! number of particles per kpc**3 c write(*,'(5x,"Dens=",g11.3,6i8)')Density,iDen,nn Return End