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
|
#!/usr/bin/env ruby
# This solves the Saha's equation, to investigate
# 'recombination' of hot cosmic plasma.
require("gsl")
include GSL
include GSL::CONST::CGSM
include Math
KB = BOLTZMANN
HBAR = PLANCKS_CONSTANT_HBAR
class Recombination
def initialize(t0, bdensity, eion)
@t = 2500
@tCMB = t0 # Present temperature of CMB [K]
@n0 = bdensity/MASS_PROTON # Barion number density [cm-3]
@eion = eion*ELECTRON_VOLT # Hydrogen ionization potential [eV]
@saha_equation = Function.alloc { |x| # x: fractional ionization
n = @n0*pow_3(@t/@tCMB)
tmp1 = pow(MASS_ELECTRON*KB*@t/2.0/PI/HBAR/HBAR, 1.5)
tmp2 = exp(-@eion/KB/@t)
x*x - (1.0 - x)*tmp1*tmp2/n
}
end
def zplus1()
@t/@tCMB
end
attr_accessor :t
attr_reader :saha_equation
end
###########################
## Cosmological parameter
TCMB = 2.725 # Temperature of cosmic microwave background at present [K]
BDENSITY = 1e-29 # Barion density [g cm-3]
EIONIZE = 13.6 # Hydrogen ionization potential [eV]
unv = Recombination.new(TCMB, BDENSITY, EIONIZE)
# FSolver
solver = Root::FSolver.alloc(Root::FSolver::BRENT)
TMIN = 2500
TMAX = 6500
begin
file = File.open("recombination.dat", "w")
t = TMIN
while t < TMAX
unv.t = t # Temperature
x, iter, status = solver.solve(unv.saha_equation, [1e-8, 1.0], [0, 1e-4])
file.printf("%e %e %e\n", unv.t, unv.zplus1, x)
t *= 1.01
end
ensure
file.close
end
system("gnuplot -persist recombination.gp")
File.delete("recombination.dat")
#puts("Try \"gnuplot -persist recombination.gp\"")
|
