rem Lyapunov Exponent ....... Rev 6.4
      rem A J Tooth // 23rd December 2004
      rem Assembly language routines included April 2005
      on error if (err=17) then quit
      
      rem ============================================================
      rem This program demonstrates the chaotic pictures generated by
      rem estimating the Lyapunov Exponent e.g. for the function F=Q[P(x)]
      rem where P(x)=px(a-x) and Q(x)=qx(b-x), for a range of values
      rem for p and q between 1 and 4.
      rem After an article by Andy Burbanks, University of Cambridge
      rem ============================================================
      
      himem=lomem + 6000000
      *FLOAT64
      
      rem Setup
      proc_setup
      
      *REFRESH OFF
      rem Main Routine
      rp&=1
      repeat
        vdu 5
        gcol 3 : move 50,1500 :print;"Press -q- for Quick, -v- for Very Quick, or -s- for Slow output."
        
        if rp&=1 then
          gcol 2 : move 50,1450 :print;"After the picture is drawn, press -";:gcol 1:\
          \print;"m";:gcol 2:print;"- to make the mouse appear."
          move 50,1400 :print;"Choose an area to zoom into. First click the left button, move the mouse and click the right button."
          rp&=2
        endif
        vdu 4
        *REFRESH
        repeat
          c$=get$
        until (c$="q" or c$="s" or c$="v")
        case c$ of
          when "q": Stp%=2
          when "s": Stp%=1
          when "v": Stp%=4
        endcase
        cls
        *REFRESH
        trk&=0
        for !i%=0 to Xnum%-1 step Stp%
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &100
          
          for !j%=0 to Ynum%-1 step Stp%
            rem Set the parameters within the range specified
            p=ps+((pe-ps)*!i%/(Xnum%-1))
            q=qs+((qe-qs)*!j%/(Ynum%-1))
            
            rem Calculate an estimate for the Lyapunov Exponent (LE)
            lya=fn_Lyapunov(p,q)
            
            rem Set colours according to whether the LE is +ve or -ve
            proc_setcol(lya)
            rem Plot on the screen
            call drw%
          next
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &20
          
          trk&+=1
          if (trk&=5 or !i%>=Xnum%-1-Stp%) then
            rem Display BMP
            sys "SetStretchBltMode", @memhdc%, 3
            command$="MDISPLAY "+str$~pic%+" 0,0,"+str$(2*xscreen%)+","+str$(2*yscreen%)
            oscli command$
            *REFRESH
            trk&=0
          endif
        next
        *REFRESH
        a$=get$
        if a$="m" then proc_zoom
      until a$<>"m"
      *REFRESH ON
      quit
      end
      rem End of Program +++++++++++++++++++++++++++++++++++++++++++++
      
      rem Set up use of Full Screen
      def proc_fullscreen(return xscreen%, return yscreen%)
      
      sys "GetSystemMetrics", 0 to xscreen%
      sys "GetSystemMetrics", 1 to yscreen%
      sys "SetWindowLong",@hwnd%,-16,&16000000
      sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0
      vdu 23,22,xscreen%;yscreen%;8,16,16,1 : rem Set fullscreen mode
      mouse off : off : rem Turns off the Mouse Pointer and the Cursor
      endproc
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Setup
      def proc_setup
      
      rem Go to fullscreen
      proc_fullscreen(xscreen%,yscreen%)
      
      rem Initialise parameters
      dim t% 0, Lim% 0
      proc_init
      
      dim drw% 500, Red% 0, Gre% 0, Blu% 0, i% 3, j% 3
      dim fPQ% 100, fPQdash% 100
      dim pic% 2359350
      
      proc_bmpheader
      
      rem Dual-pass assembly, in case of labels
      for pass&=0 to 2 step 2
        proc_drw(pass&)
        proc_fPQ(pass&)
        proc_fPQdash(pass&)
      next pass&
      endproc
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      def fn_Lyapunov(pp,qq)
      local Lya,avg,x,xp1,Res%
      
      rem Traps Computational Range Errors ONLY
      on error local goto 2410
      avg=0.0 : x=x0
      for ?t%=1 to ?Lim%
        xp1=fn_F(x,pp,qq)
        avg+=fn_lnFdash(x,pp,qq)
        x=xp1
      next
      Lya=(avg/?Lim%)
      if Lya<0 then Res=100.0*Lya else Res=250.0*Lya
      goto 1310
 1300 Res=10000.0
 1310 =Res
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Update iterated value of the combined function
      def fn_F(X,P,Q)
      local Res,Y,Z,ZZ
      
      case d$ of
        when "1": Y=fn_P(X,P) : Res=fn_Q(Y,Q)
        when "2": Y=fn_Q(X,Q) : Res=fn_P(Y,P)
        when "3": Y=fn_P(X,P) : Z=fn_Q(Y,Q) : Res=fn_Q(Z,Q)
        when "4": Y=fn_P(X,P) : Z=fn_Q(Y,Q) : Res=fn_P(Z,Q)
        when "5": Y=fn_Q(X,Q) : Z=fn_P(Y,P) : ZZ=fn_P(Z,P) : Res=fn_Q(ZZ,Q)
      endcase
      =Res
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Calculate next value of lnF'(x)
      def fn_lnFdash(X,P,Q)
      local Res,Y,Yd,Z,Zd,ZZ,ZZd,ZZZ
      
      case d$ of
        when "1": Y=fn_P(X,P) : Yd=fn_Pdash(X,P) : Res=Yd*fn_Qdash(Y,Q)
        when "2": Y=fn_Q(X,Q) : Yd=fn_Qdash(X,Q) : Res=Yd*fn_Pdash(Y,P)
        when "3":
          Y=fn_P(X,P) : Yd=fn_Pdash(X,P)
          Zd=Yd*fn_Qdash(Y,Q) : Z=fn_Q(Y,Q) : Res=Zd*fn_Qdash(Z,Q)
        when "4":
          Y=fn_P(X,P) : Yd=fn_Pdash(X,P)
          Zd=Yd*fn_Qdash(Y,Q) : Z=fn_Q(Y,Q) : Res=Zd*fn_Pdash(Z,Q)
        when "5":
          Y=fn_Q(X,Q) : Yd=fn_Qdash(X,Q) : Zd=Yd*fn_Pdash(Y,P)
          ZZ=fn_P(Y,P): ZZd=fn_Pdash(ZZ,P): ZZZ=fn_P(ZZ,P)
          Res=Zd*ZZd*fn_Qdash(ZZZ,Q)
      endcase
      if abs(Res)=<0.0000001 then Res=0.0000001
      =ln(abs(Res))
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Individual function elements
      rem ============================
      def fn_P(X,P)
      
      call fPQ%,X,P,A,Ans
      =Ans
      
      def fn_Q(X,Q)
      
      call fPQ%,X,Q,B,Ans
      =Ans
      
      def fn_Pdash(X,P)
      
      call fPQdash%,X,P,A,Ans
      =Ans
      
      def fn_Qdash(X,Q)
      
      call fPQdash%,X,Q,B,Ans
      =Ans
      rem ============================
      
      rem Enables user to zoom into a particular area
      def proc_zoom
      local x%,y%,b%,ax%,bx%,ay%,bi%,rnl,inl
      
      *REFRESH ON
      mouse on : gcol 7
      rem Select first corner of rectangular zoom area
      repeat
        mouse x%,y%,b%
        ax%=x% : ay%=y%
      until b%=4
      rem Print a "+" at the first corner
      move ax%-10,ay% : draw ax%+10,ay%
      move ax%,ay%-10 : draw ax%,ay%+10
      rem Select opposite corner of zoom area
      repeat
        mouse x%,y%,b%
        bx%=x% : bi%=y%
      until b%=1
      rem Draw a box and wait 1 second before moving on
      move ax%,ay% : draw bx%,ay% : draw bx%,bi%
      draw ax%,bi% : draw ax%,ay%
      b$=inkey$(100)
      
      rem Re-assign corner points in ascending order
      if ax%<bx% then
        lx%=ax% : ux%=bx%
      else
        lx%=bx% : ux%=ax%
      endif
      if ay%<bi% then
        ly%=ay% : uy%=bi%
      else
        ly%=bi% : uy%=ay%
      endif
      rem Re-assign viewing window limits
      pns=ps+((pe-ps)*(lx%/(2*xscreen%-1)))
      pne=ps+((pe-ps)*(ux%/(2*xscreen%-1)))
      
      qns=qs+((qe-qs)*(ly%/(2*yscreen%-1)))
      qne=qs+((qe-qs)*(uy%/(2*yscreen%-1)))
      ps=pns : pe=pne : qs=qns : qe=qne
      mouse off
      *REFRESH OFF
      endproc
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Traps Computational Range Errors ONLY
 2410 if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then goto 1300
      else restore error
      endif
      goto 1300
      rem-----------------------------------------------------------------
      
      rem Displays a backdrop
      def proc_back(Xs%,Ys%,Ws%,Hs%,Rr%,Gg%,Bb%)
      local h%,Rf%,Gf%,Bf%
      
      for h%=0 to 30
      Rf%=Rr%*h%/30 : Gf%=Gg%*h%/30 : Bf%=Bb%*h%/30
      colour 9,Rf%,Rf%,Rf% : gcol 9
      rectangle fill Xs%+h%,Ys%+h%,Ws%-2*h%,Hs%-2*h%
      next h%
      
      endproc
      rem-----------------------------------------------------------------
      
      rem Initialise parameters
      def proc_init
      
      colour 12,100,100,100
      Xnum%=1024 : Ynum%=768
      
      proc_AJTicon(10,700)
      *FONT Desdemona,18
      proc_back(400,int(1200*yscreen%/768),1000,int(200*yscreen%/768),130,130,130)
      vdu 5
      gcol 13 : move 600,int(1330*yscreen%/768) : print;"The Lyapunov Exponent"
      vdu 4
      *FONT Arial,12,b
      ps=1.0 : pe=4.0
      qs=1.0 : qe=4.0
      colour 3
      print tab(10,15);"Choose Method 1, 2, 3, 4 or 5"
      repeat
      d$=get$
      until (d$="1" or d$="2" or d$="3" or d$="4" or d$="5")
      colour 1 : print tab(40,15);"Method ";d$
      colour 3
      print tab(10,17);"Choose values for x between 0 and 1 (0.5)"
      colour 1: input tab(50,17);x0
      if (x0<=0 or x0>=1) then x0=0.5
      colour 3
      print tab(10,19);"Choose values for 2 of the parameters: A and B"
      print tab(10,20);"Enter a value for A in the range 0.8 to 1.5"
      colour 1: input tab(50,20);A
      if A=0 then A=1.0
      A=1.0*A
      colour 3: print tab(10,22);"Enter a value for B in the range 0.8 to 1.5"
      colour 1: input tab(50,22);B
      if B=0 then B=1.0
      B=1.0*B
      colour 3: print tab(10,24);"Press -l- for 5 iterations or -h- for 10 "
      repeat
      it$=get$
      until (it$="l" or it$="h")
      case it$ of
      when "l": ?Lim%=5
      when "h": ?Lim%=10
      endcase
      colour 1: print tab(65,24);"Iterations = "?Lim%
      
      colour 3: print tab(10,28);"Choose colour scheme Light (l), Dark (d), Metallic (m), Gothic (g) or Wispy (w)"
      repeat
      cl$=get$
      until (cl$="l" or cl$="d" or cl$="m" or cl$="g" or cl$="w")
      case cl$ of
      when "l": Msg$="Light" : cols&=1
      when "d": Msg$="Dark"  : cols&=2
      when "m": Msg$="Metallic" : cols&=3
      when "g": Msg$="Gothic" : cols&=4
      when "w": Msg$="Wispy" : cols&=5
      endcase
      colour 1: print tab(90,28);Msg$
      if (cols&=3 or cols&=4) then
      colour 3: print tab(10,30);"Copper (c), Steel (s) or Gold (g)? "
      repeat
        m$=get$
      until (m$="c" or m$="s" or m$="g")
      case m$ of
        when "c": Msg$="Copper" : r1%=80 : g1%=5  : b1%=0  : r2=170.0 : g2=200.0 : b2=175.0
        when "s": Msg$="Steel"  : r1%=10 : g1%=10 : b1%=30 : r2=240.0 : g2=240.0 : b2=220.0
        when "g": Msg$="Gold"   : r1%=90 : g1%=60 : b1%=15 : r2=155.0 : g2=165.0 : b2=125.0
      endcase
      colour 1: print tab(60,30);Msg$
      endif
      if cols&=1 then
      colour 3: print tab(10,32);"What rate of initial colour change Low (l), Medium (m) High (h) or Very High (v)?"
      repeat
        m$=get$
      until (m$="l" or m$="m" or m$="h" or m$="v")
      case m$ of
        when "l": Msg$="Low"    : Cut=50.0
        when "m": Msg$="Medium" : Cut=30.0
        when "h": Msg$="High"   : Cut=10.0
        when "v": Msg$="Very High"  : Cut=1.0
      endcase
      colour 1: print tab(85,32);Msg$
      endif
      colour 2
      print tab(10,34);"Once each screen has finished, press -m- to activate the mouse."
      print tab(10,35);"This will enable you to zoom in on a section of the picture."
      colour 12 : print tab(10,37);"Any area in SOLID GREY is 'out of area'"
      a$=inkey$(100)
      cls
      endproc
      rem-----------------------------------------------------------------
      
      rem Set colours according to whether the LE is +ve or -ve
      def proc_setcol(Lya)
      
      if Lya=10000.0 then
      ?Red%=100 : ?Gre%=100 : ?Blu%=100
      else
      ?Red%=0 : ?Blu%=0 : ?Gre%=0
      if cols&>2 then
        if Lya<=0.0 then
          Ly=2.0*ln(1-Lya)
          case cols& of
            when 4: fac=(1-cos(Ly))/2 : ?Red%=r1%+int(r2*fac) : ?Gre%=g1%+int(g2*fac) : ?Blu%=b1%+int(b2*fac)
            when 3: fac=(1+sin(Ly))/2 : ?Red%=r1%+int(r2*fac) : ?Gre%=g1%+int(g2*fac) : ?Blu%=b1%+int(b2*fac)
            when 5: fac=((1-cos(0.17*Ly))/1.7)^4 : if fac>1.0 then fac=1.0/fac
              ?Red%=int(254.0*fac) : ?Gre%=?Red% : ?Blu%=?Red%
            otherwise
              rem Do nothing
          endcase
        else
          Ly=0.75*ln(1+Lya)
          fac=(1.0+0.98*sin(Ly))/2
          if fac<0.0 then fac=0.0
          if fac>1.0 then fac=1.0
          ?Blu%=int(254.0*fac)
        endif
      else
        if Lya<=-255.0 and Lya>-510.0 then
          if cols&=1 then
            lam=1.0+((Lya+255.0)/255.0)
            ?Red%=int(180.0*lam)
            ?Blu%=0
          else
            Mid=-255.0-Lya : if Mid>255.0 then Mid=255.0
            ?Red%=int(Mid)
            ?Gre%=127
          endif
        endif
        if Lya<=-510.0 then
          if cols&=2 then
            ?Red%=255
            DMid=-Lya-510.0 : if DMid>255.0 then DMid=255.0
            ?Gre%=127+int(DMid/2)
          else
            DMid=-Lya-510.0 : if DMid>255.0 then DMid=255.0
            ?Gre%=int(DMid)
          endif
        endif
        
        if (Lya>=-255.0 and Lya<=0.0) then
          if cols&=1 then
            lam=-Lya/255.0 : mu=1-lam : Mid=-Lya : if Mid>Cut then Mid=Cut
            gam=Mid/Cut
            ?Red%=int(gam*(250.0*mu + 180.0*lam) + (1-gam)*80.0)
            ?Gre%=int(gam*(235.0*mu) + (1-gam)*30.0)
            ?Blu%=int((1-gam)*20.0)
          else
            ?Gre%=-int(Lya/2)
          endif
        endif
        
        if Lya>0.0 then
          Mid=Lya
          if Mid>255.0 then
            ?Blu%=255
            DMid=Mid-255.0 : if DMid>255.0 then DMid=255.0
            ?Gre%=int(DMid) : ?Red%=int(DMid)
          else
            ?Blu%=int(Mid)
          endif
        endif
      endif
      endif
      endproc
      rem-----------------------------------------------------------------
      
      rem Assembly Language Routine
      rem for Plotting
      def proc_drw(opt&)
      
      P%=drw%
      [opt opt&
      
      mov ebx,[i%]
      cmp ebx,0
      jl miss
      cmp ebx,1023
      jg miss
      mov eax,[j%]
      cmp eax,0
      jl miss
      cmp eax,767
      jg miss
      
      imul eax,3072
      add eax,[i%]
      add eax,[i%]
      add eax,[i%]
      add eax,54
      
      mov cl,[Blu%]
      mov pic%[eax],cl
      mov cl,[Gre%]
      mov pic%[eax+1],cl
      mov cl,[Red%]
      mov pic%[eax+2],cl
      
      .miss
      ret
      ]
      endproc
      rem-----------------------------------------------------------------
      
      rem Displays my Icon in .exe version
      def proc_AJTicon(i%,j%)
      
      sys "GetModuleHandle", 0 to hm%
      sys "LoadImage", hm%, "BBCWrun", 1, 32, 32, 0 to hicon%
      w% = 32
      h% = 32
      sys "DrawIconEx", @memhdc%, i%, j%, hicon%, w%, h%, 0, 0, 3
      sys "InvalidateRect", @hwnd%, 0, 0
      endproc
      rem-----------------------------------------------------------------
      
      rem Assembly Language Routines
      rem for Functions
      def proc_fPQ(opt&)
      
      P%=fPQ%
      [opt opt&
      
      finit
      mov esi,[ebp+12]
      fld qword [esi]
      mov esi,[ebp+2]
      fld qword [esi]
      fsub st1,st0
      fmulp st1,st0
      mov esi,[ebp+7]
      fld qword [esi]
      fmulp st1,st0
      mov esi,[ebp+17]
      fstp qword [esi]
      
      ret
      ]
      endproc
      rem-----------------------------------------------------------------
      
      def fn_dumP(X,P)
      =P*X*(A-X)
      def fn_dumQ(X,Q)
      =Q*X*(B-X)
      rem-----------------------------------------------------------------
      rem Assembly Language Routines
      rem for Functions
      def proc_fPQdash(opt&)
      
      P%=fPQdash%
      [opt opt&
      
      finit
      mov esi,[ebp+12]
      fld qword [esi]
      mov esi,[ebp+7]
      fld qword [esi]
      fmul st1,st0
      mov esi,[ebp+2]
      fld qword [esi]
      fmulp st1,st0
      fadd st0,st0
      fsubp st1,st0
      mov esi,[ebp+17]
      fstp qword [esi]
      
      ret
      ]
      endproc
      rem-----------------------------------------------------------------
      def fn_dumPdash(X,P)
      =A*P-2*P*X
      
      def fn_dumQdash(X,Q)
      =B*Q-2*Q*X
      rem ============================
      
      rem Standard 1024x768 Header setup
      def proc_bmpheader
      local s&, t%, a&
      
      restore
      for s&=0 to 53
      read a&
      ?(pic%+s&)=a&
      next s&
      
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      rem Setup blank bmp
      for t%=0 to 786431
      !(pic%+3*t%+54)=0
      next t%
      endproc
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      
      rem Header data
      data 66,77,54,0,36,0,0,0,0,0,54,0,0,0,40,0,0,0
      data 0,4,0,0,0,3,0,0,1,0,24,0,0,0,0,0,0,0,36,0
      data 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
      rem =======================================================================