Program XrayTrans685

Use Kinds, Only: dp
Use Kinds, Only: id
Use FileIO_Utilities, Only: Delta_Time
Implicit None
Integer(id) :: c

Call SYSTEM_CLOCK(COUNT=c)
Call Transport_Photons( '1e_5__' , 1000000 , hv0_pdf_Line_5 , (/1.E1_dp,1.E2_dp,1.E3_dp,1.E4_dp,1.E5_dp/) )
Write(*,'(ES12.3E2,A)') Real(60_id*1000000_id,dp)/Delta_Time(clock_then=c),' histories per minute'

Contains

Subroutine Transport_Photons(p_name,n_p,hv0_pdf,Rd)
    Use Kinds, Only: dp
    Use Global, Only: default_seed
    Use Global, Only: speed_of_light_SI
    Use Global, Only: m0c2
    Use Global, Only: Z_hat,Y_hat
    Use Global, Only: pi
    Use FileIO_Utilities, Only: creturn
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Character, Intent(In) :: p_name  !problem name, will start output filenames
    Integer, Intent(In) :: n_p  !number of particles to be run
    Interface
        Function hv0_pdf(r) Result(hv0)
            Use Kinds, Only: dp
            Use PRNGs, Only: ACORN_Type
            Real(dp) :: hv0
            Type(ACORN_Type), Intent(InOut) :: r
        End Function hv0_pdf
    End Interface
    Real(dp), Intent(In) :: Rd(:)  ![cm] radii of detector shells

    Real(dp), Parameter :: rho0 = 1225.0002517923321962E-6_dp  ![g/cm^3] atmospheric density at sea level
    Real(dp), Parameter :: speed_c = speed_of_light_SI * 100._dp  ![cm/s] speed of light
    Integer, Parameter :: photoelectric_absorption = 1
    Integer, Parameter :: incoherent_scatter = 2
    Integer, Parameter :: coherent_scatter = 3
    Integer, Parameter :: pair_production = 4

    Type(ACORN_Type) :: RNG
    Integer :: n_sig  !number of lines in the MAC list
    Real(dp), Allocatable :: hv_sigs(:)  ![keV] energy list for the MAC list, stored as logarithms
    Real(dp), Allocatable :: sigs(:,:)  ![cm^2/g] MAC lists, stored as logarithms
    Integer, Parameter :: n_E = 80
    Real(dp) :: E_list(0:n_E)
    Integer :: n_d
    Integer :: p_list_size
    Real(dp), Allocatable :: p_list(:,:,:)
    Integer, Allocatable :: p_listed(:)
    Real(dp), Allocatable :: f(:,:)
    Real(dp), Allocatable :: f_list(:,:,:)
    Integer :: p  !counter for particles
    Real(dp) :: r(1:3)  ![cm] position of the current photon
    Real(dp) :: omega_hat(1:3)  ![1 cm/s] direction of travel of the current photon
    Real(dp) :: hv  ![keV] energy of the current photon
    Real(dp) :: t  ![s] total flight time of the current photon
    Real(dp) :: w  !weight of the current photon
    Real(dp) :: r0  ![cm] magnitude of the photon position vector, before the most recent movement
    Real(dp) :: zeta0  !cosine of the zenith angle for the direction of travel for the current photon
    Real(dp) :: tau  ![mfp] optical thickness to the next interaction
    Real(dp) :: sig(1:4)  ![cm^2/g] mass attenuation coefficients for absorption, incoherent scat, coherent scat, and pair prod
    Real(dp) :: sigT  ![cm^2/g] total mass attenuation coefficient
    Real(dp) :: d  ![cm] geometric distance to the next interaction
    Real(dp) :: r1
    Real(dp) :: zeta
    Real(dp) :: Rca
    Real(dp) :: dca
    Real(dp) :: dz
    Real(dp) :: xi  !a scaled random number
    Integer :: interaction  !index of the next interaction
    Real(dp) :: sig_sum  !running sum of mass attenuation coefficients
    Real(dp) :: alpha  !ratio of photon energy to rest mass energy of an electron
    Real(dp) :: mu0 !cosine of the polar scattering angle
    Real(dp) :: w0  ![rad] azimuthal scattering angle
    Real(dp) :: b_hat(1:3)  ![1 cm] basis vector for correlating scattered direction to incident direction
    Real(dp) :: c_hat(1:3)  ![1 cm] basis vector for correlating scattered direction to incident direction
    Real(dp) :: nu  !sine of the polar scattering angle
    Real(dp) :: fp
    Integer :: i,j,s

    !Initialize random number generator
    Call RNG%seed(seed=default_seed)
    !Read in MAC values
    Call Read_MAC('resources/mac/07-N.txt',n_sig,hv_sigs,sigs)
    !Set up detector energy grid
    Do i = 0,n_E
        E_list(i) = 10._dp**(-4._dp + 0.1_dp*Real(i,dp))
    End Do
    !Initialize tally arrays
    n_d = Size(Rd)
    p_list_size = 2*n_d
    Allocate(p_list(1:2,1:p_list_size,1:n_d))
    p_list = 0._dp
    Allocate(p_listed(1:n_d))
    p_listed = 0
    Allocate(f(1:2,1:n_d))
    f = 0._dp
    Allocate(f_list(1:2,1:n_E,1:n_d))
    f_list = 0._dp
    !Run n_p histories
    Do p = 1,n_p
        If (Mod(p,1000) .EQ. 0) Write(*,'(A,F6.2,A)',ADVANCE='NO') creturn,100._dp*Real(p,dp)/Real(n_p,dp),'%'
        ! start a photon at the source
        r = (/ 0._dp , 0._dp , 0._dp /)
        hv = hv0_pdf(RNG)
        omega_hat = Isotropic_Omega_Hat(RNG)
        t = 0._dp
        w = 1._dp
        r0 = 0._dp
        zeta0 = 1._dp
        Do
            ! move the photon to the location of the next interaction
            tau = -Log(1._dp - RNG%r())
            sig = interp_sig(hv,n_sig,hv_sigs,sigs)
            sigT = Sum(sig)
            d = tau / (rho0 * sigT)
            r = r + d * omega_hat
            t = t + d / speed_c
            !check for and tally crossings of the detector surfaces
            r1 = Vector_Length(r)
            Do s = 1,n_d
                If (r0 .GT. Rd(s)) Then !start outside...
                    Rca = r0 * Sqrt(1._dp - zeta0**2)
                    If (r1 .LT. Rd(s)) Then !finish inside, one crossing
                        zeta = Sqrt(1._dp - (Rca/Rd(s))**2)
                        If (p_listed(s)+1 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                        p_listed(s) = p_listed(s) + 1
                        p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                    Else If (zeta0 .LT. 0._dp) Then !finish outside, but initially traveling inward...
                        If (Rca .LT. Rd(s)) Then !radius of close approach is inside, potential for two crossings exists...
                            dca = Sqrt(r0**2 - Rca**2)
                            If (dca .LT. d) Then !distance traveled exceeds that required to reach Rca, two crossings
                                zeta = Sqrt(1._dp - (Rca/Rd(s))**2)
                                If (p_listed(s)+2 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                                p_listed(s) = p_listed(s) + 1
                                p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                                dz = Rd(s) - Rca
                                zeta = Sqrt(2._dp*Rca*dz + dz**2 ) / Rd(s)
                                p_list(2,p_listed(s),s) = p_list(2,p_listed(s),s) + w*Min(1./zeta,10._dp)
                            End If
                        End If
                    End If
                Else !start inside...
                    If (r1 .GT. Rd(s)) Then !finish outside, one crossing
                        If (zeta0 .GT. 0._dp) Then
                            dz = Rd(s) - r0
                            zeta = Sqrt((r0*zeta0)**2 + 2._dp*r0*dz + dz**2) / Rd(s)
                        Else
                            Rca = r0 * Sqrt(1._dp - zeta0**2)
                            dz = Rd(s) - Rca
                            zeta = Sqrt(2._dp*Rca*dz + dz**2) / Rd(s)
                        End If
                        If (p_listed(s)+1 .GT. p_list_size) Call Resize_p_list(n_d,p_list_size,p_list)
                        p_listed(s) = p_listed(s) + 1
                        p_list(:,p_listed(s),s) = (/ hv , w*Min(1./zeta,10._dp) /)
                    End If
                End If
            End Do
            r0 = r1
            ! interact the photon in the medium
            xi = RNG%r() * sigT
            sig_sum = 0._dp
            Do interaction = 1,4
                sig_sum = sig_sum + sig(interaction)
                If (sig_sum .GT. xi) Exit
            End Do
            If (interaction .EQ. photoelectric_absorption) Exit
            Select Case (interaction)
                Case (incoherent_scatter,coherent_scatter)
                    If (interaction .EQ. incoherent_scatter) Then
                        alpha = hv / m0c2
                        mu0 = mu_KN(RNG,alpha)
                        hv = hv / (1._dp + alpha*(1._dp - mu0))
                    Else !coherent_scatter
                        mu0 = mu_Thomp(RNG)
                    End If
                    w0 = Isotropic_Azimuth(RNG)
                    If ( Abs(Dot_Product(omega_hat,Z_hat)) .LT. 0.5_dp ) Then
                        b_hat = Unit_Vector(Cross_Product(omega_hat,Z_hat))
                    Else
                        b_hat = Unit_Vector(Cross_Product(omega_hat,Y_hat))
                    End If
                    c_hat = Cross_Product(omega_hat,b_hat)
                    nu = Sqrt(1._dp - mu0**2)
                    omega_hat = mu0*omega_hat + nu*Cos(w0)*b_hat + nu*Sin(w0)*c_hat
                Case (pair_production)
                    omega_hat = Isotropic_Omega_hat(RNG)
                    hv = m0c2
                    w = w * 2._dp
            End Select
            zeta0 = Dot_Product(Unit_Vector(r),omega_hat)
        End Do
        !collect and record tallies from the history that just finished
        Do s = 1,n_d
            If (p_listed(s) .EQ. 0) Cycle
            i = p_listed(s)
            j = n_E
            fp = 0._dp
            Do
                j = Bisect_Search(p_list(1,i,s),E_list(0:j),n_E)
                Do
                    If (i .GT. 1) Then
                        If (p_list(1,i-1,s) .GT. E_list(j-1)) Then
                            p_list(1,i-1,s) = p_list(1,i-1,s) + p_list(1,i,s)
                            i = i - 1
                        Else
                            Exit
                        End If
                    Else
                        Exit
                    End If
                End Do
                fp = fp + p_list(2,i,s)
                f_list(:,j,s) = f_list(:,j,s) + (/ p_list(2,i,s) , p_list(2,i,s)**2 /)
                If (i .EQ. 1) Exit
                i = i - 1
            End Do
            f(:,s) = f(:,s) + (/ fp , fp**2 /)
        End Do
        p_listed = 0
        p_list = 0._dp
    End Do
    Write(*,*)
    Write(*,*)
    Do i = 1,n_d
        Write(*,'(2ES23.15)') f(1,i) / (n_p*4._dp*Pi*Rd(i)**2) , coeff_of_var(n_p,f(1,i),f(2,i))
    End Do
End Subroutine Transport_Photons

Subroutine Read_MAC(file_name,n,ln_hv,ln_sig)
    Use Kinds, Only: dp
    Use FileIO_Utilities, Only: Output_Message
    Implicit None
    Character(*), Intent(In) :: file_name
    Integer, Intent(Out) :: n
    Real(dp), Allocatable, Intent(Out) :: ln_hv(:)
    Real(dp), Allocatable, Intent(Out) :: ln_sig(:,:)
    Integer :: i
    Character(5) :: ctrash
    Real(dp) :: trash
    Integer :: unit, stat

    !open the file
    Open(NEWUNIT = unit , FILE = file_name , STATUS = 'OLD' , ACTION = 'READ' , IOSTAT = stat)
    If (stat .NE. 0) Call Output_Message('MACs: Read_MAC:  Error opening file '//file_name//', stat=',stat)
    !count the number of entries
    n = 0
    Do
        Read(unit,'(A5,2ES13.5E2)' , IOSTAT = stat) ctrash,trash,trash
        If (stat .LT. 0) Exit
        n = n + 1
    End Do
    Rewind(unit)
    !Allocate arrays
    Allocate(ln_hv(1:n))
    ln_hv = -1._dp
    Allocate(ln_sig(1:4,1:n))
    ln_sig = -1._dp
    !Read values
    Do i = 1,n
        Read(unit,'(7ES10.3E2)') ln_hv(i) ,    & !energy [MeV]
                               & ln_sig(3,i) , & !MAC for coherent scatter [g/cm^2]
                               & ln_sig(2,i) , & !MAC for incoherent scatter [g/cm^2]
                               & ln_sig(1,i) , & !MAC for photoelectric absorption [g/cm^2]
                               & ln_sig(4,i) , & !MAC for pair production (nuclear) [g/cm^2]
                               & trash           !MAC for pair production (atomic) [g/cm^2]
        ln_sig(4,i) = ln_sig(4,i) + trash  !merge pair production (nuclear/atomic) MACs
    End Do
    !close the file
    Close(unit)
    !convert to keV
    ln_hv = ln_hv * 1000._dp
    !convert to log scale
    ln_hv = Log(ln_hv)
    Where (ln_sig .NE. 0._dp)  !there are likely zeros in the pair production column
        ln_sig = Log(ln_sig)
    Elsewhere
        ln_sig = -33._dp
    End Where
End Subroutine

Function hv0_pdf_Line_N(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 7
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_N
Function hv0_pdf_Line_Fe(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 26
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_Fe
Function hv0_pdf_Line_Pb(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 82
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp

    hv = hv_mode
End Function hv0_pdf_Line_Pb
Function hv0_pdf_Line_5(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG

    hv = 5.E3_dp
End Function hv0_pdf_Line_5

Function Planck(hv,kT) Result(p)
    !Returns the intensity at E [keV] of a normalized Plankian blackbody with temperature kT [keV]
    Use Kinds, Only: dp
    Use Global, Only: Pi
    Implicit None
    Real(dp) :: p
    Real(dp), Intent(In) :: hv
    Real(dp), Intent(In) :: kT
    Real(dp):: b,b_sq
    Real(dp), Parameter :: c1 =  1._dp / 47900160._dp
    Real(dp), Parameter :: c2 = -1._dp / 1209600._dp
    Real(dp), Parameter :: c3 =  1._dp / 30240._dp
    Real(dp), Parameter :: c4 = -1._dp / 720._dp
    Real(dp), Parameter :: c5 =  1._dp / 12._dp
    Real(dp), Parameter :: c6 = -1._dp / 2._dp
    
    b = hv / kT
    If (b .LE. 0._dp) Then
        p = 0._dp
        Return
    Else If (b .LT. 0.1_dp) Then
        !for small values of z use a series expansion approximation
        !No doubt, a more compact or recursive form exists, but the following correctly accumulates ~20 digits of precision
        b_sq = b**2
        p = b_sq * c1 + c2
        p = p * b_sq + c3
        p = p * b_sq + c4
        p = p * b_sq + c5
        p = p * b + c6
        p = p * b + 1._dp
        p = p * b_sq * 15._dp / (kT*(Pi**4))
    Else
        !Bridgman 5-7
        p = 15._dp * hv**3 / ((Exp(b) - 1._dp) * (kT*pi)**4)
    End If
End Function Planck

Function hv0_pdf_Planck_N(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 7
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_N
Function hv0_pdf_Planck_Fe(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 26
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_Fe
Function hv0_pdf_Planck_Pb(RNG) Result(hv)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Use Global, Only: Pi
    Real(dp) :: hv
    Type(ACORN_Type), Intent(InOut) :: RNG
    Integer, Parameter :: Z = 82
    Real(dp), Parameter :: Z2_to_kT = 13.6e-3_dp * 2._dp / 3._dp
    Real(dp), Parameter :: kT = Z2_to_kT * Real(Z**2,dp)
    Real(dp), Parameter :: hv_mode = kT * 2.8214393721220789_dp
    Real(dp), Parameter :: hv_mode_P = 15._dp * hv_mode**3 / ((Exp(hv_mode/kT) - 1._dp) * (kT*pi)**4)
    Real(dp), Parameter :: hv_max = 10._dp * kT

    Do
        hv = RNG%r() * hv_max
        If (Planck(hv,kT) .LT. RNG%r() * hv_mode_P) Exit
    End Do
End Function hv0_pdf_Planck_Pb

Function Isotropic_Azimuth(RNG) Result(w)
    !Returns w (or omega), an azimuthal scattering angle uniformly distributed [0,2pi)
    Use Kinds, Only: dp
    Use Global, Only: TwoPi
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: w
    Type(ACORN_Type), Intent(InOut) :: RNG

    w = TwoPi * RNG%r()
End Function Isotropic_Azimuth

Function Isotropic_mu(RNG) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta),
    ! uniformly distributed [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG

    mu = 2._dp * RNG%r() - 1._dp
End Function Isotropic_mu

Function mu_w_2_OmegaHat(xi,w) Result(Omegahat)
    !Converts angles (polar cosine and azimuthal rotation) to a unit vector
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: OmegaHat(1:3)
    Real(dp), Intent(In) :: xi !cosine of the polar scattering angle
    Real(dp), Intent(In) :: w  !azimuthal scattering angle
    Real(dp) :: nu,mu,eta

    nu = Sqrt(1._dp - xi**2)
    mu = nu * Cos(w)
    eta = nu * Sin(w)
    OmegaHat = (/ mu, eta, xi /)
End Function mu_w_2_OmegaHat

Function Isotropic_Omega_hat(RNG) Result(OmegaHat)
    !Reutrns an isotropic random direction
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp):: OmegaHat(1:3)
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp):: xi  ! z direction cosine = cos(theta), theta is angle from zHat
    Real(dp):: w   ! angle from xHat toward yHat around zHat

    xi = Isotropic_mu(RNG)
    w = Isotropic_Azimuth(RNG)
    OmegaHat = mu_w_2_OmegaHat(xi,w)
End Function Isotropic_Omega_hat

Function Bisect_Search(x, xx, n) Result (j)
    ! Finds index j, such that x is in the interval ( xx(j-1),xx(j) ]
    ! The list xx must be sorted in acsending order
    ! j=1 or j=n+1 indicates that x is not in the range [ xx(1),xx(n) ]
    Use Kinds, Only: dp
    Implicit None
    Integer :: j
    Real(dp), Intent(In) :: x
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: xx(1:n)
    Integer :: jl,jh,jm

    !First, check for endpoints
    If (x .EQ. xx(1)) Then
        j = 2
        Return
    Else If (x .EQ. xx(n)) Then
        j = n
        Return
    End If
    !Otherwise, search the list
    jl = 0
    jh = n + 1
    Do
        If (jh - jl .LE. 1) Exit
        jm = (jh + jl) / 2  !integer arithmetic
        If (x .LT. xx(jm)) Then
            jh = jm
        Else
            jl = jm
        End If
    End Do
    j = jh
End Function Bisect_Search

Elemental Function Lin_Interp(x,x1,x2,y1,y2) Result(y)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: y
    Real(dp), Intent(In) :: x
    Real(dp), Intent(In) :: x1,x2
    Real(dp), Intent(In) :: y1,y2

    y = y1 + (y2-y1) * (x-x1) / (x2-x1)
End Function Lin_Interp

Function interp_sig(hv,n_sig,hv_sig,sig)
    Use Kinds, Only: dp
    Use Global, Only: m0c2
    Implicit None
    Real(dp) :: interp_sig(1:4)
    Real(dp), Intent(In) :: hv
    Integer, Intent(In) :: n_sig
    Real(dp), Intent(In) :: hv_sig(1:n_sig)
    Real(dp), Intent(In) :: sig(1:4,1:n_sig)
    Real(dp), Parameter :: pp_thresh = 2._dp * m0c2
    Real(dp) :: ln_hv
    Integer :: j

    ln_hv = Log(hv)
    j = Bisect_Search(ln_hv,hv_sig,n_sig)
    If (j .EQ. 0) j = 1
    interp_sig = Exp( Lin_Interp(ln_hv,hv_sig(j-1),hv_sig(j),sig(:,j-1),sig(:,j)) )
    If (hv .LT. pp_thresh) interp_sig(4) = 0._dp
End Function interp_sig

Function mu_Thomp(RNG) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta),
    ! distributed on the Thompson distribution on [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp) :: xi
    Real(dp) :: a
    Real(dp), Parameter :: one_third = 1._dp / 3._dp
    Real(dp), Parameter :: two_thirds = 2._dp / 3._dp

    xi = RNG%r()
    a = 2._dp - 4._dp * xi + Sqrt(5._dp + 16._dp*xi*(xi-1._dp))
    mu = (1._dp - a**two_thirds) / (a**one_third)
End Function mu_Thomp

Function mu_KN(RNG,alpha) Result(mu)
    !Returns mu, cosine of the polar scattering angle (theta), 
    ! distributed on the Klein-Nishina Distribution on [1,-1) corresponding to angles [0,pi)
    Use Kinds, Only: dp
    Use PRNGs, Only: ACORN_Type
    Implicit None
    Real(dp) :: mu
    Type(ACORN_Type), Intent(InOut) :: RNG
    Real(dp), Intent(In) :: alpha ![] energy of the incident photon, in rest mass energy of an electron
    Real(dp) :: kn  !Klein Nishina distribution evaluated at a proposed mu
    Real(dp) :: mu2p1 !(mu**2 + 1._dp) intermediate value for computing kn
    Real(dp) :: amamup1 !( 1._dp + a * (1._dp - mu) ) intermediate value for computing kn

    Do  !sampling by geometric rejection
        !sample mu from Thompson scattering distribution
        mu = mu_Thomp(RNG)
        !compute the value of the KN distribution at this mu
        !This is the KN distributon scaled on 0 < f(mu) < TWO  (using this scaled value eliminates a multiplication here and a
        !                                                       multiplication in the following check for acceptance)
        mu2p1 = 1._dp + mu**2
        amamup1 = 1._dp + alpha * (1._dp - mu)
        kn = mu2p1 * ( 1._dp + (alpha * (1._dp - mu))**2 / (mu2p1 * amamup1) ) / amamup1**2
        !If the probability of KN at this mu exceeds the probability of the sampled dist, then accept the value and exit
        If ( kn .GT. mu2p1*RNG%r() ) Exit
        !Otherwise, cycle and repeat the sampling process
    End Do
End Function mu_KN

Pure Function Vector_Length(v) Result(d)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: d
    Real(dp), Intent(In) :: v(1:3)

    d = Sqrt(Sum(v*v))
End Function Vector_Length

Pure Function Unit_Vector(v) Result(hat)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: hat(1:3)
    Real(dp), Intent(In) :: v(1:3)

    hat = v / Vector_Length(v)
End Function Unit_Vector

Pure Function Cross_Product(v,w) Result(u)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: u(1:3)
    Real(dp), Intent(In) :: v(1:3),w(1:3)

    u = (/ v(2)*w(3) - v(3)*w(2) , &
         & v(3)*w(1) - v(1)*w(3) , &
         & v(1)*w(2) - v(2)*w(1) /)
End Function Cross_Product

Subroutine Resize_p_list(nd,n,p_list)
    Use Kinds, Only: dp
    Implicit None
    Integer, Intent(In) :: nd
    Integer, Intent(InOut) :: n
    Real(dp), Allocatable, Intent(InOut) :: p_list(:,:,:)
    Real(dp), Allocatable :: p_swap(:,:,:)

    Allocate(p_swap(1:2,1:n,1:nd))
    p_swap = p_list
    Deallocate(p_list)
    Allocate(p_list(1:2,1:2*n,1:nd))
    p_list = 0._dp
    p_list(1:2,1:n,1:nd) = p_swap
    Deallocate(p_swap)
    n = 2 * n
End Subroutine Resize_p_list

Function std_err( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: std_err
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    If (n .EQ. 1) Then
        std_err = mu
    Else
        std_err = Sqrt( (mu_sq - mu**2 / Real(n,dp)) / Real(n-1,dp) )
    End If
End Function std_err

Function std_err_otm( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: std_err_otm
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    std_err_otm = std_err(n,mu,mu_sq) / Sqrt(Real(n,dp))
End Function std_err_otm

Function coeff_of_var( n , mu , mu_sq )
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: coeff_of_var
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: mu
    Real(dp), Intent(In) :: mu_sq

    If (n .EQ. 1) Then
        coeff_of_var = 1._dp
    Else
        coeff_of_var = std_err_otm(n,mu,mu_sq) / (mu/Real(n,dp))
    End If
End Function coeff_of_var

End Program