Nbody6
assess.f
Go to the documentation of this file.
1  SUBROUTINE assess(IPAIR,IM,ECC,SEMI,ITERM)
2 *
3 *
4 * Assessment of hierarchical stability.
5 * -------------------------------------
6 *
7  include 'common6.h'
8  common/binary/ cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
9  & hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
10  & namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
11  REAL*8 xx(3,3),vv(3,3)
12  SAVE itime
13  DATA itime /0/
14 *
15 *
16 * Define indices for components and c.m.
17  i1 = 2*ipair - 1
18  i2 = i1 + 1
19  icm = n + ipair
20 *
21 * Determine inclination between inner relative motion and outer orbit.
22  rv = 0.0
23  DO 10 k = 1,3
24  xx(k,1) = xrel(k,im)
25  xx(k,2) = 0.0
26  xx(k,3) = x(k,i2)
27  vv(k,1) = vrel(k,im)
28  vv(k,2) = 0.0
29  vv(k,3) = xdot(k,i2)
30  rv = rv + xrel(k,im)*vrel(k,im)
31  10 CONTINUE
32 *
33  CALL inclin(xx,vv,x(1,icm),xdot(1,icm),angle)
34 *
35 * Form inner eccentricity (neglect radial velocity for minimum).
36  rin = sqrt(xrel(1,im)**2 + xrel(2,im)**2 + xrel(3,im)**2)
37  semi0 = -0.5*body(i1)/hm(im)
38  ecc2 = (1.0 - rin/semi0)**2 + rv**2/(body(i1)*semi0)
39  ecc0 = sqrt(ecc2)
40  pmin = semi*(1.0 - ecc)
41 *
42 * Evaluate the general stability function.
43  IF (ecc.LT.1.0) THEN
44  eout = ecc
45 * Modify outer eccentricity for consistency with acceptance.
46  IF (eout.GT.0.80) THEN
47  de = 0.5*(1.0 - eout)
48  de = min(de,0.01d0)
49  eout = ecc - de
50  pmin = semi*(1.0 - eout)
51  END IF
52  nst = nstab(semi0,semi,ecc0,eout,angle,cm(1,im),
53  & cm(2,im),body(i2))
54  IF (nst.EQ.0) THEN
55  pcrit = 0.98*pmin
56  pcr = stability(cm(1,im),cm(2,im),body(i2),
57  & ecc0,ecc,angle)*semi0
58 * Reduce termination distance if old criterion < PMIN/2.
59  IF (pcr.LT.0.5*pmin) THEN
60  pcrit = 0.75*pmin
61  END IF
62  itime = itime + 1
63  IF (itime.GT.2000000000) itime = 0
64  IF (mod(itime,1000).EQ.0) THEN
65  alph = 360.0*angle/twopi
66  WRITE (6,20) ecc0, ecc, alph, semi, pcrit, pcr,
67  & r0(ipair)
68  20 FORMAT (' ASSESS E0 E1 INC A1 PCR PC0 R0 ',
69  & 2f7.3,f7.1,1p,4e9.1)
70  END IF
71  ELSE
72  pcrit = 1.01*pmin
73  END IF
74  ELSE
75  pcrit = stability(cm(1,im),cm(2,im),body(i2),
76  & ecc0,ecc,angle)*semi0
77  END IF
78 *
79 * Set new stability distance or define termination.
80  IF (pcrit.LT.pmin) THEN
81  iterm = 0
82  r0(ipair) = pcrit
83  ELSE
84  iterm = 1
85  END IF
86 *
87  RETURN
88 *
89  END