Nbody6
 All Files Functions Variables
xtrnl0.f
Go to the documentation of this file.
1  SUBROUTINE xtrnl0
2 *
3 *
4 * External force initialization.
5 * ------------------------------
6 *
7  include 'common6.h'
8  common/galaxy/ gmg,rg(3),vg(3),fg(3),fgd(3),tg,
9  & omega,disk,a,b,v02,rl2,gmb,ar,gam,zdum(7)
10 *
11 *
12 * Check option for cluster in circular galactic orbit.
13  IF (kz(14).NE.1) go to 20
14 *
15 * Specify Oort's constants (units of km/sec/kpc).
16  a = 14.4
17  b = -12.0
18 * Adopt local density from Gilmore & Kuijken (solar mass/pc**3).
19  rho = 0.11
20 * Convert rotation constants to units of cm/sec/pc.
21  a = 100.0*a
22  b = 100.0*b
23 *
24 * Specify the tidal term in star cluster units (solar mass & pc).
25  tidal(1) = 4.0*a*(a - b)*(pc/gm)
26 *
27 * Initialize the Y-component to zero.
28  tidal(2) = 0.0
29 *
30 * Specify the vertical force gradient.
31  tidal(3) = -(2.0*twopi*rho + 2.0*(a - b)*(a + b)*(pc/gm))
32 *
33  fac = 1.0e-10/(pc/gm)
34  WRITE (6,5) zmbar*zmass, fac*tidal(1), fac*tidal(3), pc/gm
35  5 FORMAT (/,12x,'TOTAL MASS =',f8.1,' TIDAL(1&3) =',1p,2e10.2,
36  & ' PC/GM =',e10.2)
37 *
38 * Adopt twice the angular velocity for Coriolis terms.
39  tidal(4) = 2.0*(a - b)*sqrt(pc/gm)
40 *
41 * Define time scale in seconds and velocity scale in km/sec.
42  tscale = sqrt(pc/gm)*pc
43 *
44 * Convert time scale from units of seconds to million years.
45  tscale = tscale/(3.15d+07*1.0d+06)
46 *
47 * Scale to working units of RBAR in pc & ZMBAR in solar masses.
48  DO 10 k = 1,3
49  tidal(k) = tidal(k)*rbar**3/zmbar
50  10 CONTINUE
51  tidal(4) = tidal(4)*sqrt(rbar**3/zmbar)
52  tscale = tscale*sqrt(rbar**3/(zmass*zmbar))
53 *
54 * Perform an energy scaling iteration to include ETIDE (once is OK).
55  e0 = -0.25
56  CALL xtrnlv(1,n)
57 * Form total energy based on current values and ETIDE (cf. SCALE).
58  etot = zkin - pot + etide
59  sx = e0/etot
60 * Exclude re-scaling for positive energy or binary Plummer cluster.
61  IF (etot.GE.0.0) sx = 1.0
62  IF (kz(5).EQ.2) sx = 1.0
63 *
64 * Scale coordinates & velocities to yield ETOT = -0.25.
65  DO 15 i = 1,n
66  DO 12 k = 1,3
67  x(k,i) = x(k,i)/sx
68  xdot(k,i) = xdot(k,i)*sqrt(sx)
69  12 CONTINUE
70  15 CONTINUE
71 *
72 * Consider alternatives: circular point-mass orbit or 3D galaxy model.
73  20 zmtot = zmass*zmbar
74  IF (kz(14).EQ.2) THEN
75 *
76 * Read galaxy mass and central distance (solar units and kpc).
77  READ (5,*) gmg, rg0
78 *
79 * Set circular velocity in km/sec and angular velocity in cgs units.
80  vg0 = 1.0d-05*sqrt(gmg/(1000.0*rg0))*sqrt(gm/pc)
81  omega = 100.0*vg0/rg0
82 *
83 * Obtain King tidal radius in pc (eq. (9) of Fukushige & Heggie, 1995).
84  rt = (zmtot/(3.0*gmg))**0.3333*(1000.0*rg0)
85 *
86  IF (rtide.GT.0.0) THEN
87 * Determine RBAR (N-body units) from RT (pc) and King model (see SCALE).
88  IF(kz(22).EQ.2) rbar = rt/rtide
89  ELSE
90  rtide = rt/rbar
91  END IF
92 *
93 * Convert from cgs to N-body units.
94  omega = omega*sqrt(pc/gm)*sqrt(rbar**3/zmbar)
95 *
96 * Specify the galactic parameters for equations of motion.
97  tidal(1) = 3.0*omega**2
98  tidal(2) = 0.0d0
99  tidal(3) = -omega**2
100  tidal(4) = 2.0*omega
101  gmg = gmg/zmtot
102 *
103 * Check re-scaling units to current RBAR (i.e. TSTAR & VSTAR).
104  IF (kz(22).EQ.2) THEN
105  CALL units
106  END IF
107 *
108  WRITE (6,35) gmg, rg0, omega, rtide, rbar
109 *
110 * Perform an energy scaling iteration to include ETIDE (once is OK).
111  e0 = -0.25
112  CALL xtrnlv(1,n)
113 * Form total energy based on current values and ETIDE (cf. SCALE).
114  etot = zkin - pot + etide
115  sx = e0/etot
116  IF (etot.GE.0.0) sx = 1.0
117 *
118 * Scale coordinates & velocities to yield ETOT = -0.25.
119  DO 24 i = 1,n
120  DO 22 k = 1,3
121  x(k,i) = x(k,i)/sx
122  xdot(k,i) = xdot(k,i)*sqrt(sx)
123  22 CONTINUE
124  24 CONTINUE
125 *
126 * Treat the case of 3D orbit for point-mass, bulge, disk and/or halo.
127  ELSE IF (kz(14).EQ.3) THEN
128 *
129 * Read all parameters (NB! Do not confuse with Oort's constants A, B).
130  READ (5,*) gmg, disk, a, b, vcirc, rcirc, gmb, ar, gam
131 * Gamma/eta model: GMB = mass, AR = scale radius, GAM = exponent.
132  READ (5,*) (rg(k),k=1,3), (vg(k),k=1,3)
133 *
134 * Specify planar motion from SEMI & ECC for no disk & halo if VZ = 0.
135  IF (disk + vcirc + gmb.EQ.0.0.AND.vg(3).EQ.0.0d0) THEN
136  rap = rg(1)
137  ecc = rg(2)
138  semi = rap/(1.0 + ecc)
139  vg2 = gmg/(1000.0*semi)*(1.0 - ecc)/(1.0 + ecc)
140  DO 25 k = 1,3
141  rg(k) = 0.0
142  vg(k) = 0.0
143  25 CONTINUE
144 * Initialize 2D orbit with given eccentricity at apocentre.
145  rg(1) = rap
146  vg(2) = 1.0d-05*sqrt(vg2)*sqrt(gm/pc)
147  END IF
148 *
149 * Convert from kpc and km/sec to N-body units.
150  DO 30 k = 1,3
151  rg(k) = 1000.0*rg(k)/rbar
152  vg(k) = vg(k)/vstar
153  30 CONTINUE
154 *
155 * Define the angular velocity (z-component) and mass in N-body units.
156  r02 = rg(1)**2 + rg(2)**2
157  omega = (rg(1)*vg(2) - rg(2)*vg(1))/r02
158  tidal(4) = 2.0*omega
159  gmg = gmg/zmtot
160  gmb = gmb/zmtot
161  ar = 1000.0*ar/rbar
162 *
163  IF (gmg.GT.0.0) THEN
164  WRITE (6,35) gmg, sqrt(r02), omega, rtide, rbar
165  35 FORMAT (/,12x,'POINT-MASS MODEL: GMG =',1p,e9.1,
166  & ' RG =',e9.1,' OMEGA =',e9.1,
167  & ' RT =',0p,f6.2,' RBAR =',f6.2)
168  END IF
169  IF (gmb.GT.0.0d0) THEN
170  WRITE (6,36) gmb, ar, gam
171  36 FORMAT (/,12x,'GAMMA/ETA MODEL: GMB =',1p,e9.1,
172  & ' AR =',e9.1,' GAM =',e9.1)
173  END IF
174 *
175 * Define disk and/or logarithmic halo parameters in N-body units.
176  IF (disk.GT.0.0d0) THEN
177  disk = disk/zmtot
178  a = 1000.0*a/rbar
179  b = 1000.0*b/rbar
180  WRITE (6,40) disk, a, b
181  40 FORMAT (/,12x,'DISK MODEL: MD =',1p,e9.1,
182  & ' A =',e9.1,' B =',e9.1)
183  END IF
184 *
185 * Determine local halo velocity from total circular velocity.
186  IF (vcirc.GT.0.0d0) THEN
187  vcirc = vcirc/vstar
188  rcirc = 1000.0*rcirc/rbar
189  a2 = rcirc**2 + (a + b)**2
190  v02 = vcirc**2 - (gmg/rcirc + disk*rcirc**2/a2**1.5)
191 * Include any contribution from bulge potential (V2 = R*F).
192  IF (gmb.GT.0.0d0) THEN
193  vb2 = gmb/rcirc*(1.0 + ar/rcirc)**(gam-3.0)
194  v02 = v02 - vb2
195  END IF
196  IF (v02.LT.0.0d0) THEN
197  WRITE (6,45) v02, 0.001*rcirc*rbar
198  45 FORMAT (' ',' NEGATIVE HALO VELOCITY! V02 RCIRC ',
199  & 1p,2e10.2)
200  stop
201  END IF
202 * Specify the corresponding scale length of logarithmic halo.
203  rl2 = rcirc**2*(vcirc**2 - v02)/v02
204 * Define the asymptotic circular velocity due to halo.
205  v02 = vcirc**2
206 *
207 * Include table of circular velocity on unit #51 (km/sec & kpc).
208  ri = 1000.0/rbar
209  dr = 1000.0/rbar
210  DO 60 k = 1,30
211  ri2 = ri**2
212  a2 = ri2 + (a + b)**2
213  vb2 = gmb/ri*(1.0 + ar/ri)**(gam-3.0)
214  vcirc2 = gmg/sqrt(ri2) + disk*ri2/a2**1.5 +
215  & v02*ri2/(rl2 + ri2) + vb2
216  WRITE (52,50) sqrt(vcirc2)*vstar, ri*rbar/1000.0
217  50 FORMAT (' CIRCULAR VELOCITY: VC R ',f7.1,f7.2)
218  ri = ri + dr
219  60 CONTINUE
220  CALL flush(52)
221 *
222  a2 = r02 + (a + b)**2
223  vb2 = gmb/sqrt(r02)*(1.0 + ar/sqrt(r02))**(gam-3.0)
224  vcirc2 = gmg/sqrt(r02) + disk*r02/a2**1.5 +
225  & v02*r02/(rl2 + r02) + vb2
226  vcirc = sqrt(vcirc2)*vstar
227  WRITE (6,62) vcirc, sqrt(r02)/1000.0, sqrt(rl2)/1000.0
228  62 FORMAT (/,12x,'CIRCULAR VELOCITY: VC RG RL',f7.1,2f7.2)
229  ELSE
230  v02 = 0.0
231  END IF
232 *
233 * Initialize F & FDOT of reference frame (point-mass galaxy is OK).
234  CALL gcinit
235 *
236 * Form tidal radius from circular angular velocity (assumes apocentre).
237  IF (rtide.EQ.0.0d0) rtide = (0.5/omega**2)**0.3333
238 *
239  WRITE (6,65) (rg(k),k=1,3), (vg(k),k=1,3), sqrt(v02)
240  65 FORMAT (/,12x,'SCALED ORBIT: RG = ',3f7.2,
241  & ' VG = ',3f7.1,' V0 =',f6.1)
242  END IF
243 *
244 * Include Plummer potential for 2D and 3D (use MP = 0 if not needed).
245  IF (kz(14).EQ.3.OR.kz(14).EQ.4) THEN
246 * Check input for Plummer potential.
247  READ (5,*) mp, ap2, mpdot, tdelay
248  WRITE (6,70) mp, ap2, mpdot, tdelay
249  70 FORMAT (/,12x,'PLUMMER POTENTIAL: MP =',f7.3,' AP =',f6.2,
250  & ' MPDOT =',f7.3,' TDELAY =',f5.1)
251  mp0 = mp
252  ap2 = ap2**2
253 * Rescale velocities by including the Plummer & galactic virial energy.
254  IF (zkin.GT.0.0d0) THEN
255 * Note that QVIR = Q is saved in routine SCALE and VIR < 0 with GPU.
256  CALL energy
257  vir = pot - vir
258  qv = sqrt(qvir*vir/zkin)
259  DO 80 i = 1,n
260  DO 78 k = 1,3
261  xdot(k,i) = xdot(k,i)*qv
262  78 CONTINUE
263  80 CONTINUE
264  END IF
265  IF (rtide.EQ.0.0d0) rtide = 10.0*rscale
266  ELSE
267  mp = 0.0
268  mpdot = 0.0
269  END IF
270  rtide0 = rtide
271  tscale = tstar
272 *
273 * Define tidal radius in scaled units for linearized field.
274  IF (kz(14).LE.2) THEN
275  rtide = (zmass/tidal(1))**0.3333
276  WRITE (6,90) (tidal(k),k=1,4), tscale, rtide
277  90 FORMAT (/,12x,'TIDAL PARAMETERS: ',1p,4e10.2,
278  & ' TSCALE =',e9.2,' RTIDE =',0p,f6.2,/)
279  END IF
280 *
281  RETURN
282 *
283  END