Nuclear Talent
Course 7
NT4A

group2_code

Code documentation for group $^{\rm BSM}_{\rm\ \ BBN}$

Wrapper script to batch-modify a reaction rate

First of all, we parametrize the reaction rate and the Q-value (for the back-reaction) in the corresponding subroutine in the main code file:

      subroutine rate_png(temp,den,fr,dfrdt,dfrdd,rr,drrdt,drrdd)
include 'implno.dek'
include 'tfactors.dek'
include 'helper_deuteronrate.dek'

! declare the pass
double precision temp,den,fr,dfrdt,dfrdd,rr,drrdt,drrdd

! locals
double precision term,dtermdt,rev,drevdt,aa,daa

! p(n,g)d
! smith,kawano,malany 1992

aa      = 1.0d0 - npgd_a1*t912 + npgd_a2*t9 &
- npgd_a3*t932 + npgd_a4*t92 &
- npgd_a5*t952

daa     =  -0.5d0*npgd_a1*t9i12 + npgd_a2 &
- 1.5d0*npgd_a3*t912 + 2.0d0*npgd_a4*t9 &
- 2.5d0*npgd_a5*t932

term    = npgd_a0 * aa
dtermdt = npgd_a0 * daa

! rates
fr    = den * term
dfrdt = den * dtermdt * 1.0d-9
dfrdd = term

rev      = 4.71d+09 * t932 * exp(-npgd_Eb*t9i)
drevdt   = rev*(1.5d0*t9i + npgd_Eb*t9i2)

rr    = rev * term
drrdt = (drevdt * term + rev * dtermdt) * 1.0d-9
drrdd = 0.0d0

return
end

These parameters should be read in in another subroutine, so they have to be global variables. Hence we create a helper file helper_deuteronrate.dek:

! helper file to create global variables for the n(p,g)d reaction

real*8 :: npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7
common /helpdeut/ npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7 

This helper file needs of course to be added to all subroutines working with these parameters.

We then read in these variables in the net_input subroutine after querying the cosmological variables:

      write(6,02) 'enter a0 to a7, mn and mp => '
read(5,*) npgd_Eb, npgd_a0, npgd_a1, npgd_a2, npgd_a3, npgd_a4, npgd_a5, npgd_a6, npgd_a7, mn, mp

We then create a wrapper script to read in the parameters from a file, run the BBN code and save the final composition of each run to a new file:

#!/bin/bash

# wrapper script to batch-feed mn, mp and the p(n,g)d rate coefficients into bbn code

### config section ###
prefix="test_"
endtime=1e6
inputfile="coefficients.in" #input file containing a column used as label as well as a0-a7, mn and mp for the p(n,g)d reaction rate
variablelabel="Lambda_QCD" #just a label used in the final.dat file

rm $prefix* n=0 while read line; do if [${line:0:1} != "#" ]; then #skip comment lines
prefixn=$prefix$n"_"

var=echo $line | awk '{ print$1; }'
params=echo $line | cut -f2- -d ' ' echo -e "\n$params\n$endtime\n$prefixn\n1\n" | ./bigbang_deuteron.out

#copy header for first file
if [ $n -eq 0 ]; then head -n 2$prefixn'0002_z1.dat' > $prefix"final.dat" echo " * variable varied:$variablelabel" >> $prefix"final.dat" awk '{$1=$2="";print;}'$prefixn'0003_z1.dat'  | paste $prefixn'0002_z1.dat' - | sed '3q;d' >>$prefix"final.dat"
fi

result1=tail -n 1 $prefixn'0002_z1.dat' result2=tail -n 1$prefixn'0003_z1.dat' | awk '{$1=$2="";print;}' # | cut -f3-
echo $var" "$result1" "$result2 >>$prefix"final.dat"
n=$((n+1)) fi done <$inputfile

This script expects an $inputfile as follows (made-up data for now!) #Lambda Eb a0 a1 a2 a3 a4 a5 a6 a7 mn mp 217 25.82 4.742e4 0.8504 0.4895 0.09623 8.471e-3 2.80e-4 0 0 1.67492721184d-24 1.67262163783d-24 217 27 4.742e4 0.8504 0.4895 0.09623 8.471e-3 2.80e-4 0 0 1.67492721184d-24 1.67262163783d-24 220 27 7.342e4 0.8504 0.4895 0.09623 8.471e-3 2.80e-4 0 0 1.67492721184d-24 1.67492721184d-24 Now, “all” which remains is to figure out how the binding energy, the a's and the nucleon masses depend on$\Lambda_{QCD}\$ and to calculate the above table.