-
Notifications
You must be signed in to change notification settings - Fork 26
/
kick.f
139 lines (139 loc) · 4.17 KB
/
kick.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
***
SUBROUTINE kick(kw,m1,m1n,m2,ecc,sep,jorb,vs)
implicit none
*
integer kw,k
INTEGER idum
COMMON /VALUE3/ idum
INTEGER idum2,iy,ir(32)
COMMON /RAND3/ idum2,iy,ir
integer bhflag
real*8 m1,m2,m1n,ecc,sep,jorb,ecc2
real*8 pi,twopi,gmrkm,yearsc,rsunkm
parameter(yearsc=3.1557d+07,rsunkm=6.96d+05)
real*8 mm,em,dif,der,del,r
real*8 u1,u2,vk,v(4),s,theta,phi
real*8 sphi,cphi,stheta,ctheta,salpha,calpha
real*8 vr,vr2,vk2,vn2,hn2
real*8 mu,cmu,vs(3),v1,v2,mx1,mx2
real*8 sigma
COMMON /VALUE4/ sigma,bhflag
real ran3,xx
external ran3
*
do k = 1,3
vs(k) = 0.d0
enddo
* if(kw.eq.14.and.bhflag.eq.0) goto 95
*
pi = ACOS(-1.d0)
twopi = 2.d0*pi
* Conversion factor to ensure velocities are in km/s using mass and
* radius in solar units.
gmrkm = 1.906125d+05
*
* Find the initial separation by randomly choosing a mean anomaly.
if(sep.gt.0.d0.and.ecc.ge.0.d0)then
xx = RAN3(idum)
mm = xx*twopi
em = mm
2 dif = em - ecc*SIN(em) - mm
if(ABS(dif/mm).le.1.0d-04) goto 3
der = 1.d0 - ecc*COS(em)
del = dif/der
em = em - del
goto 2
3 continue
r = sep*(1.d0 - ecc*COS(em))
*
* Find the initial relative velocity vector.
salpha = SQRT((sep*sep*(1.d0-ecc*ecc))/(r*(2.d0*sep-r)))
calpha = (-1.d0*ecc*SIN(em))/SQRT(1.d0-ecc*ecc*COS(em)*COS(em))
vr2 = gmrkm*(m1+m2)*(2.d0/r - 1.d0/sep)
vr = SQRT(vr2)
else
vr = 0.d0
vr2 = 0.d0
salpha = 0.d0
calpha = 0.d0
endif
*
* Generate Kick Velocity using Maxwellian Distribution (Phinney 1992).
* Use Henon's method for pairwise components (Douglas Heggie 22/5/97).
do 20 k = 1,2
u1 = RAN3(idum)
u2 = RAN3(idum)
* Generate two velocities from polar coordinates S & THETA.
s = sigma*SQRT(-2.d0*LOG(1.d0 - u1))
theta = twopi*u2
v(2*k-1) = s*COS(theta)
v(2*k) = s*SIN(theta)
20 continue
vk2 = v(1)**2 + v(2)**2 + v(3)**2
vk = SQRT(vk2)
if((kw.eq.14.and.bhflag.eq.0).or.kw.lt.0)then
vk2 = 0.d0
vk = 0.d0
if(kw.lt.0) kw = 13
endif
sphi = -1.d0 + 2.d0*u1
phi = ASIN(sphi)
cphi = COS(phi)
stheta = SIN(theta)
ctheta = COS(theta)
* WRITE(66,*)' KICK VK PHI THETA ',vk,phi,theta
if(sep.le.0.d0.or.ecc.lt.0.d0) goto 90
*
* Determine the magnitude of the new relative velocity.
vn2 = vk2+vr2-2.d0*vk*vr*(ctheta*cphi*salpha-stheta*cphi*calpha)
* Calculate the new semi-major axis.
sep = 2.d0/r - vn2/(gmrkm*(m1n+m2))
sep = 1.d0/sep
* if(sep.le.0.d0)then
* ecc = 1.1d0
* goto 90
* endif
* Determine the magnitude of the cross product of the separation vector
* and the new relative velocity.
v1 = vk2*sphi*sphi
v2 = (vk*ctheta*cphi-vr*salpha)**2
hn2 = r*r*(v1 + v2)
* Calculate the new eccentricity.
ecc2 = 1.d0 - hn2/(gmrkm*sep*(m1n+m2))
ecc2 = MAX(ecc2,0.d0)
ecc = SQRT(ecc2)
* Calculate the new orbital angular momentum taking care to convert
* hn to units of Rsun^2/yr.
jorb = (m1n*m2/(m1n+m2))*SQRT(hn2)*(yearsc/rsunkm)
* Determine the angle between the new and old orbital angular
* momentum vectors.
cmu = (vr*salpha-vk*ctheta*cphi)/SQRT(v1 + v2)
mu = ACOS(cmu)
* Calculate the components of the velocity of the new centre-of-mass.
90 continue
if(ecc.le.1.0)then
* Calculate the components of the velocity of the new centre-of-mass.
mx1 = vk*m1n/(m1n+m2)
mx2 = vr*(m1-m1n)*m2/((m1n+m2)*(m1+m2))
vs(1) = mx1*ctheta*cphi + mx2*salpha
vs(2) = mx1*stheta*cphi + mx2*calpha
vs(3) = mx1*sphi
else
* Calculate the relative hyperbolic velocity at infinity (simple method).
sep = r/(ecc-1.d0)
* cmu = SQRT(ecc-1.d0)
* mu = ATAN(cmu)
mu = ACOS(1.d0/ecc)
vr2 = gmrkm*(m1n+m2)/sep
vr = SQRT(vr2)
vs(1) = vr*SIN(mu)
vs(2) = vr*COS(mu)
vs(3) = 0.d0
ecc = MIN(ecc,99.99d0)
endif
*
95 continue
*
RETURN
END
***