Nbody6
newhut.f
Go to the documentation of this file.
1  subroutine hut2(spin10,spin20,spin1,spin2,nsteps,dtau)
2 *
3 * spin evolution of circular binary.
4 * ----------------------------------
5
6  implicit real*8 (a-h,o-z)
7  real*8 u(2),udot(2)
8
9  u(1)=spin10
10  u(2)=spin20
11
12  call deriv3(u,udot)
13
14 * include step reduction for large derivatives.
15  it = 0
16  1 IF (abs(udot(1))*dtau.GT.0.01*spin10) THEN
17  dtau = 0.5d0*dtau
18  nsteps = 2*nsteps
19  it = it + 1
20  IF (it.LT.5) go to 1
21  END IF
22
23  5 IF (abs(udot(2))*dtau.GT.0.01*spin20) THEN
24  dtau = 0.5d0*dtau
25  nsteps = 2*nsteps
26  it = it + 1
27  IF (it.LT.5) go to 5
28  END IF
29
30  do i=1,nsteps
31  call rk4c(dtau,u)
32  enddo
33
34  spin1=u(1)
35  spin2=u(2)
36
37  end
38
39  subroutine rk4c(dt,u)
40
41 * runge-kutta integrator.
42 * -----------------------
43
44 * author: rosemary mardling(3/98).
45
46  parameter(n=2)
47  implicit real*8 (a-h,o-z)
48  real*8 u0(n),ut(n),du(n),u(n)
49  real*8 a(4),b(4)
50
51  a(1)=dt/2.0d0
52  a(2)=a(1)
53  a(3)=dt
54  a(4)=0.0d0
55
56  b(1)=dt/6.0d0
57  b(2)=dt/3.0d0
58  b(3)=b(2)
59  b(4)=b(1)
60
61  do i=1,n
62  u0(i)=u(i)
63  ut(i)=0.0d0
64  enddo
65
66  do j=1,4
67  call deriv3(u,du)
68
69  do i=1,n
70  u(i)=u0(i)+a(j)*du(i)
71  ut(i)=ut(i)+b(j)*du(i)
72  enddo
73  enddo
74
75  do i=1,n
76  u(i)=u0(i)+ut(i)
77  enddo
78
79  end
80
81  subroutine deriv3(u,udot)
82
83  implicit real*8 (a-h,m,o-z)
84  real*8 u(2),udot(2)
85  common/spins/angmom0,rg2(2),m21,r21,semi0,c1,c2,c3,c4,c5,semi
86  SAVE ic
87  DATA ic /0/
88
89 * assume e=0.
90  spin1=u(1)
91  spin2=u(2)
92
93 * e2=e**2
94 * e4=e2**2
95 * e6=e4*e2
96 * fac=1-e2
97
98 * f2=1+7.5*e2+5.625*e4+0.3125*e6
99 * f3=1+3.75*e2+1.875*e4+0.078125*e6
100 * f4=1+1.5*e2+0.125*e4
101 * f5=1+3*e2+0.375*e4
102
103
104  semi=angmom0-rg2(1)*spin1-m21*r21**2*rg2(2)*spin2
105 * semi=(semi*(1+m21)/m21/semi0**2)**2/fac
106  semi=(semi*(1+m21)/m21/semi0**2)**2
107  oa = 1.0/semi
108
109 * udot(1)=(oa/fac)**6*c3*(oa**1.5*f2-fac**1.5*f5*spin1)
110 * udot(2)=(oa/fac)**6*c4*(oa**1.5*f2-fac**1.5*f5*spin2)
111  udot(1)=oa**6*c3*(oa**1.5-spin1) - c5*spin1
112  udot(2)=oa**6*c4*(oa**1.5-spin2)
113
114 * ic = ic + 1
115 * IF (ic.LT.10) THEN
116 * WRITE (6,1) ic, spin1, spin2, udot
117 * 1 FORMAT (' HUT DERIV # s1 s2 udot ',i8,1p,6e9.1)
118 * END IF
119
120  end