edot.f

Go to the documentation of this file.

      SUBROUTINE edot(I1,I2,J,SEMI,ECC,ECCDOT)
*
*
*       Eccentricity derivative due to perturber.
*       -----------------------------------------
*
      include 'common6.h'
      REAL*8  xrel(3),vrel(3),fp(3)
*
*
*       Set relative coordinates & velocities and initialize perturbation.
      ri2 = 0.0
      DO 1 k = 1,3
          xrel(k) = x(k,i1) - x(k,i2)
          vrel(k) = xdot(k,i1) - xdot(k,i2)
          fp(k) = 0.0d0
          ri2 = ri2 + xrel(k)**2
    1 CONTINUE
*
*       Obtain perturbation on KS components due to body #J.
      l = i1
    5 a1 = x(1,j) - x(1,l)
      a2 = x(2,j) - x(2,l)
      a3 = x(3,j) - x(3,l)
      rij2 = a1**2 + a2**2 + a3**2
      a4 = body(j)/(rij2*sqrt(rij2))
      IF (l.EQ.i2) a4 = -a4
      fp(1) = fp(1) + a1*a4
      fp(2) = fp(2) + a2*a4
      fp(3) = fp(3) + a3*a4
      IF (l.EQ.i1) THEN
          l = i2
          go to 5
      END IF
*
*       Form scalar products R*V, R*F and V*F.
      rv = 0.0
      rf = 0.0
      vf = 0.0
      DO 10 k = 1,3
          rv = rv + xrel(k)*vrel(k)
          rf = rf + xrel(k)*fp(k)
          vf = vf + vrel(k)*fp(k)
   10 CONTINUE
*
*       Evaluate time derivative of eccentricity (Douglas Heggie 30/8/96).
      eccdot = (semi**2*(1.0 - ecc**2) - ri2)*vf + rv*rf
      eccdot = eccdot/((body(i1) + body(i2))*semi*ecc)
*
*     WRITE (6,15)  NAME(I1), NAME(I2), ECC, ECCDOT, RV, RF, VF
*  15 FORMAT (' EDOT:    NAM12 E ED RV RF VF ',2I6,F7.3,1P,4E10.2)
*
      RETURN
*
      END