c-------------------------------------------------------------------- c Routines for analysis of simulations of barred galaxies c-------------------------------------------------------------------- C Get basic 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 Real*8 Xx2,Yy2,Zz2 EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) 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. 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(jp.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(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/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 (37,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 (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 write (18,150) AEXPN,time/1.e9,ISTEP,aL*Ndisk, & bar_angle*180./3.1415926, & bar_length,bar_amp, & 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,i5,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 'PMgalaxy.h' Real*8 MassTot EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. MassTot =0. Do jp=1,N MassTot = MassTot + iWeight(jp) EndDo MassTot = MassTot * ScaleM write(20,100) write(20,110) N,Ndisk,Nbulge,time,ISTEP, & ScaleM, ScaleM*Ndisk, MassTot write(20,115) c CALL GetDensity(Box) Do jp=1,Ndisk+Nbulge c Do jp=1,1e6 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) !,Dens(jp) c xx*VYp(jp)-yy*VXp(jp) c & INT(iWeight(jp)+0.1) EndDo 100 format(' TotalPart N: Disk Bulge', & 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', & T46,'Vx(kms) Vy Vz Density') 110 format(i10,i7,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 particles positions: c finds center of the system, tilt of the disk, c orientation of the bar c Code finds orientation of the bar and rotates the reference c frame so the bar is along x-axis 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 Real*8 Xx,Yy,Zz,XY,dd,Xc0,Yc0,Zc0 EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Ncenter = N/2. deg = 180./pi ! convert to degrees 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.550 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,f8.3,g11.3, & ' X(Mpc)=',3g12.6,/ & 25x,'shift(kpc) =',3g12.6,/ & 25x,' average Velocity(km/s)=',3g11.4) Xc0 = Xc Yc0 = Yc Zc0 = Zc Ncount = Ncount +1 Rselect = Rselect/1.05 R2 = Rselect**2 dn = nn/Rselect**3 If(Ncount.lt.1000.and.Rselect.gt.3.e-3.and.nn.gt.3e3)goto300 write (37,100) nn,Rselect*1000.,nn/Rselect**3,Xc,Yc,Zc, & 1000.*xshift,1000.*yshift,1000.*zshift, & Vcx,Vcy,Vcz C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (jp) 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 EndDo Call BarPosition(Xcc,Ycc,Zcc) write (*,*) ' bar_angle=',bar_angle*deg Xcos = cos(bar_angle) Ysin = sin(bar_angle) C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (jp,xold,yold,vxold,vyold) Do jp=1,N xold = Xp(jp)-Xcc yold = Yp(jp)-Ycc Xp(jp) =Xcc +xold*Xcos +yold*Ysin Yp(jp) =Ycc -xold*Ysin +yold*Xcos vxold = VXp(jp) vyold = VYp(jp) VXp(jp) = vxold*Xcos +vyold*Ysin VYp(jp) = -vxold*Ysin +vyold*Xcos c xx = (Xp(jp)-Xcc)*1000. c yy = (Yp(jp)-Ycc)*1000. c zz = (Zp(jp)-Zcc)*1000. c write (25,130) xx,yy,zz, c & VXp(jp) ,VYp(jp),VZp(jp) c 130 format(3g13.5,2x,4g12.4,i3) EndDo Do jp=1,N 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(jp.le.Ndisk)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 c Xcos = cos(bar_angle) c Ysin = sin(bar_angle) Lx =Lx/Ndisk Ly =Ly/Ndisk Lz =Lz/Ndisk V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/Ndisk write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/Ndisk), & V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*Ndisk 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. C$OMP PARALLEL DO DEFAULT(SHARED) C$OMP+PRIVATE (jp,Xnew,Ynew,Znew,VXnew,VYnew,VZnew) 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 EndDo Do jp=1,N 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(jp.le.Ndisk)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/Ndisk Ly =Ly/Ndisk Lz =Lz/Ndisk V2 =sqrt(V2/N) aL = sqrt(Lx**2 +Ly**2 +Lz**2) dR = dR/Ndisk write (*,*) '--- After rotation of the coordinates' write (*,110) xmin,xmax,ymin,ymax,zmin,zmax, & 1000.*dR,1000.*sqrt(dR2/Ndisk),V2, & Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*Ndisk 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 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(jp.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(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 If(jp.le.Ndisk+N1.and.jp.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 (37,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(37,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 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 'PMgalaxy.h' Parameter (pi180 = 180./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 accumulate statistics for bulge : C SUBROUTINE ProfileBulge(N,Box,ScaleM,drlog) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. If(Nbulge.eq.0)return Do jp=Ndisk+1,Ndisk+Nbulge 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) 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) Nradg(ind) = Nradg(ind) + 1 Radg(ind) = Radg(ind) + ds Vfig2(ind) = Vfig2(ind) + Vtan2 Vradg(ind) = Vradg(ind) + Vradial Vradg2(ind) = Vradg2(ind) + Vradial**2 ind = INT(LOG10(rs)/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Nradu(ind) = Nradu(ind) + 1 Vfig(ind) = Vfig(ind) + Vrotat EndDo V1 = 0. Din = 0. write(37,190) 190 format(1x,20('-'),' Bulge profile: Spherical shells',/ & ' Rout Rbulge Vcbulge', & ' Density ',' Vrad',5x,'sigRadial', & 4x,'V_tan',6x,'Vrot',6x,'Beta',11x,'Nbin Ninside') Do i=-Nrad,Nrad rad_out = 10.**((i+1)*drlog) V0 = V1 V1 = 4.1888*rad_out**3 wn = Nradg(i) nn =wn Din = Din + wn w0 = max(wn,1.e-9) densD = wn/(V1-V0)*ScaleM*1.e-9 Radbulge = Radg(i)/w0 If(wn.eq.0)Radbulge=rad_out Vtan = sqrt(Vfig2(i)/w0) Vradial = Vradg(i)/w0 Vrotat = Vfig(i)/max(1,Nradu(i)) Vradial2 =sqrt(max(Vradg2(i)/w0-Vradial**2,1.e-10)) Vtot = sqrt((Vradg2(i)+Vfig2(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,120)rad_out,Radbulge,Vcirc,densD, & Vradial,Vradial2,Vtan,Vrotat,beta,nn,INT(Din) EndIf EndDo 120 format(2g10.3,f6.2, & g11.3,2x,5g11.4,2x,i7,i9) Return End c-------------------------------------------------------------------- C accumulate statistics for disk: C SUBROUTINE ProfileDisk(N,Box,ScaleM,drlog) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. Do jp=1,Ndisk 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) c Vradial = (VXp(jp)*X+VYp(jp)*Y+VZp(jp)*Z) / c & max(ds,1.e-5) Vradial = (VXp(jp)*X+VYp(jp)*Y) / & max(rs,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) IF(ABS(z).lt.dZ_slice.or.dZ_slice.lt.1.e-10)Then c IF(ABS(Vvert).lt.20..and.ABS(Vradial).lt.20.)Then 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 c EndIf EndIf EndDo Return End c-------------------------------------------------------------------- C random velocities for disk: C SUBROUTINE RandV(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Logical inside real*8 vx2,vy2,vz2 EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) scaleX = 1./0.4995 ! rescale all variables ScaleMass = 1./0.03047 ! Mass_new=Mass_old/ScaleMass c scaleX = 1. c ScaleMass = 1. scaleT = sqrt(scaleX**3/ScaleMass) ! t_new= t_old/scaleT; Omega_new=Om_old*scaleT ScaleV = ScaleX/scaleT ! V_new= V_old/ScaleV write (*,'(" scaling factors: radius=",f8.4)')scaleX write (*,'(" scaling factors: time =",f8.4)')scaleT write (*,'(" scaling factors: mass =",f8.4)')scaleMass write (*,'(" scaling factors: veloc =",f8.4)')scaleV write (20,'(" scaling factors: radius=",f8.4)')scaleX write (20,'(" scaling factors: time =",f8.4)')scaleT write (20,'(" scaling factors: mass =",f8.4)')scaleMass write (20,'(" scaling factors: veloc =",f8.4)')scaleV Xcc = Box/2. Ycc = Box/2. Zcc = Box/2. drlog =0.05 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. Vzzp(i) = 0. Vzz2p(i) = 0. Vfip(i) = 0. Vfi2p(i) = 0. EndDo Do jp=1,Ndisk ! rescale coordinates and Vs Xp(jp) = (Xp(jp)- Xcc) *1000./scaleX ! scale to kpc inits Yp(jp) = (Yp(jp)- Ycc) *1000./scaleX Zp(jp) = (Zp(jp)- Zcc) *1000./scaleX VXp(jp)= VXp(jp)/scaleV VYp(jp)= VYp(jp)/scaleV VZp(jp)= VZp(jp)/scaleV EndDo Do jp=1,Ndisk ! find velocity deviations X = Xp(jp) Y = Yp(jp) Z = Zp(jp) Vxx = VXp(jp) Vyy = VYp(jp) Wzz = VZp(jp) if(abs(Z).lt.0.2.and.abs(Wzz).lt.10.)Then dc = X**2 + Y**2 dd = dc + Z**2 ds = sqrt(dd) dc = sqrt(dc) dr2 = (0.07*ds**0.7)**2 c goto 25 niter = 0 30 nn = 0 vx2 = 0. vy2 = 0. vz2 = 0. Do kp=1,Ndisk ! find average velocity of neigbours Xk = Xp(kp) Yk = Yp(kp) Zk = Zp(kp) dd = (Xk-X)**2 + (Yk-Y)**2 + (Zk-Z)**2 If(dd.lt.dr2)Then nn = nn + 1 vx2 = VXp(kp) + vx2 vy2 = VYp(kp) + vy2 vz2 = VZp(kp) + vz2 EndIf EndDo if(nn.lt.50)Then dr2 =dr2*1.5 niter =niter +1 goto 30 EndIf if(nn.gt.70)Then dr2 =dr2/1.4 niter =niter +1 goto 30 EndIf if(mod(jp,1001).eq.0) & write (*,180) jp,ds,niter,nn,sqrt(dr2),sqrt(dr2)/(0.1*ds**0.7) 180 format(i8,f8.3,i4,i7,3f8.3) 25 vx2 = vx2/nn vy2 = vy2/nn vz2 = vz2/nn dvv = sqrt((Vxx-vx2)**2+(Vyy-vy2)**2+(Wzz-vz2)**2) if(dvv.lt.15.)Then dvx = Vxx !-vx2 dvy = Vyy !-vy2 dvz = Wzz !-vz2 c dvx = vx2 c dvy = vy2 c dvz = vz2 Vradial = (dvx*X+dvy*Y)/max(dc,1.e-5) dvv = sqrt(dvx**2+dvy**2+dvz**2) dvp = (dvx*Y-dvy*X)/max(dc,1.e-5) ind = INT(LOG10(ABS(dc))/drlog+1000.) -1000 ind = min(max(ind,-Nrad),Nrad) Nradp(ind) = Nradp(ind) + 1 Radp(ind) = Radp(ind) + abs(dc) Vradp(ind) = Vradp(ind) + Vradial Vrad2p(ind) = Vrad2p(ind) + Vradial**2 Vfip(ind) = Vfip(ind) + dvp Vfi2p(ind) = Vfi2p(ind) + dvp**2 Vzzp(ind) = Vzzp(ind) + dvz Vzz2p(ind) = Vzz2p(ind) + dvz**2 endif endif EndDo iFlag = 0 write (20,190) Do ind =-Nrad,Nrad nn = max(Nradp(ind),1) Radp(ind) = Radp(ind)/nn Vradp(ind) = Vradp(ind)/nn Vrad2p(ind)= sqrt(Vrad2p(ind)/nn -Vradp(ind)**2) Vzzp(ind) = Vzzp(ind)/nn Vzz2p(ind) = sqrt(Vzz2p(ind)/nn -Vzzp(ind)**2) Vfip(ind) = Vfip(ind)/nn Vfi2p(ind) = sqrt(Vfi2p(ind)/nn-Vfip(ind)**2) If(nn.gt.10.or.iFlag.eq.1)Then iFlag =1 write (20,200) nn,Radp(ind),Vradp(ind),Vrad2p(ind), & Vzzp(ind),Vzz2p(ind),Vfip(ind),Vfi2p(ind) if(Nradp(ind).eq.0.and.Nradp(ind-1).eq.0)GoTo 210 EndIf EndDo 200 format(i5,f8.4,T15,2f8.2,T35,2f8.2,T55,2f8.2) 190 format(' N ','R(kpc)',T15,' Vrad Sig_rad', & T35,' Vz Sig_z', & T55,' Vphi Sig_phi') 210 Return End c-------------------------------------------------------------------- C accumulate statistics for disk: C SUBROUTINE PV(N,Box,ScaleM) C------------------------------------------ INCLUDE 'PMparameters.h' INCLUDE 'PMgalaxy.h' Logical inside EQUIVALENCE (Ndisk,extras(90)), & (Nbulge,extras(97)),(aMassDisk,extras(93)) r0 = 8. ! distance from the center dr0 = 0.25 ! radius of sphere to find velocity of observer angle = -20. ! position angle of observer devV = 20. ! condition for 'cold' particles: |V|