SUBROUTINE STABL4(ITERM)
*
*
*       Stability test of four-body system.
*       -----------------------------------
*
      IMPLICIT REAL*8  (A-H,M,O-Z)
      LOGICAL  SWITCH,GTYPE,GTYPE0,AUX
      COMMON/CREG/  M(4),X(12),XD(12),P(12),Q(12),TIME4,ENERGY,EPSR2,
     &              XR(9),W(9),R(6),TA(6),MIJ(6),CM(10),RMAX4,TMAX,
     &              DS,TSTEP,EPS,NSTEP4,NAME4(4),KZ15,KZ27,NREG,NFN
      COMMON/TPR/   SWITCH,GTYPE,GTYPE0
      COMMON/KSAVE/  K1,K2
      REAL*8  RC(3),VC(3),RC0(3),VC0(3)
*
*
*       Transform to physical variables (retain value of SWITCH).
      AUX = SWITCH
      SWITCH = .FALSE.
      CALL ENDREG
      SWITCH = AUX
*
*       Specify indices of two least dominant bodies (denoted K3 & K4).
      K3 = 0
      DO 1 L = 1,4
          IF (L.EQ.K1.OR.L.EQ.K2) GO TO 1
          IF (K3.EQ.0) THEN
              K3 = L
          ELSE
              K4 = L
          END IF
    1 CONTINUE
*
*       Initialize scalars for orbital elements.
      RB0 = 0.0D0
      RB = 0.0D0
      RB1 = 0.0D0
      RB2 = 0.0D0
      RDOT = 0.0D0
      RDOT1 = 0.0D0
      RDOT2 = 0.0D0
      VREL2 = 0.0D0
      VREL20 = 0.0D0
      VREL21 = 0.0D0
      VREL22 = 0.0D0
*
*       Define binary masses of smallest & widest pair (K1 & K2 and K3 & K4).
      MB0 = M(K1) + M(K2)
      MB = M(K3) + M(K4)
*
*       Form separations & velocities of MB0, MB and their relative orbit.
      DO 10 K = 1,3
          J1 = 3*(K1-1) + K
          J2 = 3*(K2-1) + K
          J3 = 3*(K3-1) + K
          J4 = 3*(K4-1) + K
          RC0(K) = (M(K1)*X(J1) + M(K2)*X(J2))/MB0
          RC(K) = (M(K3)*X(J3) + M(K4)*X(J4))/MB
          VC0(K) = (M(K1)*XD(J1) + M(K2)*XD(J2))/MB0
          VC(K) = (M(K3)*XD(J3) + M(K4)*XD(J4))/MB
          RB0 = RB0 + (X(J1) - X(J2))**2
          RB = RB + (X(J3) - X(J4))**2
          RB1 = RB1 + (RC(K) - RC0(K))**2
          RB2 = RB2 + (RC0(K) - X(J3))**2
          RDOT = RDOT + (X(J3) - X(J4))*(XD(J3) - XD(J4))
          RDOT1 = RDOT1 + (RC(K) - RC0(K))*(VC(K) - VC0(K))
          RDOT2 = RDOT2 + (RC0(K) - X(J3))*(VC0(K) - XD(J3))
          VREL2 = VREL2 + (XD(J3) - XD(J4))**2
          VREL20 = VREL20 + (XD(J1) - XD(J2))**2
          VREL21 = VREL21 + (VC(K) - VC0(K))**2
          VREL22 = VREL22 + (VC0(K) - XD(J3))**2
   10 CONTINUE
*
*       Determine semi-major axis of inner binary.
      RB0 = SQRT(RB0)
      SEMI0 = 2.0/RB0 - VREL20/MB0
      SEMI0 = 1.0/SEMI0
      E0 = 1.0 - RB0/SEMI0
*
*       Form semi-major axis & eccentricity of outer pair.
      RB = SQRT(RB)
      SEMI = 2.0D0/RB - VREL2/MB
      SEMI = 1.0/SEMI
