!****************************************************************************************
! Module 'damageConfigurations' contains two subroutines representing the two major steps
! in the damage generation algorithm: 1) randomly distribute elementary damages in a
! bi-stranded segment and 2) group elementary damages into lesions.
    MODULE damageConfigurations

    IMPLICIT NONE

    PRIVATE ::
    PUBLIC :: fstStrand, sndStrand, nd, segment, distributeEDs, detectLesions

    INTEGER, PARAMETER ::    fstStrand = 1,        &    ! 1st DNA strand
                            sndStrand = 2,        &    ! 2nd DNA strand

                            nd = 0,                &    ! not damaged
                            sb = 1,                &    ! strand break
                            bd = 2,                &    ! base damage

                            ! parameters for low-LET radiation (electrons)
                            SigSb = 900,        &
                            SigBd = 2700,        &
                            Nseg = 70000,        &
                            Nmin = 8

    INTEGER, DIMENSION(Nseg,2) :: segment


    CONTAINS
!****************************************************************************************
! Distributes elementary damages uniformly throughout two columns of an array of length Nseg
        SUBROUTINE distributeEDs
        INTEGER :: i, edsPos, strand
        REAL :: rn
        ! distribute strand breaks
        DO i = 1, SigSb
            DO
                ! select an integer number in the range [1, Nseg]
                CALL RANDOM_NUMBER( rn )
                edsPos = 1 + INT( Nseg * rn )
                ! select strand
                CALL RANDOM_NUMBER( rn )
                IF ( rn < 0.5 ) THEN
                    strand = fstStrand
                ELSE
                    strand = sndStrand
                END IF
                ! if not damaged, exit and record 'sb' at the current location
                IF ( segment(edsPos,strand) == nd ) EXIT
            END DO
            segment(edsPos,strand) = sb
        END DO
        ! distribute base damages
        DO i = 1, SigBd
            DO
                ! select an integer number in the range [1, Nseg]
                CALL RANDOM_NUMBER( rn )
                edsPos = 1 + INT( Nseg * rn )
                ! select strand
                CALL RANDOM_NUMBER( rn )
                IF ( rn < 0.5 ) THEN
                    strand = fstStrand
                ELSE
                    strand = sndStrand
                END IF
                ! if not damaged, exit and record 'bd' at the current location
                IF ( segment(edsPos,strand) == nd ) EXIT
            END DO
            segment(edsPos,strand) = bd
        END DO
        END SUBROUTINE distributeEDs


!****************************************************************************************
! Detects leions. Two elementary damages are recognized as separate lesions if there are
! at least Nmin undamaged base pairs between them. Outputs the initial (firstLine) and final
! (lastLine) positions of the next detected lesion. When the last lesion in the segment is
! detected, the logical variable 'more' is set to FALSE
        SUBROUTINE detectLesions( firstLine, lastLine, more )
        INTEGER, INTENT(OUT) :: firstLine, lastLine
        LOGICAL, INTENT(OUT) :: more
        INTEGER :: i, counter
        INTEGER, SAVE :: startingLine = 1

        more = .TRUE.
        i = startingLine
        ! scan the DNA segment from one end and locate the first ED on either (or both) strand
        DO WHILE ( (segment(i,fstStrand) == nd).AND.(segment(i,sndStrand) == nd) )
            i = i + 1
            IF ( i > Nseg ) THEN    ! so that array boundary is not exceeded
                more = .FALSE.
                startingLine = 1
                RETURN
            END IF
        END DO
        ! set the start of the lesion to the location of the first ED(s)
        firstLine = i
        DO
            lastLine = i    ! set the end of the lesion to the location of the upstream ED
            i = i + 1        ! start counting undamaged base pairs from the next position
            IF ( i > Nseg ) THEN    ! so that array boundary is not exceeded
                more = .FALSE.
                startingLine = 1
                RETURN
            END IF
            ! count undamaged base pairs
            counter = 0
            DO WHILE ( (segment(i,fstStrand) == nd).AND.(segment(i,sndStrand) == nd) )
                i = i + 1
                IF ( i > Nseg ) THEN    ! so that array boundary is not exceeded
                    more = .FALSE.
                    startingLine = 1
                    RETURN
                END IF
                counter = counter + 1
            END DO
            IF ( counter >= Nmin ) EXIT    ! the current lesion is over; return to the main program
        END DO
        startingLine = i

        END SUBROUTINE detectLesions


    END MODULE damageConfigurations