c ----------------------------------------------------------------------
c This program finds the size of the restriction enzyme partial digests.
c
c   icut(i)... ith restriction sites (plus end points)
c   iseg   ... length of the segment
c   ntotal ... total number of bases
c   ncut   ... number of cuts or restriction sites
c
c Instructor: Nam Sun Wang
c ----------------------------------------------------------------------
      dimension icut(100)
      character resp*1

c program header -------------------------------------------------------
      print *, 'This program reads the restriction sites and         '
      print *, '  finds the size of the complete and partial digests.'
      print *, ' '

c Enter restriction sites ----------------------------------------------
      print *, 'Circular DNA (Y/n)? '
      read(*,500) resp
      print *, 'Number of bases in this DNA: '
      read(*,*) ntotal

      print *, 'Enter the restriction sites in ascending order.'
      print *, 'Enter 0 if there is no more.'

      do i=1, 100
        read(*,*) icut(i)
        if(icut(i) .le. 0) goto 11
      enddo

11    ncut = i-1

      print *, ' '
      print *, '--------------------------------'
      print *, 'Segment     Start   End  # bases'
      print *, '--------------------------------'

c Find partial digest --------------------------------------------------
      if (resp .eq. 'n' .or. resp .eq. 'N')then

c Linear DNA -----------------------------------------------------------
        do j=1, ncut+1
          do i=j, ncut+1
            ilo1  = i-j+1
            icutlo = icut(ilo1-1)
            if(ilo1 .eq. 1)      icutlo = 1
            icuthi = icut(i)-1
            if(i .eq. ncut+1) icuthi = ntotal
            iseg = icuthi - icutlo + 1
            print 650, ilo1, i, icutlo, icuthi, iseg
          enddo
        enddo

c Circular DNA ---------------------------------------------------------
      else
        do j=1, ncut
          do i=1, ncut
            ihi  = mod(i+j, ncut)
            ihi1 = mod(i+j-1, ncut)
            if (ihi  .eq. 0) ihi =ncut
            if (ihi1 .eq. 0) ihi1=ncut
            iseg = (i+j-1)/ncut*ntotal + icut(ihi) - icut(i)
            print 650, i, ihi1, icut(i), icut(ihi)-1, iseg
          enddo
        enddo
      endif

c Some formats ---------------------------------------------------------
500   format(a)
650   format(i5, '-', i2, 3x, 3i7)

      end