*     E = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
*
*       Evaluate orbital elements of relative c.m. motion.
      RB1 = SQRT(RB1)
      SEMI1 = 2.0D0/RB1 - VREL21/CM(7)
      SEMI1 = 1.0/SEMI1
      E1 = SQRT((1.0D0 - RB1/SEMI1)**2 + RDOT1**2/(SEMI1*CM(7)))
*
*       Consider the inner triple.
      MB2 = MB0 + M(K3)
      RB2 = SQRT(RB2)
      SEMI2 = 2.0D0/RB2 - VREL22/MB2
      SEMI2 = 1.0/SEMI2
      E2 = SQRT((1.0D0 - RB2/SEMI2)**2 + RDOT2**2/(SEMI2*MB2))
*
*       Obtain standard stability ratio (outer pericentre / inner apocentre).
      RATIO = SEMI1*(1.0D0 - E1)/(SEMI0*(1.0D0 + E0))
*
*       Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
*     AM = (2.65 + E0)*(1.0 + MB0/MB)**0.3333
*     FM = (2.0*MB0 - MB)/(3.0*MB)
*
*       Expand natural logarithm for small arguments.
*     IF (ABS(FM).LT.0.67) THEN
*         BM = FM*(1.0 - (0.5 - 0.3333*FM)*FM)
*     ELSE
*         BM = LOG(1.0D0 + FM)
*     END IF
*
*       Adopt mass dependent criterion of Harrington (A.J. 80) & Bailyn.
*     PCRIT = AM*(1.0 + 0.7*BM)*SEMI0
*
*       Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
*     Q0 = MB/MB0
*     Q1 = MAX(M(K2)/M(K1),M(K1)/M(K2))
*     Q3 = Q0**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*     PCRIT = AR*SEMI0*(1.0D0 + E0)
*
*       Check stability (AM 1997; inner triple or well separated quadruple).
      ITERM = 0
      IF (RB1.GT.5.0*RB2.AND.E2.LT.1.0) THEN
          Q1 = M(K3)/MB0
          XFAC = (1.0 + Q1)*(1.0 + E2)/SQRT(1.0 - E2)
          PCRIT = 2.8*XFAC**0.4*SEMI0
          PMIN = SEMI2*(1.0 - E2)
          IF (PCRIT.LT.PMIN) THEN
              ITERM = -1
              RATIO = SEMI2*(1.0D0 - E2)/(SEMI0*(1.0D0 + E0))
              WRITE (6,15)  SEMI0, SEMI2, E0, E2, RATIO, RB0, RB2,
     &                      PCRIT, PMIN
   15         FORMAT ('  STABT:    A0 A2 E0 E2 RATIO R0 R2 PCR PM ',
     &                             1P,2E10.2,0P,2F7.3,F6.2,1P,4E9.1)
          END IF
      ELSE IF (RB1.GT.5.0*MAX(RB0,RB).AND.E1.LT.1.0.AND.
     &         MIN(SEMI0,SEMI).GT.0.0) THEN
*       Choose smallest binary as third body and ignore fudge factor.
          IF (SEMI.GT.SEMI0) THEN
              Q1 = MB0/MB
              AIN = SEMI
          ELSE
              Q1 = MB/MB0
              AIN = SEMI0
          END IF
          XFAC = (1.0 + Q1)*(1.0 + E1)/SQRT(1.0 - E1)
          PCRIT = 2.8*XFAC**0.4*AIN
          PMIN = SEMI1*(1.0 - E1)
          IF (PCRIT.LT.PMIN) THEN
              ITERM = -1
              WRITE (6,20)  AIN, SEMI1, E0, E1, RATIO, RB0, RB1,
     &                      PCRIT, PMIN
   20         FORMAT ('  STABQ:    AIN A1 E0 E1 RATIO R0 R1 PCR PM ',
     &                             1P,2E10.2,0P,2F7.3,F6.2,1P,4E9.1)
          END IF
      END IF
*
      RETURN
*
      END