[mesa-users] Change reaction rate of the triple alpha process
Byeong-Chan Park
pbc1918 at gmail.com
Mon Aug 1 06:15:09 EDT 2016
Dears.
I need to change the rate to investigate energy generation process in a
star. Hence I need to find the code that means nuclear reaction rates as a
function of temperature. I checked 'ratelib.f90', 'rates_def_mic.f90', and
'rates_def.f90' in the mesa/rates/privates and mesa/rates/public directory.
However I can't distinguish what code says the reaction rate formular of
the triple alpha process in these 3 files. I think that the next in the
mesa/rates/private/ratelib.f90 includes reaction formular;
! r3a, triple alpha
subroutine rate_tripalf_jina(tf, temp, fr, rr)
type (T_Factors), pointer :: tf
real(dp), intent(in) :: temp
real(dp), intent(out) :: fr, rr
real(dp) :: fr1, rr1
include 'formats.dek'
! he4 he4 he4 c12 fy05r 7.27500d+00
call jina_reaclib_3_1(ihe4, ihe4, ihe4, ic12, tf, fr, rr,
'rate_tripalf_jina')
return
call rate_tripalf_reaclib(tf, temp, fr1, rr1)
write(*,1) 'fr', fr
write(*,1) 'fr1', fr1
write(*,1) 'rr', rr
write(*,1) 'rr1', rr1
write(*,*)
stop 'rate_tripalf_jina'
end subroutine rate_tripalf_jina
subroutine rate_tripalf_reaclib(tf, temp, fr, rr)
! HE4(2A,G)C12 reaclib JINA - Fynbo et al. 2005 Nature 433,
136-139
type (T_Factors), pointer :: tf
real(dp), intent(in) :: temp
real(dp), intent(out) :: fr, rr
real(dp) :: fr1, fr2, fr3, rev
rr = 0
if (temp < 1d7) then
fr = 0; return
end if
call do_reaclib(tf, &
-9.710520d-01, 0.000000d+00, -3.706000d+01,
2.934930d+01, &
-1.155070d+02, -1.000000d+01,
-1.333330d+00, &
fr1)
call do_reaclib(tf, &
-2.435050d+01, -4.126560d+00, -1.349000d+01,
2.142590d+01, &
-1.347690d+00, 8.798160d-02,
-1.316530d+01, &
fr2)
call do_reaclib(tf, &
-1.178840d+01, -1.024460d+00, -2.357000d+01,
2.048860d+01, &
-1.298820d+01, -2.000000d+01,
-2.166670d+00, &
fr3)
fr = fr1 + fr2 + fr3
! use the fxt reverse rate term
rev = 2.00d+20*(tf% t93)*exp_cr(-84.424*(tf% t9i))
rr = fr * rev
end subroutine rate_tripalf_reaclib
subroutine rate_tripalf_nacre(tf, temp, fr, rr)
type (T_Factors), pointer :: tf
real(dp), intent(in) :: temp
real(dp), intent(out) :: fr, rr
real(dp) :: r2abe, rbeac, bb, term, rev
! he4(a, g)be8
! a0 T9i23 exp_cr(-a1 T9i13 - (T9*a2)^2)
! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
! + c0 T9i32 exp_cr(-c1/T9)
! + d0 T9i32 exp_cr(-d1/T9)
! + e0 T9^e1 exp_cr(-e2/T9)
if (tf% t9 < lowT9_cutoff) then
fr = 0; rr = 0; return
end if
call rnacre(tf, &
2.43d9, 13.490d0, 1d0/0.15d0, & ! a0, a1, a2
74.5d0, 0d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
6.09d5, 1.054d0, & ! c0, c1
0d0, 0d0, & ! d0, d1
0d0, 0d0, 0d0, & ! e0, e1, e2
r2abe)
! be8(a, g)c12
! a0 T9i23 exp_cr(-a1 T9i13 - (T9*a2)^2)
! * (1 + b0 T9 + b1 T92 + b2 T93 + b3 T94 + b4 T95)
! + c0 T9i32 exp_cr(-c1/T9)
! + d0 T9i32 exp_cr(-d1/T9)
! + e0 T9^e1 exp_cr(-e2/T9)
call rnacre(tf, &
2.76d7, 23.570d0, 1d0/0.4d0, & ! a0, a1, a2
5.47d0, 326d0, 0d0, 0d0, 0d0, & ! b0, b1, b2, b3, b4
130.7d0, 3.338d0, & ! c0, c1
2.51d4, 20.307d0, & ! d0, d1
0d0, 0d0, 0d0, & ! e0, e1, e2
rbeac)
if (tf% T9 <= 0.03d0) then
bb = 3.07d-16*(1 - 29.1d0*(tf% T9) + 1308d0*(tf% T92))
if (bb < 0) then
bb = 0
end if
else
bb = 3.44d-16*(1 + 0.0158d0*pow_cr(tf% T9,-0.65d0))
end if
term = r2abe * rbeac * bb
call rnacre_rev(tf, & ! a0 T932 exp_cr(-a1/T9)
2.003d20, 84.415d0, & ! a0, a1
rev)
fr = term
rr = rev * term
end subroutine rate_tripalf_nacre
subroutine rate_tripalf_fxt(tf, temp, fr, rr)
type (T_Factors), pointer :: tf
real(dp), intent(in) :: temp
real(dp), intent(out) :: fr, rr
real(dp) term, dtermdt, rev, drevdt, r2abe, dr2abedt, rbeac, &
drbeacdt, aa, daa, bb, dbb, cc, dcc, dd, ddd, ee,
dee, &
ff, dff, xx, dxx, yy, dyy, zz, dzz, uu, vv, f1, df1,
rc28, &
q1, q2
parameter (rc28 = 0.1d0, &
q1 = 1.0d0/0.009604d0, &
q2 = 1.0d0/0.055225d0)
if (tf% t9 < lowT9_cutoff) then
fr = 0; rr = 0; return
end if
! this is a(a, g)be8
aa = 7.40d+05 * (tf% t9i32) * exp_cr(-1.0663*(tf% t9i))
bb = 4.164d+09 * (tf% t9i23) * exp_cr(-13.49*(tf% t9i13) - (tf%
t92)*q1)
cc = 1.0d0 + 0.031*(tf% t913) + 8.009*(tf% t923) + 1.732*(tf%
t9) &
+ 49.883*(tf% t943) + 27.426*(tf% t953)
r2abe = aa + bb * cc
! this is be8(a, g)c12
dd = 130.0d0 * (tf% t9i32) * exp_cr(-3.3364*(tf% t9i))
ee = 2.510d+07 * (tf% t9i23) * exp_cr(-23.57*(tf% t9i13) - (tf%
t92)*q2)
ff = 1.0d0 + 0.018*(tf% t913) + 5.249*(tf% t923) + 0.650*(tf%
t9) + &
19.176*(tf% t943) + 6.034*(tf% t953)
rbeac = dd + ee * ff
! a factor
xx = rc28 * 1.35d-07 * (tf% t9i32) * exp_cr(-24.811*(tf% t9i))
! high temperature rate
if ((tf% t9).gt.0.08) then
term = 2.90d-16 * r2abe * rbeac + xx
! low temperature rate
else
uu = 0.8d0*exp_cr(-pow_cr(0.025*(tf% t9i),3.263d0))
yy = 0.2d0 + uu ! fixes a typo in Frank's original
vv = 4.0d0*exp_cr(-pow_cr((tf% t9)/0.025d0,9.227d0))
zz = 1.0d0 + vv
aa = 1.0d0/zz
f1 = 0.01 + yy * aa ! fixes a typo in Frank's original
term = 2.90d-16 * r2abe * rbeac * f1 + xx
end if
rev = 2.00d+20*(tf% t93)*exp_cr(-84.424d0*(tf% t9i))
! term = 1.2d0 * term
! dtermdt = 1.2d0 * term
fr = term
rr = rev * term
end subroutine rate_tripalf_fxt
but I don't know what code that I should change. In addition, it makes me
confused that form of the 'subroutine rate_tripalf_jina' is different with
the others.
I need your help.
Best regards
Park.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20160801/39cd321e/attachment.html>
More information about the Mesa-users
mailing list