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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
#!/usr/bin/env ruby
# This example solves Lane-Emden equation to calculate white dwarf
# structure, using Chandrasekhar's model, as a polytrope gas sphere
# of self-gravitating body supported by eletron degeneration pressure.
require("gsl")
include Math
include GSL::CONST::CGSM
include GSL::Odeiv
# Degenerated electron/neutron gas
# Reference: "Black holes, white dwarfs, and neutron stars"
# S.L. Shapiro & S.A. Teukolsky, Willey, 1983
module Degenerate
hbar = PLANCKS_CONSTANT_HBAR # Planck's constant
me = MASS_ELECTRON # Electron mass
mn = MASS_NEUTRON # Neutron mass
mu = UNIFIED_ATOMIC_MASS
c = SPEED_OF_LIGHT
ELambda = hbar/me/c # Compton length of electron
NLambda = hbar/mn/c
MeC2 = me*c*c # Electron rest mass energy
MnC2 = mn*c*c
Factor_xe = GSL::pow(3.0*PI*PI/mu, 1.0/3.0)*ELambda
Factor_xn = GSL::pow(3.0*PI*PI/mu, 1.0/3.0)*NLambda
# See Shapiro&Teukolsky Chapter 2
def phi(x)
tmp = sqrt(1.0 + x*x)
(x*tmp*(2.0*x*x/3.0 - 1.0) + log(x + tmp))/8/PI/PI
end
def chi(x)
tmp = sqrt(1.0 + x*x)
(x*tmp*(1.0 + 2*x*x) - log(x + tmp))/8/PI/PI
end
def xe_rho_mue(rho, mue)
Factor_xe*GSL::pow(rho/mue, 1.0/3.0)
end
def xn_rho(rho)
Factor_xn*GSL::pow(rho, 1.0/3.0)
end
end
# Polytrope gas sphere
module Polytrope
# Lane-Emden equation
# n: polytrope index
EmdenEq = Proc.new { |x, y, dydx, n|
dydx[0] = y[1]
dydx[1] = -GSL::pow(y[0], n) - 2.0/x*y[1]
}
def emden_xi(n)
dim = 2
y = GSL::Vector[1.0, 0.0]
dydx = GSL::Vector.alloc(dim)
solver = Solver.alloc(Step::RKF45, [1e-6, 0], EmdenEq, dim)
solver.set_params(n)
solver.reset
vx = GSL::Vector.alloc(10000)
vy = GSL::Vector.alloc(10000)
vdy = GSL::Vector.alloc(10000)
x = 0.0001
xend = 10.0
h = 1e-6
file = File.open("polytrope.dat", "w")
i = 0
while x < xend
x, h, status = solver.apply(x, xend, h, y)
break if GSL::isnan?(y[0])
break if GSL::isnan?(y[1])
vx[i] = x
vy[i] = y[0]
vdy[i] = y[1]
file.printf("%e %e %e %e\n", x, y[0], y[1], GSL::pow_3(y[0]))
i += 1
break if status != GSL::SUCCESS
break if y[0] <= -0.1
end
file.close
# p vx.size
# p i
vx2 = vx.subvector(0, i)
vy2 = vy.subvector(0, i)
vdy2 = vdy.subvector(0, i)
spline = GSL::Spline.alloc(GSL::Interp::AKIMA, i)
spline.init(vy2.reverse, vx2.reverse)
# Find star surface:
# Star sufrace is defined as the zero point of density structure function
x1 = spline.eval(0.0)
spline.init(vx2, vdy2)
yx2 = spline.eval(x1).abs
return [x1, yx2*x1*x1]
end
end
# Chandrasekhar white dwarf model:
# * Polytrope gas sphere
# * Support its self gravity by electron degeneration pressure
class WhiteDwarf
include Degenerate
include Polytrope
G = GRAVITATIONAL_CONSTANT
@@x1 = nil
@@x12 = nil
def initialize(n, rhoc, mue)
@n = n
@gam = 1.0 + 1.0/n
@rhoc = rhoc
@mue = mue
x = xe_rho_mue(@rhoc, mue)
phix = phi(x)
@K = MeC2*phix/GSL::pow_3(ELambda)/GSL::pow(@rhoc, @gam)
@a = sqrt((@n + 1)*@K*GSL::pow(@rhoc, 1.0/@n - 1.0)/4/PI/G)
if !@@x1
@@x1, @@x12 = emden_xi(@n)
end
@Pc = MeC2*phix/GSL::pow_3(ELambda)*phix # Electron Fermi pressure
@radius = @a*@@x1 # White dwarf radius
@mass = 4.0*PI*GSL::pow_3(@a)*@rhoc*@@x12 # White dwarf mass
end
attr_accessor :radius
attr_accessor :mass
attr_accessor :rhoc # Central density
attr_accessor :Pc # Central pressure
attr_accessor :gam # Adiabatic index
attr_accessor :n # Polytrope index
attr_accessor :mue # Mean eletron number per barion
attr_accessor :K # Normalization constant
attr_accessor :a
end
file = File.open("whitedwarf.dat", "w")
rho = 1e6 # central density
while rho < 1e11
wd = WhiteDwarf.new(3, rho, 2)
file.printf("%e %e %e %e\n", rho, wd.Pc, wd.radius/1000/100, wd.mass/SOLAR_MASS)
rho *= 1.2
end
file.close
system("gnuplot -persist whitedwarf.gp")
File.delete("polytrope.dat")
File.delete("whitedwarf.dat")
|
