Nbody6
 All Files Functions Variables
scale.f
Go to the documentation of this file.
1  SUBROUTINE scale
2 *
3 *
4 * Scaling to new units.
5 * ---------------------
6 *
7  include 'common6.h'
8  REAL*8 rsave(3,nmax),vsave(3,nmax),bsave(nmax)
9  real*8 g, m_sun, r_sun, pc, km, kmps
10  real*8 mscale, lscale, vscale
11  SAVE rsave,vsave,bsave
12 *
13 * Define physical constants (consistent with IAU & routine UNITS).
14  g = 6.6743d-08
15  m_sun = 1.9884d+33
16  r_sun = 6.955d+10
17  pc = 3.0856776d+18
18  km = 1.0d+05
19  kmps = 1.0d+05
20 *
21 * Read virial ratio, rotation scaling factors, tidal radius & SMAX.
22  READ (5,*) q, vxrot, vzrot, rtide, smax
23 * Note RTIDE should be non-zero for isolated systems (cf. CALL LAGR).
24  rsph2 = rtide
25  qvir = q
26 *
27  zmass = 0.0d0
28  DO 10 k = 1,3
29  cmr(k) = 0.0d0
30  cmrdot(k) = 0.0d0
31  10 CONTINUE
32 *
33 * Form total mass and centre of mass displacements.
34  DO 30 i = 1,n
35  zmass = zmass + body(i)
36  DO 25 k = 1,3
37  cmr(k) = cmr(k) + body(i)*x(k,i)
38  cmrdot(k) = cmrdot(k) + body(i)*xdot(k,i)
39  25 CONTINUE
40  30 CONTINUE
41 *
42 * Adjust coordinates and velocities to c.m. rest frame.
43  DO 40 i = 1,n
44  DO 35 k = 1,3
45  x(k,i) = x(k,i) - cmr(k)/zmass
46  xdot(k,i) = xdot(k,i) - cmrdot(k)/zmass
47  35 CONTINUE
48  40 CONTINUE
49 *
50 * Check optional generation of BH binary/single and other BHs.
51  IF (kz(11).GT.0) THEN
52 * Read binary elements and component masses (SEMI = 0 for single BH).
53  READ (5,*) semi, ecc, body1, body2
54 * Convert from M_sun to pre-scaled N-body units (ZMASS = N up to now).
55  body(1) = (body1 + body2)/zmbar
56  kstar(1) = 14
57 * Place a stationary heavy binary/single BH at the centre.
58  DO 41 k = 1,3
59  x(k,1) = 0.0
60  xdot(k,1) = 0.0
61  41 CONTINUE
62  IF (semi.GT.0.0) THEN
63 * Enforce initial KS for massive binary.
64  nbin0 = 1
65  nbh0 = 2
66  kstar(2) = 14
67  ilast = kz(11) + 1
68  ELSE
69  nbh0 = 1
70  ilast = kz(11)
71  END IF
72 * Add one or more optional free-floating BH.
73  IF (kz(11).GT.1) THEN
74  DO 42 i = nbh0+1,ilast
75  READ (5,*) body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
76  body(i) = body(i)/zmbar
77  kstar(i) = 14
78  42 CONTINUE
79  nbh0 = ilast
80  END IF
81  END IF
82 *
83 * Include procedure for primordial binaries & singles from fort.10.
84  IF ((kz(22).EQ.4.OR.kz(22).EQ.-1).AND.(nbin0.GT.0)) THEN
85 * Save the binaries for re-creating two-body motion after scaling.
86  DO 45 i = 1,2*nbin0
87  bsave(i) = body(i)
88  DO 44 k = 1,3
89  rsave(k,i) = x(k,i)
90  vsave(k,i) = xdot(k,i)
91  44 CONTINUE
92  45 CONTINUE
93 * Form the c.m. of each binary for standard scaling to E = -0.25.
94  DO 47 i = 1,nbin0
95  i1 = 2*i - 1
96  i2 = 2*i
97 * Prevent over-writing BODY(I) when I = 1 and 2*I - 1 = 1.
98  b1 = body(i1)
99  b2 = body(i2)
100  zmb = b1 + b2
101  body(i) = zmb
102  DO 46 k = 1,3
103  x(k,i) = (b1*x(k,i1) + b2*x(k,i2))/zmb
104  xdot(k,i) = (b1*xdot(k,i1) +
105  & b2*xdot(k,i2))/zmb
106  46 CONTINUE
107  47 CONTINUE
108 *
109 * Move any single stars up for scaling together with c.m. particles.
110  nsing = n - 2*nbin0
111 * Note assumption ZMASS = 1 including any singles (skips DO 50 loop).
112  DO 49 i = 1,nsing
113  body(nbin0+i) = body(2*nbin0+i)
114  DO 48 k = 1,3
115  x(k,nbin0+i) = x(k,2*nbin0+i)
116  xdot(k,nbin0+i) = xdot(k,2*nbin0+i)
117  48 CONTINUE
118  49 CONTINUE
119 * Re-define N & NTOT for total energy calculation (NTOT = N for ZKIN).
120  n = nbin0 + nsing
121  ntot = n
122 *
123  write (6,491) nbin0
124  491 format (/,12x,'SCALE: ', i5,
125  & ' Binaries converted to barycentres')
126  END IF
127 * Save total number of C.Ms.
128  ncm = ntot
129 *
130 * Skip scaling of masses for unscaled upload or planetesimal disk.
131  IF (kz(22).GT.2.OR.kz(22).EQ.-1.OR.kz(5).EQ.3) go to 52
132 *
133 * Scale masses to standard units of <M> = 1/N and set total mass.
134  DO 50 i = 1,n
135  body(i) = body(i)/zmass
136  50 CONTINUE
137  zmass = 1.0
138 *
139 * Obtain the total kinetic & potential energy.
140  52 CALL energy
141 *
142 * Check option for astrophysical units.
143  if (kz(22).EQ.-1) then
144 * Save K.E. and P.E. of C.Ms for TCR & TRH (units of M_sun, pc & km/s).
145  zcm = zkin
146  pcm = pot
147 *
148 * Evaluate total energy in physical units (M_sun pc^2 s^-2).
149 * G (g^-1 cm^3 s^-2) ==> G (M_sun^-1 pc^3 s^-2)
150  g = g/(pc**3/m_sun)
151 * Form kinetic energy as KE (M_sun Km^2 s^-2) ==> KE (M_sun pc^2 s^-2).
152  zkin = zkin*km**2/pc**2
153 * PE term = (M_sun^-1 pc^3 s^-2)*(M_sun^2 pc^-1) = (M_sun pc^2 s^-2).
154  eph = zkin - g*pot
155  IF (eph.GT.0.0) eph = -eph
156 * Set conversion factors (N-body to M_sun, pc, km/s; Heggie & Mathieu).
157  mscale = zmass
158  lscale = g*mscale**2*(-0.25/eph)
159  vscale = sqrt(g*mscale/lscale)*(pc/km)
160 * Specify virial radius and mean mass (pc & M_sun).
161  rbar = lscale
162  zmbar = zmass/(n+nbin0)
163 *
164  write (6,53) mscale, lscale, vscale
165  53 format (/,12x,'SCALE: Conversion factors: MSCALE = ',f10.3,
166  & ' M_sun; LSCALE = ',f6.3,' pc; VSCALE = ',
167  & f6.3,' km/sec')
168  end if
169 *
170 * Use generalized virial theorem for external tidal field.
171  IF (kz(14).GT.0) THEN
172  az = 0.0d0
173  DO 55 i = 1,n
174  az = az + body(i)*(x(1,i)*xdot(2,i) - x(2,i)*xdot(1,i))
175  55 CONTINUE
176  IF (kz(14).EQ.1) THEN
177 * Use Chandrasekhar eq. (5.535) for virial ratio (rotating frame only).
178  vir = pot - 2.0*(etide + 0.5*tidal(4)*az)
179  ELSE
180  vir = pot - 2.0*etide
181  END IF
182  ELSE
183  vir = pot
184  END IF
185 *
186 * Allow two optional ways of skipping standard velocity scaling.
187  IF (kz(22).EQ.3.OR.kz(5).EQ.2.OR.kz(5).EQ.3) THEN
188  qv = sqrt(q*vir/zkin)
189  e0 = zkin*qv**2 - pot + etide
190  sx = 1.0
191 * Rescale velocities to new masses for two Plummer spheres.
192  IF (kz(5).EQ.2) THEN
193  zkin = 0.0
194  DO 57 i = 1,n
195  DO 56 k = 1,3
196  xdot(k,i) = xdot(k,i)*qv
197  zkin = zkin + 0.5*body(i)*xdot(k,i)**2
198  56 CONTINUE
199  57 CONTINUE
200  e0 = zkin - pot + etide
201  q = zkin/pot
202  WRITE (6,59) e0, zkin/pot
203  59 FORMAT (/,12x,'UNSCALED ENERGY E =',f10.6,
204  & ' Q =',f6.2)
205  ELSE
206  IF (kz(5).EQ.3) e0 = zkin - pot
207  WRITE (6,54) e0
208  54 FORMAT (/,12x,'UNSCALED ENERGY E =',f10.6)
209  END IF
210  ELSE IF(kz(22).NE.-1) THEN
211 * Scale non-zero velocities by virial theorem ratio.
212  IF (zkin.GT.0.0d0) THEN
213  qv = sqrt(q*vir/zkin)
214  DO 60 i = 1,n
215  DO 58 k = 1,3
216  xdot(k,i) = xdot(k,i)*qv
217  58 CONTINUE
218  60 CONTINUE
219  ELSE
220  qv = 1.0
221  END IF
222 *
223 * Scale total energy to standard units (E = -0.25 for Q < 1).
224  e0 = -0.25
225 * Include case of hot system inside reflecting boundary.
226  IF (kz(29).GT.0.AND.q.GT.1.0) e0 = 0.25
227 * ETOT = (Q - 1.0)*POT
228  etot = zkin*qv**2 - pot + etide
229 * Note that final ETOT will differ from -0.25 since ETIDE = 0.
230  IF (q.LT.1.0) THEN
231  sx = e0/etot
232  ELSE
233  sx = 1.0
234  END IF
235 *
236  WRITE (6,65) sx, etot, body(1), body(n), zmass/float(n)
237  65 FORMAT (/,12x,'SCALING: SX =',f6.2,' E =',1pe10.2,
238  & ' M(1) =',e9.2,' M(N) =',e9.2,' <M> =',e9.2)
239 *
240 * Scale coordinates & velocities to the new units.
241  DO 70 i = 1,n
242  DO 68 k = 1,3
243  x(k,i) = x(k,i)/sx
244  xdot(k,i) = xdot(k,i)*sqrt(sx)
245  68 CONTINUE
246  70 CONTINUE
247 *
248 * Set current energies to consistent values (may be needed in XTRNL0).
249  zkin = zkin*qv**2*sx
250  pot = pot*sx
251  END IF
252 *
253 * Perform second stage of the optional binary uploading procedure.
254  IF ((kz(22).EQ.4.OR.kz(22).EQ.-1).AND.(nbin0.GT.0)) THEN
255 * Place any singles last and re-create the original binaries.
256  DO 72 i = nsing,1,-1
257  l1 = nbin0
258  l2 = 2*nbin0
259  body(l2+i) = body(l1+i)
260  DO 71 k = 1,3
261  x(k,l2+i) = x(k,l1+i)
262  xdot(k,l2+i) = xdot(k,l1+i)
263  71 CONTINUE
264  72 CONTINUE
265 * Split each c.m. into binary components using saved quantities.
266  DO 74 i = 1,nbin0
267  ip = nbin0 - i + 1
268  body(2*ip-1) = bsave(2*ip-1)
269  body(2*ip) = bsave(2*ip)
270  zmb = body(2*ip-1) + body(2*ip)
271 * Use the original relative motion for unscaled two-body elements.
272  DO 73 k = 1,3
273  xrel = rsave(k,2*ip-1) - rsave(k,2*ip)
274  vrel = vsave(k,2*ip-1) - vsave(k,2*ip)
275 * Prevent over-writing of first location when IP = 1 and 2*IP-1 = 1.
276  x1 = x(k,ip)
277  v1 = xdot(k,ip)
278  x(k,2*ip-1) = x(k,ip) + bsave(2*ip)*xrel/zmb
279  x(k,2*ip) = x1 - bsave(2*ip-1)*xrel/zmb
280  xdot(k,2*ip-1) = xdot(k,ip) + bsave(2*ip)*vrel/zmb
281  xdot(k,2*ip) = v1 - bsave(2*ip-1)*vrel/zmb
282  73 CONTINUE
283  74 CONTINUE
284 * Set final values of the particle number (increase from N + NBIN0).
285  n = n + nbin0
286  nzero = n
287  ntot = n
288 *
289  write (6,741) nbin0
290  741 format (/,12x,'SCALE: ', i5, ' Binaries restored')
291  END IF
292 *
293 * Scale masses, coordinates and velocities for input in physical units.
294  if (kz(22).EQ.-1) then
295 * Re-calculate the scaled total mass.
296  zmass = 0.0d0
297  DO 743 i = 1,ntot
298  body(i) = body(i)/mscale
299  zmass = zmass + body(i)
300  DO 742 k = 1,3
301  x(k,i) = x(k,i)/lscale
302  xdot(k,i) = xdot(k,i)/vscale
303  742 continue
304  743 continue
305 *
306 * Perform energy check after scaling.
307  CALL energy
308  write (6,745) zkin - pot
309  745 format (/,12x,'SCALE: Scaled total energy: ',f8.5)
310 * Set energy of singles & C.Ms only in N-body units for TCR & TRH.
311 * escale = mscale*(vscale**2)
312  sx = 1.0
313  zkin = zcm/(mscale*vscale**2)
314  pot = pcm*lscale/mscale**2
315  e0 = zkin - pot
316  write (6,746) e0
317  746 format (/,12x,'SCALE: Scaled C.M. energy: ',f8.5)
318  end if
319 *
320 * Introduce optional BH binary (SEMI in scaled N-body units).
321  IF (kz(11).GT.0) THEN
322  IF (semi.GT.0.0) THEN
323  zmb = body(1)
324 * Split c.m. into components with given mass ratio.
325  body(1) = body1/(body1 + body2)*zmb
326  body(2) = body2/(body1 + body2)*zmb
327  rap = semi*(1.0 + ecc)
328  vap = sqrt(zmb*(1.0 - ecc)/(semi*(1.0 + ecc)))
329 * Specify two-body motion for SEMI & ECC in x-y plane.
330  x(1,2) = x(1,1) - body(1)*rap/zmb
331  x(1,1) = x(1,1) + body(2)*rap/zmb
332  xdot(2,2) = xdot(2,1) - body(1)*vap/zmb
333  xdot(2,1) = xdot(2,1) + body(2)*vap/zmb
334  DO 75 i = 1,2
335  x(2,i) = 0.0
336  x(3,i) = 0.0
337  xdot(1,i) = 0.0
338  xdot(3,i) = 0.0
339  75 CONTINUE
340  END IF
341  DO 77 i = 1,ilast
342  WRITE (6,76) body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
343  76 FORMAT (12x,'BH COMPONENTS M X XDOT ',
344  & 1p,e10.2,0p,2(2x,3f7.3))
345  77 CONTINUE
346  END IF
347 *
348 * Check whether to include rotation (VXROT = 0 in standard case).
349  IF (vxrot.GT.0.0d0) THEN
350 *
351 * Set angular velocity for retrograde motion (i.e. star clusters).
352  omega = -sx*sqrt(zmass*sx)
353  WRITE (6,78) vxrot, vzrot, omega
354  78 FORMAT (/,12x,'VXROT =',f6.2,' VZROT =',f6.2,
355  & ' OMEGA =',f7.2)
356 *
357 * Add solid-body rotation about Z-axis (reduce random velocities).
358  DO 80 i = 1,n
359  xdot(1,i) = xdot(1,i)*vxrot - x(2,i)*omega
360  xdot(2,i) = xdot(2,i)*vxrot + x(1,i)*omega
361  xdot(3,i) = xdot(3,i)*vzrot
362  80 CONTINUE
363  END IF
364 *
365 * Check option for writing the initial conditions on unit #10.
366  IF (kz(22).EQ.1) THEN
367  DO 85 i = 1,n
368  WRITE (10,84) body(i), (x(k,i),k=1,3), (xdot(k,i),k=1,3)
369  84 FORMAT (1p,7e14.6)
370  85 CONTINUE
371  END IF
372 *
373 * Check for writing initial conditions (physical units) on unit #99.
374  IF (kz(22).EQ.-1.OR.kz(22).EQ.5) THEN
375  DO 86 i = 1,n
376  write (99,84) body(i)*mscale, (x(k,i)*lscale,k=1,3),
377  & (xdot(k,i)*vscale,k=1,3)
378  86 CONTINUE
379  END IF
380 *
381 * Check option for reading initial subsystems (solar masses).
382  IF (kz(24).GT.0) THEN
383  k = kz(24)
384  DO 90 i = 1,k
385  READ (5,*) body(i), (x(j,i),j=1,3), (xdot(j,i),j=1,3)
386  body(i) = body(i)/(zmbar*float(n))
387  90 CONTINUE
388  END IF
389 *
390 * Set initial crossing time in scaled units.
391  tcr = zmass**2.5/(2.0d0*abs(e0))**1.5
392  tcr0 = tcr
393 *
394 * Obtain approximate half-mass radius after scaling.
395  rscale = 0.5*zmass**2/pot
396 * Set square radius of reflecting sphere (used with option 29).
397  rsph2 = (rsph2*rscale)**2
398 * Form equilibrium rms velocity (temporarily defined as VC).
399  vc = sqrt(2.0d0*abs(e0)/zmass)
400 *
401 * Check for general binary search of initial condition.
402  IF (kz(4).GT.0) THEN
403  CALL evolve(0,0)
404  END IF
405 *
406 * Print half-mass relaxation time & equilibrium crossing time.
407  a1 = float(ncm)
408  trh = 4.0*twopi/3.0*(vc*rscale)**3/(15.4*zmass**2*log(a1)/a1)
409  WRITE (6,95) trh, tcr, 2.0*rscale/vc, smax
410  95 FORMAT (/,12x,'TIME SCALES: TRH =',1p,e8.1,' TCR =',0p,f6.2,
411  & ' 2<R>/<V> =',f6.2,' SMAX =',f7.3,/)
412 *
413  RETURN
414 *
415  END