Nbody6
gcint.f
Go to the documentation of this file.
1  SUBROUTINE gcint
2 *
3 *
4 * Integration of 3D cluster orbit.
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  REAL*8 fm(3),fd(3),fs(3),fsd(3)
11 *
12 *
13 * Predict coordinates and velocities to order FDOT.
14  dt = time + toff - tg
15 * Note: integration step may exceed STEPX (depends on block-step).
16  dt2 = 0.5*dt
17  dt3 = one3*dt
18  rg2 = 0.0
19  rgvg = 0.0
20  DO 10 k = 1,3
21  rg(k) = ((fgd(k)*dt3 + fg(k))*dt2 + vg(k))*dt + rg(k)
22  vg(k) = (fgd(k)*dt2 + fg(k))*dt + vg(k)
23  rg2 = rg2 + rg(k)**2
24  rgvg = rgvg + rg(k)*vg(k)
25  10 CONTINUE
26 *
27 * Obtain force and first derivative of point-mass galaxy.
28  IF (kz(14).EQ.3.AND.gmg.GT.0.0d0) THEN
29  CALL fnuc(rg,vg,fm,fd)
30  ELSE
31  DO 15 k = 1,3
32  fm(k) = 0.0
33  fd(k) = 0.0
34  15 CONTINUE
35  END IF
36 *
37 * Check bulge force.
38  IF (gmb.GT.0.0d0) THEN
39  CALL fbulge(rg,vg,fs,fsd)
40  DO 20 k = 1,3
41  fm(k) = fm(k) + fs(k)
42  fd(k) = fd(k) + fsd(k)
43  20 CONTINUE
44  END IF
45 *
46 * Include optional Miyamoto disk component (Book eq. 8.52).
47  IF (disk.GT.0.0d0) THEN
48  CALL fdisk(rg,vg,fs,fsd)
49  DO 25 k = 1,3
50  fm(k) = fm(k) + fs(k)
51  fd(k) = fd(k) + fsd(k)
52  25 CONTINUE
53  END IF
54 *
55 * Check addition of logarithmic galaxy potential (not linearized).
56  IF (v02.GT.0.0d0) THEN
57  CALL fhalo(rg,vg,fs,fsd)
58  DO 30 k = 1,3
59  fm(k) = fm(k) + fs(k)
60  fd(k) = fd(k) + fsd(k)
61  30 CONTINUE
62  END IF
63 *
64 * Set time factors for corrector.
65  dt13 = one3*dt
66  dtsq12 = one12*dt**2
67  tg = tg + dt
68 *
69 * Include the Hermite corrector and update F & FDOT.
70  DO 40 k = 1,3
71  df = fg(k) - fm(k)
72  sum = fgd(k) + fd(k)
73  at3 = 2.0*df + sum*dt
74  bt2 = -3.0*df - (sum + fgd(k))*dt
75  rg(k) = rg(k) + (0.6*at3 + bt2)*dtsq12
76  vg(k) = vg(k) + (0.75*at3 + bt2)*dt13
77  fg(k) = fm(k)
78  fgd(k) = fd(k)
79  40 CONTINUE
80 *
81 * Update angular velocity in case of non-circular orbit.
82 * OM1 = (RG(2)*VG(3) - RG(3)*VG(2))/RG2
83 * OM2 = (RG(3)*VG(1) - RG(1)*VG(3))/RG2
84  om3 = (rg(1)*vg(2) - rg(2)*vg(1))/rg2
85 * OMEGA2 = OM1**2 + OM2**2 + OM3**2
86  omega2 = om3**2
87  omega = sqrt(omega2)
88  tidal(4) = 2.0*omega
89 *
90  RETURN
91 *
92  END