get_ddpt_anisotropy Subroutine

subroutine get_ddpt_anisotropy(plasma, v1, v3, kappa)
    !+ Gets d(d,p)T anisotropy defecit/enhancement factor for a beam interacting with a target plasma
    type(LocalProfiles), intent(in)         :: plasma
        !+ Plasma Paramters
    real(Float64), dimension(3), intent(in) :: v1
        !+ Beam velocity [cm/s]
    real(Float64), dimension(3), intent(in) :: v3
        !+ Charged Fusion Product velocity [cm/s]
    real(Float64), intent(out)              :: kappa
        !+ Anisotropy factor
    !+ Reference: Eq. (1) and (3) of NIM A236 (1985) 380

    real(Float64), dimension(3,12) :: abc
    real(Float64), dimension(13)   :: bhcor !!! 13?
    real(Float64), dimension(12)   :: e, a, b, c
    real(Float64), dimension(3)    :: vcm, v3cm, vrel
    type(InterpolCoeffs1D) :: c1D

    real(Float64) :: ai, bi, ci, b1, b2, b11, b12, b21, b22, cos_phi, sin_phi
    real(Float64) :: eb, e1com, vnet_square, cos_theta, k, KE, Q, mp, k0, JMeV
    integer :: ei, i, err_status

    !! Calculate effective beam energy
    vrel = v1-plasma%vrot ![cm/s]
    vnet_square=dot_product(vrel, vrel) ![(cm/s)**2]
    eb = v2_to_E_per_amu*fbm%A*vnet_square ![keV]

    !!Calculate anisotropy enhancement/deficit factor
    JMeV = 1.60218d-13 ! Conversion factor from MeV to Joules
    mp = H1_amu*mass_u  ![kg]
    Q = 4.04*JMeV ![J]

    vcm = 0.5*(v1+plasma%vrot) ![cm/s]
    KE = 0.5*mp*vnet_square*1.d-4  ! [J] C-O-M kinetic energy
    k0 = norm2(vcm) * sqrt(2*mp/(3*(Q+KE)))*100.d0 ![(cm/s)**2]
    if ((norm2(vcm)*norm2(v3)).gt.0.d0) then
        cos_phi = dot_product(vcm, v3) / (norm2(vcm)*norm2(v3))
        sin_phi = sin(acos(cos_phi))

        if (abs(k0*sin_phi).le.1) then
            cos_theta = cos_phi*sqrt(1-(k0*sin_phi)**2) - k0*sin_phi**2
        else
            cos_theta = 0.d0
        endif

    else
        cos_theta = 0.d0
    endif

    !Brown-Jarmie coefficients in Table I with prepended isotropic low-energy extrapolated point
    e = [0.0,19.944,29.935,39.927,49.922,59.917,69.914,79.912,89.911,99.909,109.909,116.909]
    a = [0.0,0.0208,0.0886,0.1882,0.3215,0.4636,0.6055,0.7528,0.8976,1.041,1.166,1.243]
    b = [0.0,0.0023,0.0184,0.0659,0.076,0.168,0.251,0.250,0.378,0.367,0.521,0.55]
    c = [0.0,0.0,0.0,0.0,0.048,0.021,0.039,0.162,0.142,0.276,0.203,0.34]
    e(1) = 10.0
    a(1) = 0.00903775/(4*pi)
    abc(1,:) = a ; abc(2,:) = b ; abc(3,:) = c

    !Correction factor to make it consistent with Bosch & Hale
    bhcor=[1.0,1.02614,.98497,1.006,1.015,1.012,1.0115,1.021,1.012,1.012,1.016,1.018,1.012]
    do i=1,12
        abc(:,i) = abc(:,i)*bhcor(i)
    enddo

    e1com=0.5d0*eb
    call interpol_coeff(e, e1com, c1D, err_status)

    ei = c1D%i
    b1 = c1D%b1
    b2 = c1D%b2

    ai = b1*abc(1,ei) + b2*abc(1,ei+1)
    bi = b1*abc(2,ei) + b2*abc(2,ei+1)
    ci = b1*abc(3,ei) + b2*abc(3,ei+1)

    kappa = (ai + bi*cos_theta**2 + ci*cos_theta**4) / (ai+bi/3.d0+ci/5.d0)

end subroutine get_ddpt_anisotropy