      PROGRAM period
*     find fixed point of f(x) iterated p times
      IMPLICIT NONE
      REAL*8 epsilon, gleft, gright, r, x, xleft, xright, f
      INTEGER i, p
      LOGICAL Ltmp1_
      CALL initial(r,p,epsilon,xleft,xright,gleft,gright)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL bisection(r,p,xleft,xright,gleft,gright)
         Ltmp1_ = abs(xleft - xright) .GE. epsilon
      END DO
      x = 0.5D0*(xleft + xright)
      WRITE(*,*) 'explicit demonstration of period ', p, ' behavior'
      WRITE(*,*)
      WRITE(*,'(I12,F20.13)') 0, x          ! result
      DO i = 1, 2*p + 1
         x = f(x,r,1)
         WRITE(*,'(I12,F20.13)') i, x
      END DO
      END
C
      SUBROUTINE initial(r,p,epsilon,xleft,xright,gleft,gright)
      IMPLICIT NONE
      REAL*8 r, epsilon, xleft, xright, gleft, gright, f
      INTEGER p
      CHARACTER*80 done
      WRITE(*,'(A,$)') 'control parameter r = '
      READ(*,*) r
      WRITE(*,'(A,$)') 'period = '
      READ(*,*) p
      WRITE(*,'(A,$)') 'desired precision = '
      READ(*,*) epsilon
      done = 'no'
      DO WHILE(done .EQ. 'no')              ! do until zero between xleft and xright
         WRITE(*,'(A,$)') 'guess for xleft = '
         READ(*,*) xleft
         WRITE(*,'(A,$)') 'guess for xright = '
         READ(*,*) xright
         gleft = f(xleft,r,p) - xleft
         gright = f(xright,r,p) - xright
         IF(gleft*gright .LT. 0) THEN
            done = 'yes'
         ELSE
            WRITE(*,*) 'range does not necessarily enclose a root'
         END IF
      END DO
      END
C
      SUBROUTINE bisection(r,p,xleft,xright,gleft,gright)
*     midpoint between xleft and xright
      IMPLICIT NONE
      REAL*8 gmid, xmid, r, xleft, xright, gleft, gright, f
      INTEGER p
      xmid = 0.5D0*(xleft + xright)
      gmid = f(xmid,r,p) - xmid
      IF(gmid*gleft .GT. 0) THEN
         xleft = xmid                       ! change xleft
         gleft = gmid
      ELSE
         xright = xmid                      ! change xright
         gright = gmid
      END IF
      END
