Revision dd8daa98413988b11c594e5a3bff8f5659d7ecf7 authored by Florent Renaud on 26 June 2014, 09:41:19 UTC, committed by Florent Renaud on 26 June 2014, 09:41:19 UTC
Originally released on 18 March 2013
Based on Nbody6 downloaded on 29 January 2013
1 parent edab1b7
Raw File
erel.f
      SUBROUTINE EREL(IM,EB,SEMI)
*
*
*       Dominant two-body energy in chain system.
*       -----------------------------------------
*
*     INCLUDE 'COMARC.CH'
      INCLUDE 'ARCCOM2e2.CH'
*
*
*       Determine indices of dominant motion.
      RX = 0.0
      DO 1 K = 1,N-1
          IF (RINV(K).GT.RX) THEN
              RX = RINV(K)
              IM = K
          END IF
    1 CONTINUE
      K1 = INAME(IM)
      K2 = INAME(IM+1)
*
*       Form kinetic energy terms of dominant c.m. (K1 + K2) and the rest.
      ZK = 0.0D0
      P0 = 0.0D0
      DO 10 K = 1,3
          PS = 0.0
*       Exclude dominant bodies using names of chain members.
          DO 5 I = 1,N
              IF (INAME(I).NE.K1.AND.INAME(I).NE.K2) THEN
                  J = 3*(I - 1)
                  PS = PS + M(I)*V(J+K)
                  ZK = ZK + M(I)*V(J+K)**2
              END IF
    5     CONTINUE
          P0 = P0 + PS**2
   10 CONTINUE
*
*       Obtain potential energy due to non-dominant chain distances.
      POT = 0.0D0
      DO 20 I = 1,N-1
          LI = 3*(I - 1)
          DO 15 J = I+1,N
              IF (I.EQ.K1.AND.J.EQ.K2) GO TO 15
              LJ = 3*(J - 1)
              RIJ2 = 0.0
              DO 12 K = 1,3
                  RIJ2 = RIJ2 + (X(LI+K) - X(LJ+K))**2
   12         CONTINUE
              POT = POT + M(I)*M(J)/SQRT(RIJ2)
   15     CONTINUE
   20 CONTINUE
*
*       Derive binding energy from total energy and perturbing function.
      MB = M(K1) + M(K2)
      VP = 0.5D0*(P0/MB + ZK) - POT
      EB = ENERGY - VP
*
*       Set semi-major axis.
      SEMI = -M(K1)*M(K2)/(2.0D0*EB)
*
      RETURN
*
      END
back to top