orbit.f

Go to the documentation of this file.

      SUBROUTINE orbit(I,J,SEMI,ECC,GI)
*
*
*       Two-body elements.
*       -----------------
*
      include 'common6.h'
*
*
*       Find the dominant neighbour if body #J not specified.
      gi = 0.0
      IF (j.LE.0) THEN
          nnb = list(1,i)
          fmax = 1.0d-10
          jm = list(2,i)
          DO 2 l = 1,nnb
              jj = list(l+1,i)
              jlist(l) = jj
              rij2 = 0.0
              DO 1 k = 1,3
                  rij2 = rij2 + (x(k,i) - x(k,jj))**2
    1         CONTINUE
              fij = (body(i) + body(jj))/rij2
*       Exclude any c.m. bodies from dominant motion.
              IF (fij.GT.fmax.AND.jj.LE.n) THEN
                  fmax = fij
                  jm = jj
              END IF
    2     CONTINUE
*       Avoid rare case of halo orbit with zero neighbour number.
          IF (jm.EQ.i) jm = i - 1
*
*       Obtain the relative perturbation for decision-making.
          CALL fpert(i,jm,nnb,pert)
          gi = pert/fmax
          j = jm
      END IF
*
*       Determine the semi-major axis and eccentricity.
      rij2 = 0.0
      vij2 = 0.0
      rdot = 0.0
      DO 5 k = 1,3
          rij2 = rij2 + (x(k,i) - x(k,j))**2
          vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
          rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) - xdot(k,j))
    5 CONTINUE
      rij = sqrt(rij2)
      a1 = 2.0/rij - vij2/(body(i) + body(j))
      semi = 1.0/a1
      ecc2 = (1.0 - rij/semi)**2 + rdot**2/(semi*(body(i) + body(j)))
      ecc = sqrt(ecc2)
*
      RETURN
*
      END