assess.f

Go to the documentation of this file.

      SUBROUTINE assess(IPAIR,IM,ECC,SEMI,ITERM)
*
*
*       Assessment of hierarchical stability.
*       -------------------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
      REAL*8  xx(3,3),vv(3,3)
      SAVE itime
      DATA itime /0/
*
*
*       Define indices for components and c.m.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      icm = n + ipair
*
*       Determine inclination between inner relative motion and outer orbit.
      rv = 0.0
      DO 10 k = 1,3
          xx(k,1) = xrel(k,im)
          xx(k,2) = 0.0
          xx(k,3) = x(k,i2)
          vv(k,1) = vrel(k,im)
          vv(k,2) = 0.0
          vv(k,3) = xdot(k,i2)
          rv = rv + xrel(k,im)*vrel(k,im)
   10 CONTINUE
*
      CALL inclin(xx,vv,x(1,icm),xdot(1,icm),angle)
*
*       Form inner eccentricity (neglect radial velocity for minimum).
      rin = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
      semi0 = -0.5*body(i1)/hm(im)
      ecc2 = (1.0 - rin/semi0)**2 + rv**2/(body(i1)*semi0)
      ecc0 = sqrt(ecc2)
      pmin = semi*(1.0 - ecc)
*
*       Evaluate the general stability function.
      IF (ecc.LT.1.0) THEN
          eout = ecc
*       Modify outer eccentricity for consistency with acceptance.
          IF (eout.GT.0.80) THEN
              de = 0.5*(1.0 - eout)
              de = min(de,0.01d0)
              eout = ecc - de
              pmin = semi*(1.0 - eout)
          END IF
          nst = nstab(semi0,semi,ecc0,eout,angle,cm(1,im),
     &                                    cm(2,im),body(i2))
          IF (nst.EQ.0) THEN
              pcrit = 0.98*pmin
              pcr = stability(cm(1,im),cm(2,im),body(i2),
     &                                          ecc0,ecc,angle)*semi0
*        Reduce termination distance if old criterion < PMIN/2.
              IF (pcr.LT.0.5*pmin) THEN
                  pcrit = 0.75*pmin
              END IF
              itime = itime + 1
              IF (itime.GT.2000000000) itime = 0
              IF (mod(itime,1000).EQ.0) THEN
                  alph = 360.0*angle/twopi
                  WRITE (6,20)  ecc0, ecc, alph, semi, pcrit, pcr,
     &                          r0(ipair)
   20             FORMAT (' ASSESS    E0 E1 INC A1 PCR PC0 R0 ',
     &                                2f7.3,f7.1,1p,4e9.1)
              END IF
          ELSE
              pcrit = 1.01*pmin
          END IF
      ELSE
          pcrit = stability(cm(1,im),cm(2,im),body(i2),
     &                                        ecc0,ecc,angle)*semi0
      END IF
*
*       Set new stability distance or define termination.
      IF (pcrit.LT.pmin) THEN
          iterm = 0
          r0(ipair) = pcrit
      ELSE
          iterm = 1
      END IF
*
      RETURN
*
      END