calc_mixed_start_contr_approx Subroutine

    subroutine calc_mixed_start_contr_approx(t, csf_i, excitInfo, integral)
        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_start_contr_approx"

        integer :: st, se, en, elecInd, holeInd, sw, step, i
        real(dp) :: bot_cont, mat_ele, stay_mat, start_mat
        logical :: below_flag

        st = excitInfo%fullStart
        se = excitInfo%firstEnd
        en = excitInfo%fullEnd
        ! depending on the type of excitaiton, calculation of orbital pgens
        ! change
        if (excitInfo%typ == excit_type%fullStart_L_to_R) then
            elecInd = en
            holeInd = se
        else if (excitInfo%typ == excit_type%fullstart_R_to_L) then
            elecInd = se
            holeInd = en
        else
            call stop_all(this_routine, "should not be here!")
        end if

        ! I am sure first switch is at full-start
        sw = st

        ! what can i precalculate beforehand?
        step = csf_i%stepvector(st)

        integral = h_cast(0.0_dp)

        if (step == 1) then

            ASSERT(isTwo(t, st))

            bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) - 1.0_dp) * &
                                       (csf_i%B_real(st) + 1.0_dp)))

        else

            ASSERT(isOne(t, st))

            bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) + 1.0_dp) * &
                                       (csf_i%B_real(st) + 3.0_dp)))

        end if

        if (.not. near_zero(bot_cont)) then

            mat_ele = 1.0_dp
            below_flag = .false.

            do i = st - 1, 1, -1
                if (csf_i%Occ_int(i) /= 1) cycle

                ! then check if thats the last stepvalue to consider
                if (csf_i%stepvector(i) == 1 .and. csf_i%B_int(i) == 1) then
                    below_flag = .true.
                end if

                ! then deal with the matrix element and branching probabilities
                step = csf_i%stepvector(i)

                ! get both start and staying matrix elements -> and update
                ! matrix element contributions on the fly to avoid second loop!
                call getDoubleMatrixElement(step, step, -1, gen_type%R, gen_type%L, csf_i%B_real(i), &
                                            1.0_dp, x1_element=start_mat)

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

                if (near_zero(stay_mat)) below_flag = .true.
                ! "normally" matrix element shouldnt be 0 anymore... still check

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

                ! also update matrix element on the fly
                mat_ele = stay_mat * mat_ele

            end do
        end if

        ! maybe I need to deal with the pgens here, if they are not
        ! correctly considered outside..

    end subroutine calc_mixed_start_contr_approx