PROGRAM Lyapunov

        PROGRAM Lyapunov
C
C
C Plots points in Lyapunov space
C
C J.J.Bunn 1992
C
        IMPLICIT INTEGER (A-Z)
C
        INCLUDE 'SYS$LIBRARY:DECW$XLIBDEF'
C
        RECORD /X$VISUAL/ VISUAL            ! visual type
        RECORD /X$SET_WIN_ATTRIBUTES/ XSWDA ! window attributes
        RECORD /X$GC_VALUES/ XGCVL          ! gc values
        RECORD /X$SIZE_HINTS/ XSZHN         ! hints
        RECORD /X$EVENT/ EVENT              ! input event
        RECORD /X$COLOR/ COLOUR             ! colours
C
        INTEGER GC
C
C       Initialize display id and screen id
C
        DPY = X$OPEN_DISPLAY()
        IF (DPY .EQ. 0) then
           status = lib$put_output('Error opening display')
           status = sys$exit(%val(1))
        endif
        SCREEN = X$DEFAULT_SCREEN_OF_DISPLAY(DPY)
C
C       Create the window
C
        DEPTH = X$DEFAULT_DEPTH_OF_SCREEN(SCREEN)
C
        CALL X$DEFAULT_VISUAL_OF_SCREEN(SCREEN,VISUAL)
C
        ATTR_MASK = X$M_CW_EVENT_MASK .OR. X$M_CW_BACK_PIXEL
C
        XSWDA.X$L_SWDA_EVENT_MASK = X$M_EXPOSURE.OR.X$M_BUTTON_PRESS
     &                              .OR.X$M_BUTTON_RELEASE
        XSWDA.X$L_SWDA_BACKGROUND_PIXEL =
     &      X$BLACK_PIXEL_OF_SCREEN(SCREEN)
C
        WIDTH = 400
        HEIGHT = 400
        WINDOW = X$CREATE_WINDOW(DPY,
     &           X$ROOT_WINDOW_OF_SCREEN(SCREEN),
     &           0, 0, WIDTH, HEIGHT, 0,
     &           DEPTH, X$C_INPUT_OUTPUT, VISUAL, ATTR_MASK, XSWDA)
C
        XSZHN.X$L_SZHN_X = WIDTH
        XSZHN.X$L_SZHN_Y = HEIGHT
        XSZHN.X$L_SZHN_WIDTH = WIDTH
        XSZHN.X$L_SZHN_HEIGHT = HEIGHT
        XSZHN.X$L_SZHN_FLAGS = X$M_P_POSITION .OR. X$M_P_SIZE
c
        CALL X$SET_NORMAL_HINTS(DPY, WINDOW, XSZHN)
c
        CALL X$STORE_NAME(DPY,WINDOW,'Lyapunov Fractal')
C
C       Map the window
C
        CALL X$MAP_WINDOW(DPY, WINDOW)
C
        MAP = X$DEFAULT_COLORMAP_OF_SCREEN(SCREEN)
C
        XGCVL.X$L_GCVL_BACKGROUND = X$BLACK_PIXEL_OF_SCREEN(SCREEN)
        GCVLF = X$M_GC_FOREGROUND .OR. X$M_GC_BACKGROUND
        XGCVL.X$L_GCVL_FOREGROUND = X$WHITE_PIXEL_OF_SCREEN(SCREEN)
        FLAGS = X$M_DO_RED.OR.X$M_DO_GREEN.OR.X$M_DO_BLUE
        GC = X$CREATE_GC(DPY,WINDOW,GCVLF,XGCVL)
C
C       Handle events
C
        IJK = 0
        DO WHILE (IJK.eq.0)
          CALL X$NEXT_EVENT(DPY, EVENT)
          if(event.evnt_type.eq.x$c_expose) then
c
c Pass control to the real graphics stuff
c
            call graphics(dpy,window,gc,width,height)
          else IF (EVENT.EVNT_TYPE .EQ. X$C_BUTTON_PRESS) THEN
            IF(EVENT.EVNT_BUTTON.X$L_BTEV_BUTTON.EQ.
     &           X$C_BUTTON3) THEN
C
C Button 3 pressed ...
C Unmap and destroy window
C
              CALL X$UNMAP_WINDOW(DPY, WINDOW)
              CALL X$DESTROY_WINDOW(DPY, WINDOW)
              CALL X$CLOSE_DISPLAY(DPY)
              CALL SYS$EXIT(%VAL(1))
            ENDIF
          END IF
        END DO

        END

subroutine graphics

      subroutine graphics(dpy,window,gc,wp,hp)
      parameter (maxcol=32)
      parameter (nf=12,niter=1000)
      integer ivmap(maxcol),wp,hp
      integer index(nf)
      integer index_iter(niter)
      integer*2 xpo(maxcol,50000)
      integer*2 ypo(maxcol,50000)
      integer npo(maxcol)
      real fac,b,t,step
      structure /point/
        integer*2 xpoint
        integer*2 ypoint
      end structure
      record /point/ p(50000)
c
      data b,t /3.3,4.0/
      data iseed /9999/
c
c generate a reandom index
c
      n_one = max(1,nint(nf*ran(iseed)))
      do 29 j=1,nf
         index(j) = 0
   29 continue
      do 30 j=1,n_one
         i_one = max(1,nint(nf*ran(iseed)))
   28    continue
         if (index(i_one).eq.1) then
            i_one = mod(i_one+1,nf+1)
            goto 28
         endif
         index(i_one) = 1
   30 continue
      write(6,'(a,15i1)') ' Generating with ',index
      do 31 i=1,niter
         index_iter(i) = index(mod(i-1,nf)+1)
   31 continue
c
      do 5 i=1,maxcol
         npo(i) = 0
    5 continue
      o = 1./log(2.)
      fac = (t-b)/real(wp-1)
      istep_size = 1
      write(6,*) ' '
      do 3 ir1=1,wp,istep_size
         r1 = b + fac*(ir1-1)
         rt1 = r1*2.
         write(6,'(a,i3,a,i4)') '+Step ',ir1,' of ',wp/istep_size
         do 4 ir2=1,hp,istep_size
            r2 = b + fac*(ir2-1)
            rt2 = r2*2.
            y = 0.51
            s = 0.0
            do 1 i=1,niter
               if(index_iter(i).eq.0) then
                 R = r1
                 RT = rt1
               else
                 R = r2
                 RT = rt2
               endif
               y = R*y*(1.-y)
               c = abs(R-RT*y)
               if(c.le.0.) goto 1
               s = s + log(c)
   1        continue
            if(s.gt.0.) goto 4
            s = o*s/real(niter)
            e = abs(s)*50.
            icol = max(min(maxcol,nint(e)),1)
            npo(icol) = npo(icol) + 1
            xpo(icol,npo(icol)) = ir1
            ypo(icol,npo(icol)) = ir2
   4     continue
   3  continue
      do 6 i=1,maxcol
         if(npo(i).eq.0) goto 6
c         write(6,*) ' colour ',i,npo(i),' points'
         do 7 j=1,npo(i)
            p(j).xpoint = xpo(i,j)
            p(j).ypoint = ypo(i,j)
   7     continue
         call x$set_foreground(dpy,gc,i+10)
         call x$draw_points(dpy,window,gc,p,npo(i),
     &                      x$c_coord_mode_origin)
   6  continue
      end