C ______________________________ Profiles for Galaxy Simulations C C INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' REAL AMprevious(Np),AMlast(Np),Daverg(0:500) Real*8 Xc,Yc,Zc DIMENSION Iaverg(0:500) REAL INPUT Character FileASCII*50 Character File1*100,File2*100,directory*3 EQUIVALENCE (Ndisk,extras(90)) c aMassDisk =6.e+10 ! M1/2 aMassDisk =4.e+8/0.7 ! MOVIE mass of disk in real Msun aMassDisk =4.80e+10 ! MW_Large mass of disk in real Msun C................................................................... C Read data and open files OPEN(18,FILE='potentP.dat',STATUS='UNKNOWN',form='UNFORMATTED') OPEN(19,FILE='EnergyH.dat',STATUS='UNKNOWN') OPEN(20,FILE='EnergyD.dat',STATUS='UNKNOWN') write (*,*) ' RDTAPE is done' c Box =INPUT(' Enter box size in comoving Mpc/h =') Box =2.0 AEXPN = 0.8840 hubble = 0.7 Ndisk = 1299997 Om0 = 0.3 Oml0 = 0.7 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 if(Ndisk.le.0)write (*,*) ' wrong number of disk particles:',Ndisk if(Ndisk.le.0)Stop ScaleM = aMassDisk/Ndisk CALL ReadPnts(N,ScaleC,ScaleV,ScaleM) Call ChangePnts(N,Box,ScaleH) c Call DumpPnts(N,Box,ScaleM) c Call Statist(N,Box,ScaleM) c Call AngMom(N,30000,3.,Box,ScaleM) c Call AngMom(N,60000,3.,Box,ScaleM) END C------------------------------------------ C Get statistics of the system SUBROUTINE Statist(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.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 EQUIVALENCE (Ndisk,extras(90)) 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 = 0. !Box/2. Ycc = 0. !Box/2. Zcc = 0. !Box/2. Do ip=1,N X = (Xp(ip)- Xcc) *1000. ! scale to kpc inits Y = (Yp(ip)- Ycc) *1000. Z = (Zp(ip)- Zcc) *1000. aLx =Y*VZp(ip) -Z*VYp(ip) aLy =Z*VXp(ip) -X*VZp(ip) aLz =X*VYp(ip) -Y*VXp(ip) dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) If(ip.le.Ndisk)Then 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(ip) Lxh =Lxh +aLx*W Lyh =Lyh +aLy*W Lzh =Lzh +aLz*W dRh = dRh + ds*W dR2h= dR2h+ dd*W Vcx = Vcx + VXp(ip) *W Vcy = Vcy + VYp(ip) *W Vcz = Vcz + VZp(ip) *W Wh = Wh + iWeight(ip) 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(ip)**2 +VYp(ip)**2+VZp(ip)**2 EndDo Lx =Lx/Ndisk Ly =Ly/Ndisk Lz =Lz/Ndisk V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/Ndisk Lhalo = sqrt(Lxh**2 +Lyh**2 +Lzh**2) write (*,110) & dR,sqrt(dR2/Ndisk), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*Ndisk, & aL, & Lhalo,Lxh,Lyh,Lzh, & Lhalo/Wh,Lxh/Wh,Lyh/Wh,Lzh/Wh write (17,110) & dR,sqrt(dR2/Ndisk), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*Ndisk, & 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 (17,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 (17,140) write (17,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 (17,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 (17,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 (17,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) time,ISTEP,aL*Ndisk, & Rad1,Lz1/max(aM1,1.d+0),INT(Am1), & Rad2,Lz2/max(aM2,1.d+0),INT(Am2), & Rad3,Lz3/max(aM3,1.d+0),INT(Am3), & Rad5,Lz5/max(aM5,1.d+0),INT(Am5) 150 format(g11.3,i4,g12.4,4(f6.1,f7.2,i7)) 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 Change positions: center, inclination SUBROUTINE DumpPnts(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Real*8 MassTot EQUIVALENCE (Ndisk,extras(90)) Xcc = 0. !Box/2. Ycc = 0. !Box/2. Zcc = 0. !Box/2. MassTot =0. Do ip=1,N MassTot = MassTot + iWeight(ip) EndDo MassTot = MassTot * ScaleM write(20,100) write(20,110) N,Ndisk,age(AEXPN)-age(0.6),ISTEP, & ScaleM, ScaleM*Ndisk, MassTot write(30) N,Ndisk,age(AEXPN)-age(0.6),ISTEP, & ScaleM, ScaleM*Ndisk, MassTot Do ip=1,N write (20,120) (Xp(ip)-Xcc)*1000.,(Yp(ip) -Ycc)*1000., & (Zp(ip)-Zcc)*1000., & VXp(ip) ,VYp(ip),VZp(ip), & INT(iWeight(ip)+0.1) EndDo write (30) ((Xp(ip)-Xcc)*1000.,(Yp(ip) -Ycc)*1000., & (Zp(ip)-Zcc)*1000., & VXp(ip) ,VYp(ip),VZp(ip), & INT(iWeight(ip)+0.1),ip=1,N) 100 format(' TotalPart DiskParticles', & T32,'Time',T42,'Step', & T48,'MassParticle',T62,'MassDisk', & T74,'MassTotal',/ & T32,'(Gyrs)',T48,'(Msun)',T62,'(Msun)', & T74,'(Msun)') 110 format(i10,i10,T30,1P,g11.3,T41,i4,T48,g11.4,T62,g11.4,T74,g11.4) 120 format(1P,3g12.4,2x,3g12.4,i3) Return End C------------------------------------------ C Change positions: center, inclination SUBROUTINE ChangePnts(N,Box,ScaleH) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.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 EQUIVALENCE (Ndisk,extras(90)) c Ncenter = Min(1500000,N) ! find center using central particles only Ncenter = N ! find center If(Ncenter.eq.0)write (*,*) ' Error!! No particles in ChangePnts' Xcc = 0. !Box/2. Ycc = 0. !Box/2. Zcc = 0. !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 nn = 0 c Find center of mass and mean V Do ip=1,Ncenter If(Jp(ip).le.15000)Then nn = nn +1 Xc = Xc +Xp(ip) Yc = Yc +Yp(ip) Zc = Zc +Zp(ip) Vcx = Vcx + Vxp(ip) Vcy = Vcy + Vyp(ip) Vcz = Vcz + Vzp(ip) EndIf EndDo Ncenter = nn 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,'shift(kpc) =',3g12.5,/ & 30x,' average Velocity(km/s)=',3g11.4) Do ip=1,N Xp(ip) =Xp(ip) +xshift Yp(ip) =Yp(ip) +yshift Zp(ip) = Zp(ip) +zshift VXp(ip) = VXp(ip) +ScaleH*(Xp(ip) -Xcc) VYp(ip) = VYp(ip) +ScaleH*(Yp(ip) -Ycc) VZp(ip) = VZp(ip) +ScaleH*(Zp(ip) -Zcc) c IF(Xp(ip).lt.-Box/2.)Xp(ip)=Xp(ip) +Box c IF(Yp(ip).lt.-Box/2.)Yp(ip)=Yp(ip) +Box c IF(Zp(ip).lt.-Box/2.)Zp(ip)=Zp(ip) +Box c IF(Xp(ip).GE.Box/2.)Xp(ip)=Xp(ip) -Box c IF(Yp(ip).GE.Box/2.)Yp(ip)=Yp(ip) -Box c IF(Zp(ip).GE.Box/2.)Zp(ip)=Zp(ip) -Box xmax =MAX(xmax,Xp(ip)) xmin =MIN(xmin,Xp(ip)) ymax =MAX(ymax,Yp(ip)) ymin =MIN(ymin,Yp(ip)) zmax =MAX(zmax,Zp(ip)) zmin =MIN(zmin,Zp(ip)) XL(ip)=(Yp(ip)-Ycc)*VZp(ip) & -(Zp(ip)-Zcc)*VYp(ip) YL(ip)=(Zp(ip)-Zcc)*VXp(ip) & -(Xp(ip)-Xcc)*VZp(ip) ZL(ip)=(Xp(ip)-Xcc)*VYp(ip) & -(Yp(ip)-Ycc)*VXp(ip) If(Jp(ip).le.Ndisk)Then Lx =Lx +XL(ip) Ly =Ly +YL(ip) Lz =Lz +ZL(ip) dd = (Xp(ip)-Xcc)**2 + & (Yp(ip)-Ycc)**2 + (Zp(ip)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd EndIf V2 = V2 + VXp(ip)**2 +VYp(ip)**2+VZp(ip)**2 EndDo Lx =Lx/Idisk Ly =Ly/Idisk Lz =Lz/Idisk V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/Idisk write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/Idisk), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*Idisk 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.3,/ & 25x,' Disk Ang.Momentum: L/ =',g11.4,' km/s',/ & 25x,' Disk Ang.Momentum: L =',g11.4,' kpc*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.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 ip=1,N Xnew = ex1*(Xp(ip)-Xcc) +ey1*(Yp(ip)-Ycc) +ez1*(Zp(ip)-Zcc)+Xcc Ynew = ex2*(Xp(ip)-Xcc) +ey2*(Yp(ip)-Ycc) +ez2*(Zp(ip)-Zcc)+Ycc Znew = ex3*(Xp(ip)-Xcc) +ey3*(Yp(ip)-Ycc) +ez3*(Zp(ip)-Zcc)+Zcc VXnew = ex1*Vxp(ip) +ey1*Vyp(ip) +ez1*Vzp(ip) VYnew = ex2*Vxp(ip) +ey2*Vyp(ip) +ez2*Vzp(ip) VZnew = ex3*Vxp(ip) +ey3*Vyp(ip) +ez3*Vzp(ip) Xp(ip) = Xnew Yp(ip) = Ynew Zp(ip) = Znew c Periodical Boundary Conditions c IF(Xp(ip).lt.0.)Xp(ip)=Xp(ip) +Box c IF(Yp(ip).lt.0.)Yp(ip)=Yp(ip) +Box c IF(Zp(ip).lt.0.)Zp(ip)=Zp(ip) +Box c IF(Xp(ip).GE.Box)Xp(ip)=Xp(ip) -Box c IF(Yp(ip).GE.Box)Yp(ip)=Yp(ip) -Box c IF(Zp(ip).GE.Box)Zp(ip)=Zp(ip) -Box Vxp(ip) = VXnew Vyp(ip) = VYnew Vzp(ip) = VZnew XL(ip)=(Yp(ip)-Ycc)*VZp(ip) & -(Zp(ip)-Zcc)*VYp(ip) YL(ip)=(Zp(ip)-Zcc)*VXp(ip) & -(Xp(ip)-Xcc)*VZp(ip) ZL(ip)=(Xp(ip)-Xcc)*VYp(ip) & -(Yp(ip)-Ycc)*VXp(ip) If(Jp(ip).gt.Ndisk)Then write(19,'(8g12.4)')Xp(ip)*1000.,Yp(ip)*1000., & Zp(ip)*1000., & XL(ip),YL(ip),ZL(ip),Fip(ip), & Fip(ip)+(Vxp(ip)**2+ & Vyp(ip)**2+Vzp(ip)**2)/2. Else write(20,'(8g12.4)')Xp(ip)*1000.,Yp(ip)*1000., & Zp(ip)*1000., & XL(ip),YL(ip), ZL(ip),Fip(ip), & Fip(ip)+(Vxp(ip)**2+ & Vyp(ip)**2+Vzp(ip)**2)/2. EndIf xmax =MAX(xmax,Xp(ip)) xmin =MIN(xmin,Xp(ip)) ymax =MAX(ymax,Yp(ip)) ymin =MIN(ymin,Yp(ip)) zmax =MAX(zmax,Zp(ip)) zmin =MIN(zmin,Zp(ip)) If(Jp(ip).le.Ndisk)Then Lx =Lx +XL(ip) Ly =Ly +YL(ip) Lz =Lz +ZL(ip) dd = (Xp(ip)-Xcc)**2 + & (Yp(ip)-Ycc)**2 + (Zp(ip)-Zcc)**2 dR = dR + sqrt(dd) dR2= dR2+ dd EndIf V2 = V2 + VXp(ip)**2 +VYp(ip)**2+VZp(ip)**2 EndDo Lx =Lx/Idisk Ly =Ly/Idisk Lz =Lz/Idisk V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/Idisk write (*,*) '--- After rotation of the coordinates' write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/Idisk),V2, & Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*Idisk 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 'PMgalaxy.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) EQUIVALENCE (Ndisk,extras(90)) 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 ip=1,N X = (Xp(ip)- Xcc) *1000. ! scale to kpc inits Y = (Yp(ip)- Ycc) *1000. Z = (Zp(ip)- Zcc) *1000. aLx =Y*VZp(ip) -Z*VYp(ip) aLy =Z*VXp(ip) -X*VZp(ip) aLz =X*VYp(ip) -Y*VXp(ip) dd = X**2 + Y**2 + Z**2 ds = sqrt(dd) ind = INT(aLz/hz+10000.) -10000 ind = Min(Max(ind,-Nang),Nang) If(ip.le.Ndisk)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(ip.le.N1)Then Lzd1(ind) =Lzd1(ind) +aLz Mzd1(ind) =Mzd1(ind) + 1 EndIf Else W = iWeight(ip) 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 If(ip.le.Ndisk+N1.and.ip.gt.Ndisk+N1-10000)Then Lzh1(ind) =Lzh1(ind) +aLz*W Mzh1(ind) =Mzh1(ind) + W EndIf EndIf EndDo write (*,100)N1,R1,N1,R1 write (17,100)N1,R1,N1,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) bL1= Lzh1(i)/Max(Mzh1(i),1) write(*,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i), & bL,Mzh(i),bL1,Mzh1(i),bL2,Mzh2(i) write(17,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i), & bL,Mzh(i),bL1,Mzh1(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------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Parameter (pi180 = 180./3.1415926535) Parameter (pi = 3.1415926535) EQUIVALENCE (Ndisk,extras(90)) drlog =0.033 d_lim = 0.5 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 ip=1,Ndisk X = (Xp(ip)- Xcc) *1000. ! scale to kpc inits Y = (Yp(ip)- Ycc) *1000. Z = (Zp(ip)- Zcc) *1000. rd = X**2 + Z**2 r_par = X*cosa +Z*sina c if( rd .lt. r_par**2) c & write (*,*) ' error in radius: ',sqrt(rd) ,r_par dd = sqrt(max(rd -r_par**2,1.e-10)) If(dd.le.d_lim)Then ! inside the strip ind = INT(LOG10(ABS(Y))/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Vradial = VXp(ip)*cosa +VZp(ip)*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 EndDo V1 =0 write(17,90) Alpha 90 format(' 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.2)Then write(17,100)rad_out,Raddisk, & Vradial/cosa,Vradial2/cosa,nn EndIf EndDo 100 format(2g10.3,2f6.1,i7) Return End C------------------------------------------ C SUBROUTINE Profile(N,Box,ScaleM,rb_min,rb_max) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Parameter (pi180 = 180./3.1415926535) Parameter (pi = 3.1415926535) EQUIVALENCE (Ndisk,extras(90)) drlog =0.033 dphi =180./Nbp Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. Do ip=1,N X = (Xp(ip)- Xcc) *1000. ! scale to kpc inits Y = (Yp(ip)- Ycc) *1000. Z = (Zp(ip)- Zcc) *1000. aLx =Y*VZp(ip) -Z*VYp(ip) aLy =Z*VXp(ip) -X*VZp(ip) aLz =X*VYp(ip) -Y*VXp(ip) 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(ip)*X+VYp(ip)*Y+VZp(ip)*Z) / & max(ds,1.d-5) Vvert = VZp(ip) c If(Abs(Z).le.1.5.and. Abs(VZp(ip)).le.50. c & .and.abs(Vradial).lt.50.)Then ind = INT(LOG10(ds)/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) If(ip.le.Ndisk)Then c If(ip/10*10.eq.ip)write(20,500)X,Y,Z c 500 format(3g13.5) Nradd(ind) = Nradd(ind) + 1 Radd(ind) = Radd(ind) + ds Vfi(ind) = Vfi(ind) + Vrotat Vfi2(ind) = Vfi2(ind) + Vrotat**2 Vrad(ind) = Vrad(ind) + Vradial Vrad2(ind) = Vrad2(ind) + Vradial**2 Vzz(ind) = Vzz(ind) + Vvert 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 phi =pi180*acos(X/rs) ! find angle in the plane of disk if(Y.le.0.)phi =360.-phi if(phi.ge.180.)phi =phi -180. ! limit phi to 0-180 deg inf = min(INT(phi/dphi) +1,Nbp) if(inf.gt.Nbp)write (*,*) ' ERRRRor in phi=',phi,inf c if(ip.le.200)write (*,*) ip,rs,phi,ind,inf 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 Else W = iWeight(ip) Nradh(ind) = Nradh(ind) + W Radh(ind) = Radh(ind) + ds*W Vfih(ind) = Vfih(ind) + Vrotat*W Vfih2(ind) = Vfih2(ind) + Vrotat**2*W Vradh(ind) = Vradh(ind) + Vradial*W Vradh2(ind) = Vradh2(ind) + Vradial**2*W EndIf c EndIF ! end abs(z) and abs(Vz) condition EndDo c --------------------- spherical shells V1 = 0. Din = 0. Dhin =0. write(17,90) 90 format(' Rout Rdisk VcTot ', & 'Vcdisk Dens_disk Vrot sig_rot ', & ' Vrad sigradial V_z sig_z Ndisk', & ' Halo: Rhalo ', & 'Vc_halo dens_halo Vrot sig_rot', & ' V_rad sig_rad Nhalo'/, & 4x,'kpc',6x,'kpc km/s km/s Msun/pc^3 km/s ') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 =4.1888*rad_out**3 nn = Nradd(i) Din = Din + nn n0 = max(nn,1) densD = nn/(V1-V0)*ScaleM*1.e-9 Raddisk = Radd(i)/n0 If(nn.eq.0)Raddisk=rad_out Vrotat = Vfi(i)/n0 Vrot2 = sqrt(max(Vfi2(i)/n0-Vrotat**2,1.e-10)) Vradial = Vrad(i)/n0 Vradial2=sqrt(max(Vrad2(i)/n0-Vradial**2,1.e-10)) Vvert = Vzz(i)/n0 Vvert2=sqrt(max(Vzz2(i)/n0-Vvert**2,1.e-10)) Vcirc = 2.08e-3*sqrt(max(Din*ScaleM/rad_out,1.d-10)) nh = Nradh(i) Dhin = Dhin + nh nh0 = max(nh,1) densH = nh/(V1-V0)*ScaleM*1.e-9 RadHalo = Radh(i)/nh0 If(nn.eq.0)RadHalo=rad_out VrotatH = Vfih(i)/nh0 Vrot2H = sqrt(max(Vfih2(i)/nh0-VrotatH**2,1.e-10)) VradialH = Vradh(i)/nh0 Vradial2H=sqrt(max(Vradh2(i)/nh0-VradialH**2,1.e-10)) VcircH= 2.08e-3*sqrt(max(Dhin*ScaleM/rad_out,1.d-10)) VcircT = sqrt(Vcirc**2+VcircH**2) If(nn.gt.2.or.nh.gt.2)Then write(17,100)rad_out,Raddisk,VcircT,Vcirc,densD, & Vrotat,Vrot2,Vradial,Vradial2,Vvert,Vvert2,nn, & RadHalo,VcircH,densH,VrotatH,Vrot2H, & VradialH,Vradial2H,nh EndIf EndDo 100 format(2g10.3,2f6.1,g11.3,6f7.1,i7, & g10.3,f7.1,g11.3,4f7.1,i8) c --------------------- cylindrical shells ind0 =0 ind1 =0 nn =0 Do i=-Nrad,Nrad ! find the angle of the bar rad_out = 10.**((i+1)*drlog) If(rad_out.gt.rb_min .and. rad_out.lt.rb_max)Then im = 0 imn =1000000 do jj =1,Nbp if(Nradb(i,jj).gt.im)Then im =Nradb(i,jj) jmax =jj endif if(Nradb(i,jj).lt.imn)Then imn =Nradb(i,jj) jmin =jj endif enddo write (19,*) ' max/min is at bin=',jmax,jmin,' Radius=', & rad_out ind0 = ind0 + jmax ind1 = ind1 + jmin nn = nn +1 Endif Enddo jmax = float(ind0)/Float(nn)+0.5 jmin = float(ind1)/float(nn)+0.5 write (19,'(" Final Max/min of bar density is at bin=",2i4, & " angles=",2f8.1)') & jmax,jmin,(jmax-0.5)*dphi,(jmin-0.5)*dphi V0 = V1 V1 = pi*rad_out**2 nn = Nradp(i) n0 = max(nn,1) densD = nn/(V1-V0)*ScaleM*1.e-6 V1 = 0. ! print results ninside = 0 dp = dphi/180. write(19,95) 95 format(' Rout Rdisk SurfDens Np Nin Vrot ', & 'sig_rot ','sig_rad sig_z ', & 'SurfBar Vrot_bar SurfoffBar Vrot_off_bar', & /4x,'kpc',6x,'kpc',7x,'Msun/pc^2',5x,' km/s km/s ') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 = pi*rad_out**2 nn = Nradp(i) ninside = ninside + nn n0 = max(nn,1) densD = nn/(V1-V0)*ScaleM*1.e-6 Raddisk = Radp(i)/n0 If(nn.eq.0)Raddisk=rad_out Vrotat = Vfip(i)/n0 Vrot2 = sqrt(max(Vfi2p(i)/n0-Vrotat**2,1.e-10)) Vradial = Vradp(i)/n0 Vradial2=sqrt(max(Vrad2p(i)/n0-Vradial**2,1.e-10)) Vvert = Vzzp(i)/n0 Vvert2=sqrt(max(Vzz2p(i)/n0-Vvert**2,1.e-10)) Vrotat0 = Vfib(i,jmax)/max(Nradb(i,jmax),1) Vrotat1 = Vfib(i,jmin)/max(Nradb(i,jmin),1) If(nn.gt.1)Then write(19,120)rad_out,Raddisk,densD,nn,ninside, & Vrotat,Vrot2,Vradial2,Vvert2, & ScaleM*1.e-6*Nradb(i,jmax)/dp/(V1-V0),Vrotat0, & ScaleM*1.e-6*Nradb(i,jmin)/dp/(V1-V0),Vrotat1 EndIf EndDo 120 format(2g11.3,g10.3,i6,i9,4f6.1,2(2x,g10.3,f8.1)) Return End C------------------------------------------ C Read particles SUBROUTINE ReadPnts(N,ScaleC,ScaleV,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' EQUIVALENCE (Ndisk,extras(90)) Real*8 xx,yy,zz,vvx,vvy,vvz Box =ScaleC*NGRID Radius2 =Radius**2 write (*,*) ' Number of disk particles=',Ndisk N =0 ! Number of particles 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 Icount_h =0 C Loop over particles 10 read(18,end=20, err=20)xx,yy,zz,vvx,vvy,vvz,gx,gy,gz, & phi,nn,Lev N = N + 1 If(nn.gt.Ndisk)Icount_h = Icount_h + 1 X =ScaleC* (xx-1.) Y =ScaleC* (yy-1.) Z =ScaleC* (zz-1.) c write (*,'(2i9,10g12.4)') N,nn, c & X,Y,Z,xx,yy,zz,phi Vxs=ScaleV* vvx Vys=ScaleV* vvy Vzs=ScaleV* vvz vmax =MAX(vmax,Vxs,Vys,Vzs) vmin =MIN(vmin,Vxs,Vys,Vzs) 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 Fip(N) = (ScaleV)**2*phi Jp(N) = nn 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)) GoTO 10 20 write (*,*) ' Number of particles selected:=',N, + ' halo=',Icount_h,' disk=',N-Icount_h write (*,*) ' Limits =',xmin,xmax write (*,*) ' =',ymin,ymax write (*,*) ' =',zmin,zmax Ihalo =Icount_h Idisk =N-Ihalo 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 'PMgalaxy.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 Nradh/Nrad2*0/,Vfih/Nrad2*0./,Vradh/Nrad2*0./, & Vfih2/Nrad2*0./,Vradh2/Nrad2*0./,Radh/Nrad2*0./ End C--------------------------------------------------- SUBROUTINE RDTAPE(File1,File2) C--------------------------------------------------- INCLUDE 'Statparam.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,Id,Norb,Rsnfw WRITE (16,100) HEADER, + AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, + EKIN, + NROWC,NGRID,NRECL,Om0,Oml0,hubble, + Ocurv,Id,Norb,Rsnfw 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, + 1x,' Nd= ',I6,1x,' Nt= ',I8,1x,' Rsnfw= ',f7.4) 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-------------------------------------------- INCLUDE 'Statparam.h' READ (20+Ifile,REC=IROW) RECDAT RETURN END