====== 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:
<code fortran>
      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
</code>

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'':
<code fortran>
! 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 
</code>
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:
<code fortran>
      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
</code>

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:
<code bash>
#!/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
</code>

This script expects an ''$inputfile'' as follows (made-up data for now!)
<code>
#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
</code>

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.

===== Downloads =====
The complete BBN code with our changes can be found {{::bigbang_deuteronrate.zip|here}}.

The rates and other input variables have been calculated outside the BBN code.