Nbody6
 All Files Functions Variables
chaos0.f
Go to the documentation of this file.
1  SUBROUTINE chaos0(QPERI,ECC,EB0,ZJ0,M1,M2,S1,S2,W,ECRIT,AR,BR,
2  & idis)
3 *
4 *
5 * Initial chaos boundary parameters.
6 * ----------------------------------
7 *
8 * Theory of Rosemary Mardling, Ap. J. XX, YYY, 1995.
9 * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
10 *
11  IMPLICIT REAL*8 (a-h,m,o-z)
12  REAL*8 w(4)
13  DATA c0,c1,c2,c3,c4,c5
14  & /-1.763e-03,0.1277,-3.743e-03,3.525e-03,-6.179e-04,2.744e-05/
15  DATA cb /1.836/
16  DATA a0,a1,a2,a3 /0.079,-0.022,1.275,-0.670/
17 *
18 *
19 * Form chaos boundary parameter and decide which star used for scaling.
20  p = ((m2/m1)*w(2)/w(1))**0.4*(s1/s2)
21  IF (p.GT.1.0) THEN
22  s = s1
23  m = m1
24  m21 = m2/m1
25  w2 = w(1)
26  ELSE
27  s = s2
28  m = m2
29  m21 = m1/m2
30  w2 = w(2)
31  END IF
32 *
33 * Specify energy & angular momentum mass factors.
34  ce = 0.5*m1*m2
35  zmu = m1*m2/(m1 + m2)
36  cj = zmu*sqrt(m1 + m2)
37 *
38 * Define RP.
39  rp = 0.5*qperi*(1.0 + ecc)
40 *
41 * Write polynomial for E(c).
42  y = log((2*rp/(cb*s))**5*(w2/(1.0 + m21))**2)
43  ec = ((((c5*y + c4)*y + c3)*y + c2)*y + c1)*y + c0
44 *
45 * Compare EC to ECDIS (disruption value for EC).
46  iout = idis
47  idis = 0
48  ecdis = ((a3*ecc + a2)*ecc + a1)*ecc + a0
49  IF (ec.LT.ecdis) THEN
50  idis = 1
51  END IF
52 *
53 * Check whether spiral is indicated (unless collision or ECC > 1).
54  IF (ecc.LT.ec.AND.idis.EQ.0.AND.ecc.LT.1.0) THEN
55  idis = -1
56  RETURN
57  END IF
58 *
59 * Obtain R(c) & Ecrit.
60  rc = 2.0*rp/(1.0 + ec)
61  ecrit = ce*(ec - 1.0)/rc
62 *
63 * Calculate RP(b) & E(b).
64  alf2 = (2.0/sqrt(w2))*sqrt(s**3/m)
65  cfac = (2.0*ecrit - eb0)/ce
66  yfac = (zj0 + 2.0*alf2*(ecrit - eb0))/cj
67  tmp = 1.0 + cfac*yfac**2
68  tmp = max(tmp,0.0d0)
69  eb = sqrt(tmp)
70  rb = (eb - 1.0)*ce/(2.0*ecrit - eb0)
71 *
72 * Form the A & B parameters for determining minimum eccentricity.
73  ar = (rb - rc)/(eb - ec)
74  br = (eb*rc - ec*rb)/(eb - ec)
75 * Note: AR & BR in units of length and ECRIT in mass/radius.
76 *
77 * ---------------------------------------------------------------
78 * Include occasional diagnostic output (skip CHRECT/CHAOS).
79  IF (iout.GE.0) THEN
80  WRITE (6,1) m1, m2, m21, qperi/s, ec, ecdis, eb
81  1 FORMAT (' INIT: M1 M2 M21 QP/S EC ECD EB ',
82  & 1p,2e10.2,0p,2f5.1,3f6.2)
83  END IF
84 * ---------------------------------------------------------------
85  RETURN
86 *
87  END