[mesa-users] Change reaction rate of the triple alpha process
Francis Timmes
fxt44 at mac.com
Mon Aug 1 11:39:03 EDT 2016
changing code is probably not what you want;
rate multipliers can be done via the inlist:
http://mesa.sourceforge.net/star_job_defaults.html#num_special_rate_factors
the overall scaling factors
http://mesa.sourceforge.net/controls_defaults.html#eps_nuc_factor
http://mesa.sourceforge.net/controls_defaults.html#dxdt_nuc_factor
may also help you reach your goal.
fxt
> On Aug 1, 2016, at 3:15 AM, Byeong-Chan Park <pbc1918 at gmail.com> wrote:
>
> 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.
> ------------------------------------------------------------------------------
> _______________________________________________
> mesa-users mailing list
> mesa-users at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/mesa-users
More information about the Mesa-users
mailing list