-
Notifications
You must be signed in to change notification settings - Fork 115
Expand file tree
/
Copy pathplot_accretion_from_windy_star.py
More file actions
108 lines (91 loc) · 3.06 KB
/
Copy pathplot_accretion_from_windy_star.py
File metadata and controls
108 lines (91 loc) · 3.06 KB
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
"""
Visualization for simple N-body integration.
Reads particle set from file (nbody.hdf5) and prints subsequent frames.
"""
import sys
import numpy
from matplotlib import pyplot
from amuse.lab import *
def Bondi_Hoyle_Littleton_accretion_rate(Mb, vs, a0, Mdot_donor):
a = (3.e-7 | units.MSun/units.yr)
b = (Mb/(1|units.MSun))**2
c = ((10|units.kms)/vs)**4
d = ((100|units.AU)/a0)**2
e = (Mdot_donor/(1.e-4 | units.MSun/units.yr))
print(a, b, c, d, e)
Mdot_BHL = a*b*c*d*e
#Mdot_BHL = (3.e-7 | units.MSun/units.yr) * (Mb/(1|units.MSun))**2 ((10|units.kms)/vs)**4 * ((100|units.AU)/a0)**2 * (Mdot_donor/(1.e-4 | units.MSun/units.yr))
return Mdot_BHL
def read_accretion_rate(filename):
t = [] | units.yr
m = [] | units.MSun
n = []
for line in open(filename):
if "N accreted:" in line:
l = line.split()
t.append(float(l[2])|units.Myr)
n.append(float(l[4]))
m.append(float(l[6])|units.MSun)
return t, n, m
def plot_accretion_from_wind(filename):
t, dn, dm = read_accretion_rate(filename)
m = numpy.cumsum(dm)
#MMoon = 3.69145063653e-08 | units.MSun
#m /= MMoon
m /= (1.e-9|units.MSun)
pyplot.plot(t.value_in(units.yr), m)
# pyplot.show()
def v_terminal_teff(temperature):
print(numpy.log10(temperature.value_in(units.K))-3.)
t4=(numpy.log10(temperature.value_in(units.K))-4.).clip(0.,1.)
print(t4)
return (30 | units.km/units.s) + ((4000 | units.km/units.s)*t4)
def main():
"""
Mb = (2. + 1.) | units.MSun
Mb = 2. | units.MSun
a0 = 10. | units.AU
# vs = 35. | units.kms
vs = 30. | units.kms
Mdot_donor = 0.11 |units.MSun/units.Myr
Mdot = Bondi_Hoyle_Littleton_accretion_rate(Mb, vs, a0, Mdot_donor)
print Mdot.in_(units.MEarth/units.yr)
t = numpy.arange(0, 3, 0.1) | units.yr
mdot = Mdot*t
t += 4 | units.yr
"""
Mb = (2. + 1.) | units.MSun
Mb = 2. | units.MSun
a0 = 10. | units.AU
#vs = 30. | units.kms
# vs = 18. | units.kms
vs = 17.26 | units.kms
vorb = numpy.sqrt(constants.G*Mb/a0)
print("vorb=", vorb.in_(units.kms))
k = vorb/vs
m1 = 1.924785833858|units.MSun
m2 = 1.|units.MSun
mu = m2/(m1+m2)
Mdot_donor = 0.11 |units.MSun/units.Myr
cvw = 0.
print("k and mu:", k, mu)
Mdot = Mdot_donor * mu**2 * k**4/(1 + k**2 + cvw**2)**(3./2.)
# print "Mdot:", Mdot.in_(units.MEarth/units.yr)
print("Mdot:", Mdot.in_(units.MSun/units.Myr))
t = numpy.arange(0, 3, 0.1) | units.yr
mdot = Mdot*t
t += 4.2 | units.yr
pyplot.figure()
pyplot.xlabel('t [yr]')
pyplot.ylabel('M [$10^{-9}$M$_{\odot}$]')
filename = "hydro_give_or_take.data"
print(t, mdot.value_in(units.MSun))
mdot /= (1.e-9|units.MSun)
pyplot.plot(t.value_in(units.yr), mdot)
plot_accretion_from_wind(filename)
filename = "hydro_give_or_take_gravity_NoG.data"
plot_accretion_from_wind(filename)
pyplot.savefig("hydro_accretion_from_windy_star")
# pyplot.show()
if __name__ in ('__main__', '__plot__'):
main()