Nbody6
 All Files Functions Variables
rkint.f
Go to the documentation of this file.
1  subroutine rkint(dt,u)
2 
3 * Runge-Kutta integrator.
4 * -----------------------
5 
6 * Author: Rosemary Mardling (3/98).
7 
8  parameter(n=6)
9  implicit real*8 (a-h,o-z)
10  real*8 u0(n),ut(n),du(n),u(n)
11  real*8 a(4),b(4)
12 
13  iq = 0
14  a(1)=dt/2.0d0
15  a(2)=a(1)
16  a(3)=dt
17  a(4)=0.0d0
18 
19  b(1)=dt/6.0d0
20  b(2)=dt/3.0d0
21  b(3)=b(2)
22  b(4)=b(1)
23 
24  do i=1,n
25  u0(i)=u(i)
26  ut(i)=0.0d0
27  enddo
28 
29  do j=1,4
30  call deriv(u,du,iq)
31 
32  do i=1,n
33  u(i)=u0(i)+a(j)*du(i)
34  ut(i)=ut(i)+b(j)*du(i)
35  enddo
36  enddo
37 
38  if(iq.gt.0)then
39  do i=1,n
40  ut(i)=0.0d0
41  enddo
42  endif
43 
44  do i=1,n
45  u(i)=u0(i)+ut(i)
46  enddo
47 
48  end
49 
50 
51  subroutine deriv(u,du,iq)
52 
53 * Differential equations for hierarchical binary.
54 * ----------------------------------------------
55 
56 * Author: Rosemary Mardling (3/98).
57 
58  implicit real*8 (a-h,m,o-z)
59  common/rksave/ coeff,hohat(3),e,a,hh,mb
60  common/tidal/ cq(2),ct(2),cgr,dedt
61  real*8 ehat(3),hhat(3),qhat(3),edot(3),hdot(3),u(6),du(6)
62  SAVE itime,igr
63  DATA itime,igr /0,0/
64 
65 
66 * Save initial basic elements.
67  e0 = e
68  a0 = a
69  h0 = hh
70 
71 * Update eccentricity and angular momentum.
72  e=sqrt(u(1)**2+u(2)**2+u(3)**2)
73  if(e.ge.1.0)then
74  iq = 1
75  goto 10
76  endif
77  hh=sqrt(u(4)**2+u(5)**2+u(6)**2)
78  a=hh**2/mb/(1.0-e**2)
79 
80 * Define unit vectors.
81  do i=1,3
82  ehat(i)=u(i)/e
83  hhat(i)=u(i+3)/hh
84  enddo
85 
86 * Calculate unit vector orthogonal to ehat and hhat (Peter's q).
87 
88  call cross(hhat,ehat,qhat)
89 
90 * Calculate components of Peter's Sij tensor (third body average).
91 
92  s11=-coeff*(3*(dot(hohat,ehat))**2 - 1)
93  s12=-3*coeff*dot(hohat,ehat)*dot(hohat,qhat)
94  s13=-3*coeff*dot(hohat,ehat)*dot(hohat,hhat)
95  s22=-coeff*(3*(dot(hohat,qhat))**2 - 1)
96  s23=-3*coeff*dot(hohat,qhat)*dot(hohat,hhat)
97 
98 * Calculate rate of change of evec and hvec.
99 
100  do i=1,3
101  edot(i)=(e*a*hh/2/mb)*
102  & (-5*s12*ehat(i)+(4*s11-s22)*qhat(i)-s23*hhat(i))
103 
104  hdot(i)=(a**2/2)*((1-e**2)*s23*ehat(i)
105  & -(1+4*e**2)*s13*qhat(i)+5*e**2*s12*hhat(i))
106  enddo
107 
108  ekd = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
109  tav = e/sqrt(edot(1)**2+edot(2)**2+edot(3)**2)
110 *
111 * Adopt eccentricity functions from Hut theory (A & A 99, 126, 1981).
112  ge=30*e+45*e**3+3.75*e**5
113  ge=0.5*ge
114  he=(1.0+3.75*e**2+1.875*e**4+5.0/64.0*e**6)/(1.0-e**2)**6.5
115 
116 * Specify scaling factors for quadrupole and tidal interaction.
117  slr=a*(1.0-e**2)
118  zq=ge/(slr**5*a*sqrt(a))
119  zz=he/a**8
120  zt=9.0*e*zz*(ct(1) + ct(2))
121  zq=zq*(cq(1) + cq(2))
122 * Correct for wrong eccentricity dependence.
123  zq=zq*(1.0 - e**2)
124 * Note that angular momentum derivative is of same form as de/dt.
125 * zh=hh*zz*(ct(1) + ct(2))
126 *
127 * Form scaled expression for GR precession.
128  zg = cgr/(slr*a*sqrt(a))
129 
130 * zq = 0.0
131 * zg = 0.0
132 * zt = 0.0
133 * zh = 0.0
134 * IF (e.GT.0.995) THEN
135 * WRITE (6,9) (edot(i),i=1,3),(zt*ehat(i),i=1,3)
136 * 9 FORMAT (' edot zt*ehat ',1P,3E10.2,2X,3E10.2)
137 * END IF
138 *
139 * Include tidal and quadrupole terms for each star and also GR.
140  do i=1,3
141  edot(i) = edot(i) + zt*ehat(i) + (zq + zg)*qhat(i)
142 * hdot(i) = hdot(i) + zh*hhat(i)
143  enddo
144 * Save scalar value of eccentricity derivative for decision-making.
145  dedt = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
146 *
147 * IF (e.GT.0.995) THEN
148 * EDF = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
149 * WRITE (6,3) e, EKD, (EKD - EDF)/EKD, e/zt, e/zq, e/zg
150 * 3 FORMAT (' E EKD DED/ED TT TQ TGR ',F9.5,1P,5E9.1)
151 * END IF
152 * IF (IGR.LT.10.AND.zg.GT.zq*(cq(1) + cq(2)).AND.e.GT.0.99) THEN
153 * IF (e.GT.0.995) THEN
154 * ZD = ABS(EKD)
155 * ZD = MAX(ZD,zq*(cq(1)+cq(2)))
156 * if (e.GT.0.05.AND.zg.GT.ZD) THEN
157 * WRITE (6,4) e, a, a*(1.0-e), zg, ZD, 6.28/zg
158 * 4 FORMAT (' GR E A QP ZGR QUAD P_gr ',
159 * & F8.4,1P,5E10.2)
160 * IGR = IGR + 1
161 * END IF
162 *
163 * TAU = e/sqrt(edot(1)**2+edot(2)**2+edot(3)**2)
164 * WRITE (6,5) e, TAV, TAU, edot, zt, zq
165 * 5 FORMAT (' DERIV e TAV TAU ed zt zt ',F9.5,2F9.3,1P,5E10.2)
166  itime = itime + 1
167  IF (e.GT.0.980.AND.itime.EQ.-4) THEN
168  edt = zt
169  edq = zq
170  tf = e/edt
171  tq = e/edq
172  edtot = edot(1)*ehat(1)+edot(2)*ehat(2)+edot(3)*ehat(3)
173  IF (itime.EQ.4) WRITE (6,6) e0,tf,tq,edt,edq,edtot,ekd
174  6 FORMAT (' DERIV e TF TQ EDT EDQ ED EK ',f9.5,1p,6e10.2)
175  END IF
176  IF (itime.GE.4) itime = 0
177 
178 * Copy derivatives for Runge-Kutta integrator.
179  do i=1,3
180  du(i)=edot(i)
181  du(i+3)=hdot(i)
182  enddo
183  10 return
184  end