Nbody6
edot.f
Go to the documentation of this file.
1  SUBROUTINE edot(I1,I2,J,SEMI,ECC,ECCDOT)
2 *
3 *
4 * Eccentricity derivative due to perturber.
5 * -----------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 xrel(3),vrel(3),fp(3)
9 *
10 *
11 * Set relative coordinates & velocities and initialize perturbation.
12  ri2 = 0.0
13  DO 1 k = 1,3
14  xrel(k) = x(k,i1) - x(k,i2)
15  vrel(k) = xdot(k,i1) - xdot(k,i2)
16  fp(k) = 0.0d0
17  ri2 = ri2 + xrel(k)**2
18  1 CONTINUE
19 *
20 * Obtain perturbation on KS components due to body #J.
21  l = i1
22  5 a1 = x(1,j) - x(1,l)
23  a2 = x(2,j) - x(2,l)
24  a3 = x(3,j) - x(3,l)
25  rij2 = a1**2 + a2**2 + a3**2
26  a4 = body(j)/(rij2*sqrt(rij2))
27  IF (l.EQ.i2) a4 = -a4
28  fp(1) = fp(1) + a1*a4
29  fp(2) = fp(2) + a2*a4
30  fp(3) = fp(3) + a3*a4
31  IF (l.EQ.i1) THEN
32  l = i2
33  go to 5
34  END IF
35 *
36 * Form scalar products R*V, R*F and V*F.
37  rv = 0.0
38  rf = 0.0
39  vf = 0.0
40  DO 10 k = 1,3
41  rv = rv + xrel(k)*vrel(k)
42  rf = rf + xrel(k)*fp(k)
43  vf = vf + vrel(k)*fp(k)
44  10 CONTINUE
45 *
46 * Evaluate time derivative of eccentricity (Douglas Heggie 30/8/96).
47  eccdot = (semi**2*(1.0 - ecc**2) - ri2)*vf + rv*rf
48  eccdot = eccdot/((body(i1) + body(i2))*semi*ecc)
49 *
50 * WRITE (6,15) NAME(I1), NAME(I2), ECC, ECCDOT, RV, RF, VF
51 * 15 FORMAT (' EDOT: NAM12 E ED RV RF VF ',2I6,F7.3,1P,4E10.2)
52 *
53  RETURN
54 *
55  END