chaos0.f

Go to the documentation of this file.

      SUBROUTINE chaos0(QPERI,ECC,EB0,ZJ0,M1,M2,S1,S2,W,ECRIT,AR,BR,
     &                                                           idis)
*
*
*       Initial chaos boundary parameters.
*       ----------------------------------
*
*       Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      IMPLICIT REAL*8 (a-h,m,o-z)
      REAL*8  w(4)
      DATA  c0,c1,c2,c3,c4,c5
     &   /-1.763e-03,0.1277,-3.743e-03,3.525e-03,-6.179e-04,2.744e-05/
      DATA  cb /1.836/
      DATA  a0,a1,a2,a3 /0.079,-0.022,1.275,-0.670/
*
*
*       Form chaos boundary parameter and decide which star used for scaling.
      p = ((m2/m1)*w(2)/w(1))**0.4*(s1/s2)
      IF (p.GT.1.0) THEN
          s = s1
          m = m1
          m21 = m2/m1
          w2 = w(1)  
      ELSE
          s = s2
          m = m2
          m21 = m1/m2
          w2 = w(2)  
      END IF
*
*       Specify energy & angular momentum mass factors.
      ce = 0.5*m1*m2
      zmu = m1*m2/(m1 + m2)
      cj = zmu*sqrt(m1 + m2)
*
*       Define RP.
      rp = 0.5*qperi*(1.0 + ecc)
*
*       Write polynomial for E(c).
      y = log((2*rp/(cb*s))**5*(w2/(1.0 + m21))**2)
      ec = ((((c5*y + c4)*y + c3)*y + c2)*y + c1)*y + c0
*
*       Compare EC to ECDIS (disruption value for EC).
      iout = idis
      idis = 0
      ecdis = ((a3*ecc + a2)*ecc + a1)*ecc + a0
      IF (ec.LT.ecdis) THEN
          idis = 1
      END IF
*
*   Check whether spiral is indicated (unless collision or ECC > 1).
      IF (ecc.LT.ec.AND.idis.EQ.0.AND.ecc.LT.1.0) THEN
          idis = -1
      RETURN
      END IF
*
*       Obtain R(c) & Ecrit.
      rc = 2.0*rp/(1.0 + ec)
      ecrit = ce*(ec - 1.0)/rc
*
*       Calculate RP(b) & E(b).
      alf2 = (2.0/sqrt(w2))*sqrt(s**3/m)
      cfac = (2.0*ecrit - eb0)/ce
      yfac = (zj0 + 2.0*alf2*(ecrit - eb0))/cj
      tmp = 1.0 + cfac*yfac**2
      tmp = max(tmp,0.0d0)
      eb = sqrt(tmp)
      rb = (eb - 1.0)*ce/(2.0*ecrit - eb0)
*
*       Form the A & B parameters for determining minimum eccentricity.
      ar = (rb - rc)/(eb - ec)
      br = (eb*rc - ec*rb)/(eb - ec)
*       Note: AR & BR in units of length and ECRIT in mass/radius.
*
*       ---------------------------------------------------------------
*       Include occasional diagnostic output (skip CHRECT/CHAOS).
      IF (iout.GE.0) THEN
          WRITE (6,1)  m1, m2, m21, qperi/s, ec, ecdis, eb
    1     FORMAT (' INIT:    M1 M2 M21 QP/S EC ECD EB ',
     &                       1p,2e10.2,0p,2f5.1,3f6.2)
      END IF
*       ---------------------------------------------------------------
      RETURN
*
      END