calcDoubleR2L Subroutine

    subroutine calcDoubleR2L(ilut, csf_i, excitInfo, excitations, nExcits, &
                             posSwitches, negSwitches)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(inout) :: excitInfo
        integer(n_int), intent(out), allocatable :: excitations(:, :)
        integer, intent(out) :: nExcits
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        character(*), parameter :: this_routine = "calcDoubleR2L"

        integer :: iOrb, start1, start2, ende1, ende2
        type(WeightObj_t) :: weights
        real(dp) :: plusWeight, minusWeight, zeroWeight
        integer(n_int), allocatable :: tempExcits(:, :)
        !todo asserts

        ASSERT(.not. isThree(ilut, excitInfo%fullStart))
        ASSERT(.not. isZero(ilut, excitInfo%secondStart))
        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
        ASSERT(.not. isThree(ilut, excitInfo%fullEnd))

        start1 = excitInfo%fullStart
        start2 = excitInfo%secondStart
        ende1 = excitInfo%firstEnd
        ende2 = excitInfo%fullEnd

        ! todo: create correct weights:
        weights = init_fullDoubleWeight(csf_i, start2, ende1, ende2, negSwitches(start2), &
                                        negSwitches(ende1), posSwitches(start2), posSwitches(ende1), &
                                        csf_i%B_real(start2), csf_i%B_real(ende1))

        excitInfo%currentGen = excitInfo%firstGen

        ! then do single start:
        call createSingleStart(ilut, csf_i, excitInfo, posSwitches, negSwitches, &
                               weights, tempExcits, nExcits)

        ! and single update until semi start
        do iOrb = excitInfo%fullStart + 1, excitInfo%secondStart - 1
            call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, negSwitches, &
                              weights, tempExcits, nExcits)
        end do

        ! change weights... maybe need both single and double type weights
        ! maybe semistart is wrong here..
        weights = weights%ptr

        minusWeight = weights%proc%minus(negSwitches(start2), csf_i%B_real(start2), weights%dat)
        plusWeight = weights%proc%plus(posSwitches(start2), csf_i%B_real(start2), weights%dat)
        zeroWeight = weights%proc%zero(negSwitches(start2), posSwitches(start2), &
                                       csf_i%B_real(start2), weights%dat)

        ! then do lowering semi start
        call calcLoweringSemiStart(ilut, csf_i, excitInfo, &
                                   tempExcits, nExcits, plusWeight, minusWeight, zeroWeight)

        ! then do double excitation over double excitation region
        do iOrb = excitInfo%secondStart + 1, excitInfo%firstEnd - 1
            call doubleUpdate(ilut, csf_i, iOrb, excitInfo, weights, tempExcits, nExcits, &
                              negSwitches, posSwitches)
        end do

        ! update weights again:
        weights = weights%ptr
        minusWeight = weights%proc%minus(negSwitches(ende1), csf_i%B_real(ende1), weights%dat)
        plusWeight = weights%proc%plus(posSwitches(ende1), csf_i%B_real(ende1), weights%dat)

        ! then do lowering semi stop
        call calcRaisingSemiStop(ilut, csf_i, excitInfo, tempExcits, nExcits, plusWeight, &
                                 minusWeight)

        ! have to set used generators correctly
        excitInfo%currentGen = excitInfo%lastGen

        ! and then do final single region again
        do iOrb = excitInfo%firstEnd + 1, excitInfo%fullEnd - 1
            call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, negSwitches, &
                              weights, tempExcits, nExcits)
        end do

        ! and finally end step
        call singleEnd(ilut, csf_i, excitInfo, tempExcits, nExcits, excitations)

        ! that should be it...

    end subroutine calcDoubleR2L