Nbody6
himod.f
Go to the documentation of this file.
1  SUBROUTINE himod(evec,hvec,m3,mt,EOvec,HOvec,ICALL,dt0,dt,e,r,v)
2 *
3 *
4 * Modification of hierarchical binary.
5 * ------------------------------------
6 *
7 * Author: Rosemary Mardling (3/98).
8
9  IMPLICIT REAL*8 (a-h,m,o-z)
10  common/rksave/ coeff,hohat(3),e0,a,hh,mb
11  REAL*8 evec(3),hvec(3),eovec(3),hovec(3),
12  & element(6),r(3),v(3),u(6),udot(6)
13 *
14 *
15 * Given inner Runge-Lenz vector {\bf e}=evec and {\bf h}=hvec.
16 * Calculate eccentricity e, magnitude of hvec and semi-major axis a.
17
18  mb=mt-m3
19  e0=sqrt(dot(evec,evec))
20  hh=sqrt(dot(hvec,hvec))
21  a=hh**2/mb/(1-e0**2)
22
23 * Given OUTER {\bf E}=Eout and {\bf H}=Hout.
24 * Calculate outer eccentricity and semi-major axis.
25
26  eout=sqrt(dot(eovec,eovec))
27  hhout=sqrt(dot(hovec,hovec))
28  aout=hhout**2/mt/(1-eout**2)
29  do i=1,3
30  hohat(i)=hovec(i)/hhout
31  enddo
32
33 * Calculate the semi-minor axis Bout and the semi-latus rectum SLRout.
34
35  bout=aout*sqrt(1-eout**2)
36  slrout=aout*(1-eout**2)
37  coeff=m3/(2*aout*bout*slrout)
38  dtsum = 0.0
39
40  IF (icall.LT.0) THEN
41  do i=1,3
42  u(i)=evec(i)
43  u(i+3)=hvec(i)
44  enddo
45  CALL deriv(u,udot,icall)
46  tau = e0/sqrt(udot(1)**2+udot(2)**2+udot(3)**2)
47  WRITE (6,5) e0,tau,a,eout,aout,(udot(k),k=1,3)
48  5 FORMAT (' DERIV e0 TAU a E1 A1 edot ',
49  & f8.4,1p,2e10.2,0p,f8.4,1p,4e10.2)
50  CALL flush(3)
51  icall = icall + 1
52  END IF
53
54 * Initialize integration variables.
55
56  10 do i=1,3
57  u(i)=evec(i)
58  u(i+3)=hvec(i)
59  enddo
60
61 * Perform integration.
62
63  call rkint(dt,u)
64
65  if(e0.ge.1.0) goto 20
66  eh=0.0
67  do i=1,3
68  evec(i)=u(i)
69  hvec(i)=u(i+3)
70  eh=eh+evec(i)*hvec(i)
71  enddo
72
73  e0=sqrt(dot(evec,evec))
74  hh=sqrt(dot(hvec,hvec))
75
76 * Apply symmetric re-orthogonalization (small eh; Seppo Mikkola 3/98).
77
78  eps=-eh/(e0**2+hh**2)
79  eps=0.0
80  do k=1,3
81  evec(k)=evec(k)+eps*hvec(k)
82  hvec(k)=hvec(k)+eps*evec(k)
83  end do
84
85  IF (dtsum + dt.GT.dt0) dt = dt0 - dtsum + 1.0d-10
86  dtsum = dtsum + dt
87  IF (dtsum.LT.dt0) go to 10
88
89 * Calculate new orbital parameters.
90  20 e2=dot(evec,evec)
91  h2=dot(hvec,hvec)
92  a=h2/mb/(1-e2)
93  e=sqrt(e2)
94  hh=sqrt(h2)
95
96 * Calculate updated vectors r and v.
97 * Note that the vectors evec and hvec give no information about
98 * the orbital phase. Thus there is a kind of degeneracy. Since
99 * we are updating to a later time using an approximation, we can't
100 * expect to know it anyway. Hence the phase of peri is set to zero.
101
102  call transform4(evec,hvec,mb,element)
103  call transform2(element,mb,r,v)
104
105  RETURN
106  END
107
108
109  subroutine cross(u,v,w)
110
111 * Vectorial cross product.
112 * ------------------------
113  real*8 u(3),v(3),w(3)
114
115  w(1) = u(2)*v(3) - u(3)*v(2)
116  w(2) = u(3)*v(1) - u(1)*v(3)
117  w(3) = u(1)*v(2) - u(2)*v(1)
118
119  end
120
121
122  real*8 function dot(u,v)
123
124 * Scalar dot product.
125 * ------------------
126  real*8 u(3),v(3)
127
128  dot=u(1)*v(1)+u(2)*v(2)+u(3)*v(3)
129
130  end
131
132
133  subroutine transform2(element,mb,r,v)
134 *
135 * Calculates vectors r and v given 6 orbital elements.
136 * ----------------------------------------------------
137 *
138  IMPLICIT REAL*8 (a-h,o-z)
139  real*8 r(3),v(3),element(6),mb
140  real*8 tempr(3),tempv(3)
141  real*8 inc,b(3,3)
142
143
144  a=element(1)
145  e=element(2)
146  inc=element(3)
147  w=element(4)
148  om=element(5)
149  phi=element(6)
150
151  cosp=cos(phi)
152  sinp=sin(phi)
153  cosi=cos(inc)
154  sini=sin(inc)
155  cosw=cos(w)
156  sinw=sin(w)
157  cosom=cos(om)
158  sinom=sin(om)
159
160  b(1,1) = cosw*cosom - sinw*cosi*sinom
161  b(1,2) = cosw*sinom + sinw*cosi*cosom
162  b(1,3) = sinw*sini
163  b(2,1) = -sinw*cosom - cosw*cosi*sinom
164  b(2,2) = -sinw*sinom + cosw*cosi*cosom
165  b(2,3) = cosw*sini
166  b(3,1) = sini*sinom
167  b(3,2) = -sini*cosom
168  b(3,3) = cosi
169
170  h=sqrt(mb*a*(1-e**2))
171  rr=a*(1-e**2)/(1+e*cosp)
172  rd=e*h*sinp/(a*(1-e**2))
173  phid=h/rr**2
174
175  r(1)=rr*cosp
176  r(2)=rr*sinp
177  r(3)=0.0
178
179  v(1)=rd*cosp-rr*phid*sinp
180  v(2)=rd*sinp+rr*phid*cosp
181  v(3)=0.0
182
183  do i=1,3
184  sum1=0.0
185  sum2=0.0
186  do j=1,3
187  sum1=sum1 + b(j,i)*r(j)
188  sum2=sum2 + b(j,i)*v(j)
189  enddo
190  tempr(i)=sum1
191  tempv(i)=sum2
192  enddo
193  do i=1,3
194  r(i)=tempr(i)
195  v(i)=tempv(i)
196  enddo
197
198  end
199
200
201  subroutine transform4(e,h,mb,element)
202 *
203 * Calculates 5 orbital elements given vectors evec and h and mass mb.
204 * -------------------------------------------------------------------
205 *
206  IMPLICIT REAL*8 (a-h,o-z)
207  real*8 element(6)
208  real*8 h(3),e(3),n(3),hn(3),he(3)
209  real*8 eh(3)
210  real*8 ii(3),jj(3),kk(3),nn,mb
211  data ii/1.d0,0.d0,0.d0/
212  data jj/0.d0,1.d0,0.d0/
213  data kk/0.d0,0.d0,1.d0/
214
215
216  hh=sqrt(dot(h,h))
217
218  do i=1,3
219  h(i)=h(i)/hh
220  enddo
221
222  call cross(kk,h,n)
223
224  nn=sqrt(dot(n,n))
225  ecc=sqrt(dot(e,e))
226  a=hh**2/mb/(1-ecc**2)
227
228  do i=1,3
229  n(i)=n(i)/nn
230  eh(i)=e(i)/ecc
231  enddo
232
233  call cross(h,n,hn)
234  call cross(h,eh,he)
235
236  cosom=dot(ii,n)
237  sinom=dot(jj,n)
238  cosi=dot(kk,h)
239  sini=nn
240  cosw=dot(n,eh)
241  sinw=dot(eh,hn)
242
243  element(1)=a
244  element(2)=ecc
245  element(3)=atan2(sini,cosi)
246  element(4)=atan2(sinw,cosw)
247  element(5)=atan2(sinom,cosom)
248  element(6)=0.0
249 * Pericentre phase (can be specified).
250
251 * ZI = 360.0*ELEMENT(3)/(8.0*ATAN(1.0D0))
252 * WRITE (6,5) ECC,ELEMENT(3),ELEMENT(4),ELEMENT(5),ZI
253 * 5 FORMAT (' TRANSFORM e i w Om IN ',F8.4,F9.3,3F9.3)
254  end