PROGRAM graph_sol
*     graphical solution for trajectory of logistic map
      IMPLICIT NONE
      REAL*8 r, x
      INTEGER IR_, iterate
      CALL GWopen(IR_, 0)
      CALL initial(x,r,iterate)
      CALL draw_function(r,iterate)
      CALL trajectory(x,r,iterate)          ! press any key to stop
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(x0,r,iterate)
      IMPLICIT NONE
      REAL*8 x0, r
      INTEGER iterate
      WRITE(*,'(A,$)') 'control parameter r = '
      READ(*,*) r
      WRITE(*,'(A,$)') 'initial value of x = '
      READ(*,*) x0
      WRITE(*,'(A,$)') 'iterate of f(x) = '
      READ(*,*) iterate
      END
C
      SUBROUTINE draw_function(r,iterate)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 delta, x, y, r, f
      REAL margin
      INTEGER IR_, i, nplot, iterate
      CALL GWsetpen(IR_, -1, -1, -1, -1)
      nplot = 200                           ! # of points at which function computed
      delta = 1.0D0/nplot
      margin = 0.1
      CALL GWindow(IR_, -margin, -margin, 1 + margin, 1 + margin)
      WRITE(Stmp1_,*) 'r = ', r
      CALL GWputtxt(IR_, 0.0, 1.0, Stmp1_)
      CALL GWline(IR_, 0.0, 0.0, 1.0, 1.0)  ! draw diagonal line y = x
      CALL GWmove2(IR_, 0.0, 1.0)           ! draw axes
      CALL GWline2(IR_, 0.0, 0.0)
      CALL GWline2(IR_, 1.0, 0.0)
      CALL GWsetpen(IR_, -1, -1, -1, -1)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      x = 0
      DO i = 0, nplot
         y = f(x,r,iterate)
         CALL GWline2(IR_, REAL(x), REAL(y))
         x = x + delta
      END DO
      END
C
      SUBROUTINE trajectory(x,r,iterate)
      IMPLICIT NONE
      REAL*8 x0, y0, x, y, r, f
      INTEGER IR_, k, iterate
      LOGICAL TBkeyinput
      y0 = 0
      x0 = x
      CALL GWsetpen(IR_, 16, -1, -1, -1)
      DO WHILE(.NOT. TBkeyinput())
         y = f(x,r,iterate)
         CALL GWline2(IR_, REAL(x0), REAL(y0))
         CALL GWline2(IR_, REAL(x0), REAL(y))
         CALL GWline2(IR_, REAL(y), REAL(y))
         CALL GWsetpen(IR_, -1, -1, -1, -1)
         x0 = y
         y0 = y
         x = y
      END DO
      CALL TBgetkey(k)
      END