calc_mixed_end_contr_approx Subroutine

    subroutine calc_mixed_end_contr_approx(t, csf_i, excitInfo, integral)
        ! for the approx. mixed end contribution i "just" need to
        ! calculate the correct matrix element influences
        integer(n_int), intent(in) :: t(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        HElement_t(dp), intent(out) :: integral
        character(*), parameter :: this_routine = "calc_mixed_end_contr_approx"

        integer :: st, se, en, elecInd, holeInd, step, sw, i
        real(dp) :: top_cont, mat_ele, stay_mat, end_mat
        logical :: above_flag

        ! do as much stuff as possible beforehand
        st = excitInfo%fullStart
        se = excitInfo%secondStart
        en = excitInfo%fullEnd
        if (excitInfo%typ == excit_type%fullstop_L_to_R) then
            elecInd = st
            holeInd = se
        else if (excitInfo%typ == excit_type%fullstop_R_to_L) then
            elecInd = se
            holeInd = st
        else
            call stop_all(this_routine, "should not be here!")
        end if

        integral = h_cast(0.0_dp)

        step = csf_i%stepvector(en)

        ! i am sure the last switch happens at the full-stop!
        sw = en

        if (en < nSpatOrbs) then
            select case (step)
            case (1)
                if (isOne(t, en)) then
                    top_cont = -Root2 * sqrt((csf_i%B_real(en) + 2.0_dp) / &
                                             csf_i%B_real(en))

                else
                    top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))

                end if
            case (2)
                if (isOne(t, en)) then
                    top_cont = -Root2 / sqrt(csf_i%B_real(en) * (csf_i%B_real(en) + 2.0_dp))

                else
                    top_cont = Root2 * sqrt(csf_i%B_real(en) / &
                                            (csf_i%B_real(en) + 2.0_dp))
                end if

            case default
                call stop_all(this_routine, "wrong stepvalues!")

            end select

            if (.not. near_zero(top_cont)) then

                above_flag = .false.
                mat_ele = 1.0_dp

                do i = en + 1, nSpatOrbs
                    if (csf_i%Occ_int(i) /= 1) cycle

                    ! then check if thats the last step
                    if (csf_i%stepvector(i) == 2 .and. csf_i%B_int(i) == 0) then
                        above_flag = .true.
                    end if

                    ! in the other routine i check if the orbital pgen
                    ! is 0 for the above orbitals.. do I need to do that
                    ! also here?? or is this implicit if the matrix
                    ! element will be 0??

                    step = csf_i%stepvector(i)

                    call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, csf_i%B_real(i), &
                                                1.0_dp, x1_element=stay_mat)

                    call getMixedFullStop(step, step, 0, csf_i%B_real(i), &
                                          x1_element=end_mat)

                    ! this check should never be true, but just to be sure
                    if (near_zero(stay_mat)) above_flag = .true.

                    if (.not. near_zero(end_mat)) then
                        integral = integral + end_mat * mat_ele * &
                                   (get_umat_el(i, holeInd, elecInd, i) + &
                                    get_umat_el(holeInd, i, i, elecInd)) / 2.0_dp
                    end if

                    if (above_flag) exit

                    ! otherwise update your running matrix element vars
                    mat_ele = mat_ele * stay_mat

                end do

                integral = integral * top_cont
            end if
        end if

    end subroutine calc_mixed_end_contr_approx