c Program to calculate theoretical King's (1966) denity profile. c 2: Produces an N-body snapshot of an isotropic King model. Based on eqs. (4-130) - (4-132) c from Galactic Dynamics by Binney & Tremaine (1994, p. 232), and eq. (15) from c ApJS, 131, 39 (2000). Output format: NEMO ASCII (http://bima.astro.umd.edu/nemo). c 2.01: Cleaned version. c Maximum number of N-body particles: #define NN 1000000 c Maximum number of points for radial density profile (for ODE integrator): #define KMAXX 1000000 c Number of ODEs (fixed): #define Nmax 2 program King2 implicit NONE real*8 alpha, Eps,t0,t2,A(2),dt1, dxsav, xp(KMAXX), # yp(KMAXX,2),rho,C,mu,con,beta, # r_tid,dy2dxp(KMAXX),P(KMAXX),P1,I0, # erfc1,PB, PEx, m1, R, V(NN), Rx,Ry,Rz, # Vx,Vy,Vz, PE,E, zriddr, cos_t,sin_t,phi,sum,Psi, # Interpolation,Gs,Mt,sgm,r0,Rkpc,rho1,RR,sgm1 real*4 ran2 integer nok,nbad,kmax, kount,k,kt, N, TT(3), Seed, i logical Error,ex external derivs,rkqs,PE common /path/ kmax, kount, dxsav, xp, yp common /com_PE/ Psi, I0, PB, PEx include 'const.f' c Initializing the seed for random generators: inquire (file='Seed.tmp', exist=ex) if (ex) then open (2,file='Seed.tmp',status='old') read (2,*) Seed Seed = -Seed close (2) else call Itime (TT) Seed = -(TT(1)+TT(2)+TT(3)) end if c G in local units (which are assumed to be M_sun, kpc, and km/s): Gs = G *M_Sun/1d10/kpc c Accuracy for ODE solver: Eps = 1.d-9 c Maximum dimensionless radius for ODE integration: t0 = 10d0 c Initial (very close to 0) dimensionless radius for ODE integration: dt1 = 1.d-15 kmax = KMAXX write(*,*) ' Number of particles N:' read (*,*) N write (*,*) N write(*,*) ' Parameter alpha=psi/sgm^2:' read (*,*) alpha write (*,*) alpha write(*,*) ' Core radius r0 (kpc):' read (*,*) r0 write (*,*) r0 write(*,*) ' Velocity dispersion parameter sgm (km/s):' read (*,*) sgm write (*,*) sgm if (N .gt. NN) then write (*,*) 'Please recompile with NN >=',N stop end if t2 = t0 1 continue c Poisson equation solving to find radial profile of potential energy: A(1) = alpha A(2) = 0.d0 dxsav = t2 / KMAXX call odeint (A,2,dt1,t2,Eps,dt1,0.0d0,nok,nbad, # derivs,rkqs) if (A(1) .gt. 0.d0) then c In the unlikely event of r_tidal > t2: t2 = t2 * 3.d0 write (*,*) t2 go to 1 end if c Writing file with the radial profile. open (1,file='king.dat') write (1,'(A1,3(A13))') '#','r ','rho ','sgm ' C = dsqrt(dexp(alpha)*derf(dsqrt(alpha)) - # 2.d0*dsqrt(alpha/pi)*(1.d0+2.d0/3.d0*alpha)) mu = 0.d0 beta = 0.d0 r_tid = 0.d0 con = 1.d10 do k=1,kount rho = dexp(yp(k,1))*derf(dsqrt(yp(k,1))) - # 2.d0*dsqrt(yp(k,1)/pi)*(1.d0+2.d0/3.d0*yp(k,1)) dy2dxp(k) = -xp(k)**2 * rho if (rho .ge.0.d0) then kt = k R = xp(k) / 3.d0 * C con = R rho1 = rho * (3.d0*sgm/C/r0)**2 / 4.d0/Pi/Gs RR = R * r0 c 1D velocity dispersion (in km/s): if (yp(k,1) .ge. 3.d-4) then sgm1 = sgm * dsqrt(1.d0 - (2.d0*yp(k,1))**2.5d0 / # (15.d0*dexp(yp(k,1))*dsqrt(Pi/2.d0)*derf(dsqrt(yp(k,1))) # - 5.d0*dsqrt(2.d0*yp(k,1))*(3.d0+2.d0*yp(k,1)))) else c Asymptotic behavior for small yp(k,1): sgm1 = sgm * dsqrt(2.d0/7.d0*yp(k,1)) end if c Outputting the theoretical radial profiles: radius (kpc), density (M_Sun/kpc^3), c and 1D velocity dispersion (km/s): write (1,'(1PE12.6,2(1PE13.6))') RR,rho1,sgm1 end if end do close(1) write (*,*) write (*,'(A,F8.5)') 'Concentration c =',dlog10(con) write (*,'(A,1PE13.6,A)') 'Tidal radius r_t =',xp(kt)/3.d0*C*r0, # ' kpc' c N-body realization of the King model. sum = 0.d0 c Calculating P(k) - integral probability function for radial particle distribution: do k=2,kt sum = sum + (dy2dxp(k-1)+dy2dxp(k)) * (xp(k)-xp(k-1)) P(k) = sum end do P(1) = 0.d0 c Total mass in M_Sun: Mt = -C*r0*sgm**2/6.d0/Gs * P(kt) write (*,'(A,1PE13.6,A)') 'Total mass M =',Mt,' M_Sun' c Normalizing P(k) to 1 at r_tid: do k=2,kt P(k) = P(k)/P(kt) end do c Generating coordinates and isotropic velocities of particles c Mass of one particle (M_Sun): m1 = Mt / dble(N) c Writing output in the ascii format which can be transformed to snapshot NEMO c (http://bima.astro.umd.edu/nemo/) format with NEMO program "atos": c atos king_snap.dat king_snap.s open (1, file='king_snap.dat') c Number of bodies: write (1,*) N c Number of dimensions: write (1,'(I1)') 3 c Initial moment of time: write (1,*) 0.0 c Particle masses: do i=1,N write (1,'(1PE12.6)') m1 end do do i=1,N 2 continue c Generating R (dimensionless distance from the center): P1 = dble(ran2(Seed)) c 3rd order interpolation to get R R = Interpolation (P,xp,kt,P1,3,error) if (Error) go to 2 c 3rd order interpolation to get psi(R) - relative potential Psi = Interpolation (xp,yp,kt,R,3,error) c Generating E (dimensionless relative energy): if (Psi .gt. 1.d0) then erfc1 = derfc(dsqrt(Psi)) else erfc1 = 1.d0 - derf(dsqrt(Psi)) end if PB = 0.5d0*dsqrt(Pi)*dexp(Psi) I0 = dsqrt(Psi)*(2.d0/3.d0*Psi+1.d0) + PB*erfc1 PEx = dble(ran2(Seed)) c Solving non-linear equation to get E: E = zriddr(PE,0.d0,Psi,1.d-9) c Velocity vector module (km/s): V(i) = sgm * dsqrt(2.d0*(Psi-E)) c Three components of the isotropic random vectors R: cos_t = dble(ran2(Seed))*2.d0 - 1.d0 sin_t = dsqrt(1.d0-cos_t**2) phi = dble(ran2(Seed))*2.d0*Pi Rkpc = R/3.d0 * C * r0 Rx = Rkpc * sin_t * dcos(phi) Ry = Rkpc * sin_t * dsin(phi) Rz = Rkpc * cos_t write (1,'(1PE13.6,2(1PE14.6))') Rx,Ry,Rz end do do i=1,N cos_t = dble(ran2(Seed))*2.d0 - 1.d0 sin_t = dsqrt(1.d0-cos_t**2) phi = dble(ran2(Seed))*2.d0*Pi Vx = V(i) * sin_t * dcos(phi) Vy = V(i) * sin_t * dsin(phi) Vz = V(i) * cos_t write (1,'(1PE13.6,2(1PE14.6))') Vx,Vy,Vz end do close (1) open (2,file='Seed.tmp') write (2,*) Seed close (2) stop end c===================================================================================== function Interpolation (x,y,N,x1,Order,error) c Interpolation = y(x1) with a polynom of Order's order. implicit NONE integer N,Order,j,m,k real*8 Interpolation, x(N),y(N),dx,x1 logical error error = .false. call locate(x,N,x1,j) if (j.eq.0 .or. j.eq.N+1) then write (*,*) 'Out of range during interpolation!' error = .true. end if c Cubic (m-1) polynom interpolation: m = Order + 1 k = min(max(j-(m-1)/2,1),N+1-m) call polint(x(k),y(k),m,x1,Interpolation,dx) return end c===================================================================================== function PE (E) c Integral probability function for relative energy minus PEx. implicit NONE real*8 PE, E, I0, PB, Psi, PsiE, erfc1, PEx common /com_PE/ Psi, I0, PB, PEx PsiE = dsqrt(Psi - E) if (PsiE .gt. 1.d0) then erfc1 = derfc(PsiE) else erfc1 = 1.d0 - derf(PsiE) end if PE = (PsiE * (2.d0/3.d0*(Psi-E)+dexp(E)) + PB*erfc1 - I0) / # (PB-I0) - PEx return end c===================================================================================== subroutine derivs (x,y,dydx) c Right side of ODEs. implicit NONE real*8 x, y(2), dydx(2), x2, sqy c real erf include 'const.f' x2 = x**2 sqy = dsqrt(y(1)) dydx(1) = y(2)/x2 dydx(2) = -x2*(dexp(y(1))*derf(sqy) - # 2.d0*sqy/dsqrt(pi)*(1.d0+2.d0/3.d0*y(1))) return end c-------------------------------------------------------------------------------- SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs, *rkqs) INTEGER nbad,nok,nvar,MAXSTP REAL*8 eps,h1,hmin,x1,x2,ystart(nvar),TINY EXTERNAL derivs,rkqs PARAMETER (MAXSTP=10000,TINY=1.d-30) INTEGER i,kmax,kount,nstp REAL*8 dxsav,h,hdid,hnext,x,xsav,dydx(Nmax),xp(KMAXX),y(Nmax), *yp(KMAXX,Nmax),yscal(Nmax) COMMON /path/ kmax,kount,dxsav,xp,yp x=x1 h=sign(h1,x2-x1) nok=0 nbad=0 kount=0 do 11 i=1,nvar y(i)=ystart(i) 11 continue if (kmax.gt.0) xsav=x-2.*dxsav do 16 nstp=1,MAXSTP call derivs(x,y,dydx) do 12 i=1,nvar yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY 12 continue if(kmax.gt.0)then if(abs(x-xsav).gt.abs(dxsav)) then if(kount.lt.kmax-1)then kount=kount+1 xp(kount)=x do 13 i=1,nvar yp(kount,i)=y(i) 13 continue xsav=x endif endif endif if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs) if(hdid.eq.h)then nok=nok+1 else nbad=nbad+1 endif if((x-x2)*(x2-x1).ge.0.)then do 14 i=1,nvar ystart(i)=y(i) 14 continue if(kmax.ne.0)then kount=kount+1 xp(kount)=x do 15 i=1,nvar yp(kount,i)=y(i) 15 continue endif return endif if(abs(hnext).lt.hmin) pause *'stepsize smaller than minimum in odeint' h=hnext 16 continue pause 'too many steps in odeint' return END c------------------------------------------------------------------- SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) INTEGER n REAL*8 eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n) EXTERNAL derivs CU USES derivs,rkck INTEGER i REAL*8 errmax,h,htemp,xnew,yerr(Nmax),ytemp(Nmax),SAFETY,PGROW, *PSHRNK,ERRCON PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4) h=htry 1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs) errmax=0. do 11 i=1,n errmax=max(errmax,abs(yerr(i)/yscal(i))) 11 continue errmax=errmax/eps if(errmax.gt.1.)then htemp=SAFETY*h*(errmax**PSHRNK) h=sign(max(abs(htemp),0.1*abs(h)),h) xnew=x+h if(xnew.eq.x)pause 'stepsize underflow in rkqs' goto 1 else if(errmax.gt.ERRCON)then hnext=SAFETY*h*(errmax**PGROW) else hnext=5.*h endif hdid=h x=x+h do 12 i=1,n y(i)=ytemp(i) 12 continue return endif END c------------------------------------------------------------------- SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs) INTEGER n REAL*8 h,x,dydx(n),y(n),yerr(n),yout(n) EXTERNAL derivs CU USES derivs INTEGER i REAL*8 ak2(Nmax),ak3(Nmax),ak4(Nmax),ak5(Nmax),ak6(Nmax), *ytemp(Nmax),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51,B52,B53, *B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6 PARAMETER (A2=.2,A3=.3,A4=.6,A5=1.,A6=.875,B21=.2,B31=3./40., *B32=9./40.,B41=.3,B42=-.9,B43=1.2,B51=-11./54.,B52=2.5, *B53=-70./27.,B54=35./27.,B61=1631./55296.,B62=175./512., *B63=575./13824.,B64=44275./110592.,B65=253./4096.,C1=37./378., *C3=250./621.,C4=125./594.,C6=512./1771.,DC1=C1-2825./27648., *DC3=C3-18575./48384.,DC4=C4-13525./55296.,DC5=-277./14336., *DC6=C6-.25) do 11 i=1,n ytemp(i)=y(i)+B21*h*dydx(i) 11 continue call derivs(x+A2*h,ytemp,ak2) do 12 i=1,n ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i)) 12 continue call derivs(x+A3*h,ytemp,ak3) do 13 i=1,n ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i)) 13 continue call derivs(x+A4*h,ytemp,ak4) do 14 i=1,n ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+B54*ak4(i)) 14 continue call derivs(x+A5*h,ytemp,ak5) do 15 i=1,n ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+B64*ak4(i)+ *B65*ak5(i)) 15 continue call derivs(x+A6*h,ytemp,ak6) do 16 i=1,n yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i)) 16 continue do 17 i=1,n yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)+DC6* *ak6(i)) 17 continue return END c------------------------------------------------------------------- FUNCTION ran2(idum) INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END c------------------------------------------------------------------- SUBROUTINE locate(xx,n,x,j) INTEGER j,n REAL*8 x,xx(n) INTEGER jl,jm,ju jl=0 ju=n+1 10 if(ju-jl.gt.1)then jm=(ju+jl)/2 if((xx(n).ge.xx(1)).eqv.(x.ge.xx(jm)))then jl=jm else ju=jm endif goto 10 endif if(x.eq.xx(1))then j=1 else if(x.eq.xx(n))then j=n-1 else j=jl endif return END c----------------------------------------------------------------------------- SUBROUTINE polint(xa,ya,n,x,y,dy) implicit NONE INTEGER n,NMAX REAL*8 dy,x,y,xa(n),ya(n) PARAMETER (NMAX=10) INTEGER i,m,ns REAL*8 den,dif,dift,ho,hp,w,c(NMAX),d(NMAX) ns=1 dif=abs(x-xa(1)) do 11 i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(i) d(i)=ya(i) 11 continue y=ya(ns) ns=ns-1 do 13 m=1,n-1 do 12 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.0.)pause 'failure in polint' den=w/den d(i)=hp*den c(i)=ho*den 12 continue if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 endif y=y+dy 13 continue return END c------------------------------------------------------------------- FUNCTION zriddr(func,x1,x2,xacc) INTEGER MAXIT REAL*8 zriddr,x1,x2,xacc,func,UNUSED PARAMETER (MAXIT=60,UNUSED=-1.11E30) EXTERNAL func CU USES func INTEGER j REAL*8 fh,fl,fm,fnew,s,xh,xl,xm,xnew fl=func(x1) fh=func(x2) if((fl.gt.0..and.fh.lt.0.).or.(fl.lt.0..and.fh.gt.0.))then xl=x1 xh=x2 zriddr=UNUSED do 11 j=1,MAXIT xm=0.5*(xl+xh) fm=func(xm) s=sqrt(fm**2-fl*fh) if(s.eq.0.)return xnew=xm+(xm-xl)*(sign(1.,fl-fh)*fm/s) if (abs(xnew-zriddr).le.xacc) return zriddr=xnew fnew=func(zriddr) if (fnew.eq.0.) return if(sign(fm,fnew).ne.fm) then xl=xm fl=fm xh=zriddr fh=fnew else if(sign(fl,fnew).ne.fl) then xh=zriddr fh=fnew else if(sign(fh,fnew).ne.fh) then xl=zriddr fl=fnew else pause 'never get here in zriddr' endif if(abs(xh-xl).le.xacc) return 11 continue pause 'zriddr exceed maximum iterations' else if (fl.eq.0.) then zriddr=x1 else if (fh.eq.0.) then zriddr=x2 else pause 'root must be bracketed in zriddr' endif return END c-------------------------------------------------------------------