Nbody6
 All Files Functions Variables
hut.f
Go to the documentation of this file.
1  subroutine hut(es0,spin10,spin20,ecc,spin1,spin2,nsteps,dtau)
2 
3  implicit real*8 (a-h,o-z)
4  real*8 u(3),udot(3)
5 
6  u(1)=es0
7  u(2)=spin10
8  u(3)=spin20
9 
10  call deriv2(u,udot)
11 
12 * include step reduction for large de/dt or primary spin rate.
13  it = 0
14  1 IF (abs(udot(1))*dtau.GT.0.01*max(es0,0.01d0)) THEN
15  dtau = 0.5d0*dtau
16  nsteps = 2*nsteps
17  it = it + 1
18  IF (it.LT.10) go to 1
19  END IF
20 
21  2 IF (abs(udot(2))*dtau.GT.0.01*u(2).AND.u(2).GT.0.1) THEN
22  dtau = 0.5d0*dtau
23  nsteps = 2*nsteps
24  it = it + 1
25  IF (it.LT.10) go to 2
26  END IF
27 
28 * treat large eccentricity carefully for rapid spin change.
29  IF (es0.GT.0.9995.and.
30  & abs(udot(2))*dtau.GT.0.01*u(2).AND.u(2).GT.0.1) THEN
31  dtau = 0.5d0*dtau
32  nsteps = 2*nsteps
33  it = it + 1
34  END IF
35 
36  IF (it.GT.0) THEN
37  WRITE (96,3) nsteps, it, u, udot, dtau
38  3 FORMAT (' HUT REDUCE # IT u ud dt ',i6,i3,f8.4,1p,7e9.1)
39  CALL flush(96)
40  END IF
41 
42  do i=1,nsteps
43 * Save spins in case eccentricity goes negative.
44  usave2 = u(2)
45  usave3 = u(3)
46  call rk4b(dtau,u)
47 * note there are 4 calls to deriv2 when u(1) may go negative.
48  IF (u(1).lt.0.002.or.u(1).gt.0.99999) then
49 * enforce circularization by adopting e=0.00199 and copying spins.
50  IF (u(1).LT.0.0) WRITE (6,4) i, nsteps, u(1)
51  4 FORMAT (' HUT SAFETY EXIT! I nsteps u1 ',2i5,f8.4)
52  u(1) = 0.00199
53  u(2) = usave2
54  u(3) = usave3
55  go to 5
56  END IF
57  enddo
58 
59  5 ecc=u(1)
60  spin1=u(2)
61  spin2=u(3)
62 
63  end
64 
65  subroutine rk4b(dt,u)
66 
67 * runge-kutta integrator.
68 * -----------------------
69 
70 * author: rosemary mardling(3/98).
71 
72  parameter(n=3)
73  implicit real*8 (a-h,o-z)
74  real*8 u0(n),ut(n),du(n),u(n)
75  real*8 a(4),b(4)
76 
77  a(1)=dt/2.0d0
78  a(2)=a(1)
79  a(3)=dt
80  a(4)=0.0d0
81 
82  b(1)=dt/6.0d0
83  b(2)=dt/3.0d0
84  b(3)=b(2)
85  b(4)=b(1)
86 
87  do i=1,n
88  u0(i)=u(i)
89  ut(i)=0.0d0
90  enddo
91 
92  do j=1,4
93  call deriv2(u,du)
94 
95  do i=1,n
96  u(i)=u0(i)+a(j)*du(i)
97  ut(i)=ut(i)+b(j)*du(i)
98  enddo
99  enddo
100 
101  do i=1,n
102  u(i)=u0(i)+ut(i)
103  enddo
104 
105  end
106 
107  subroutine deriv2(u,udot)
108 
109  implicit real*8 (a-h,m,o-z)
110  real*8 u(3),udot(3)
111  common/spins/angmom0,rg2(2),m21,r21,semi0,c1,c2,c3,c4,c5,semi
112 * SAVE ic
113 * DATA ic /0/
114 
115  e=u(1)
116  spin1=u(2)
117  spin2=u(3)
118  semi1=semi
119 
120  e2=e**2
121  e4=e2**2
122  e6=e4*e2
123  fac=1-e2
124 
125  f2=1+7.5*e2+5.625*e4+0.3125*e6
126  f3=1+3.75*e2+1.875*e4+0.078125*e6
127  f4=1+1.5*e2+0.125*e4
128  f5=1+3*e2+0.375*e4
129 
130  semi=angmom0-rg2(1)*spin1-m21*r21**2*rg2(2)*spin2
131  semi=(semi*(1+m21)/m21/semi0**2)**2/fac
132  oa = 1.0/semi
133 
134  if(e.le.0.0.or.e.ge.1.0.or.oa.lt.0.0)then
135  udot(1) = 0.0
136  udot(2) = 0.0
137  udot(3) = 0.0
138  else
139  udot(1)=-oa**8*(e/fac**6.5)*(
140  & c1*(f3-(11./18.)*fac**1.5*f4*spin1/oa**1.5)+
141  & c2*(f3-(11./18.)*fac**1.5*f4*spin2/oa**1.5))
142  udot(2)=(oa/fac)**6*c3*(oa**1.5*f2-fac**1.5*f5*spin1)
143  udot(3)=(oa/fac)**6*c4*(oa**1.5*f2-fac**1.5*f5*spin2)
144  endif
145  if (e.lt.0.002.and.semi.lt.0.0d0) semi=semi1
146 
147 * ic = ic + 1
148 * IF (ic.EQ.1.OR.mod(ic,10000).EQ.0) THEN
149 * WRITE (6,1) ic, e, spin1, spin2, udot
150 * 1 FORMAT (' HUT DERIV # e s1 s2 udot ',i10,f8.4,1p,5e9.1)
151 * END IF
152 
153  end