elbow_coeff Function

    function elbow_coeff(r_elbow, ang_elbow, lambda, dh)

        !! Function calculates elbow bends coefficients for pressure drop
        !! calculations.
        !!
        !! author: G. Graham, CCFE
        !!
        !! References:
        !!
        !!      [Ide1969]   Idel'Cik, I. E. (1969), Memento des pertes de charge,
        !!                  Collection de la Direction des Etudes et Recherches d'Electricité de France.
        !!
        !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        use constants, only: pi

        implicit none

        !! Function return parameter !!!!!

        !! Elbow coefficient for pressure drop calculation
        real(dp) :: elbow_coeff

        !! Arguments !!!!!!!!!!!!!!!!!!!!!

        !! Pipe elbow radius (m)
        real(dp), intent(in) :: r_elbow

        !! Pipe elbow angle (degrees)
        real(dp), intent(in) :: ang_elbow

        !! Darcy Friction Factor
        real(dp), intent(in) :: lambda

        !! Hydraulic Diameter (m)
        real(dp), intent(in) :: dh

        !! Local variables !!!!!!!!!!!!!!!

        !! A and B from singularity coefficient equation
        real(dp)  :: a, b

        !! Singularity and friction coefficients
        real(dp)  :: xifn, xift, ximn, ximt

        !! Ratio of elbow radius to Darcy friction factor
        real(dp)  :: r_ratio

        !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        !! Pressure loss for bends calculated using equivilent lengths.
        !! coeff = lambda * L/D --> get equiv. of L

        if (ang_elbow==90) then
            a = 1.0D0
        else if (ang_elbow<70) then
            a = 0.9D0 * sin(ang_elbow * pi/180.0D0)
        else if (ang_elbow>100) then
            a = 0.7D0 + (0.35D0 * sin((ang_elbow/90.0D0) * (pi/180.0D0)))
        else
            write(*,*) 'No formula for 70 <= elbow angle(deg) <= 100, only 90 deg option available in this range.'
            stop 1
        endif

        r_ratio = r_elbow/dh
        if (r_ratio > 1) b = 0.21D0/r_ratio**0.5
        if (r_ratio < 1) b = 0.21D0/r_ratio**2.5
        if (r_ratio == 1) b = 0.21D0

        !! Singularity
        ximt = a * b

        !! Friction
        xift = 0.0175D0 * lambda * (r_elbow/dh) * ang_elbow

        !! Elbow Coefficient
        elbow_coeff = ximt + xift

    end function elbow_coeff