rem Catacaustic Curves ........ Rev 3.2
      rem A J Tooth // 29th December 2003
      on error if (err=17) then quit
      
      rem ===============================================================
      rem Displays Catacaustic curves. These are curves formed when light
      rem reflects from the surface of the given curve.
      rem ===============================================================
      
      rem Preamble
      *FLOAT64
      proc_fullscreen
      proc_theta : rem Defines a bit pattern for a "theta" character
      xmax%=1024 : ymax%=768 : rmax%=ymax%
      
      Ctr&=1 : rem Controls where to return after computational errors
      
      *FONT Verdana,10,B
      
      rem Initial Choices
      proc_init
      
      repeat
        rem Draw the Axes
        proc_axes
        
        rem Identify Preset Curve details
        proc_iden
        
        rem Draw the Curve
        proc_curv
        
        rem Enables selection of Light Source using the mouse.
        proc_light
        
        rem Draw the Catacaustic Curve
        proc_rays
        
        rem Repeat Section
        mouse on 3
        Flg%=0 : a$=""
        repeat
          mouse j%,k%,b%
          if (abs(j%-px%)>15) or (abs(k%-py%)>15) then Flg%=1
          a$=inkey$(10)
        until (a$<>"" or inkey(-10) or Flg%=1)
        cls
      until (a$<>" " and a$<>"")
      
      if a$="r" or a$="R" then run
      
      quit
      end
      rem End of Program ++++++++++++++++++++++++++++++++++++++++++++++++
      rem ===============================================================
      
      rem Set up use of Full Screen
      def proc_fullscreen
      
      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 Sets up Axes
      def proc_axes
      
      if Ctr&=1 then
        colour 3 : print tab(10,35);" Choose Centred (c) or Positive Half-Plane (h) or 1st Quadrant only (q) "
        ch$=get$ : cls
      endif
      case ch$ of
        when "c","C" : ffx%=0 : ffy%=0
        when "h","H" : ffx%=1000 : ffy%=0
        when "q","Q" : ffx%=1000 : ffy%=750
        otherwise
          ffx%=0 : ffy%=0
      endcase
      
      setx%=(xmax%-ffx%) : sety%=(ymax%-ffy%)
      origin setx%,sety% : gcol 3
      move (-xmax%+ffx%),0 : draw (xmax%+ffx%),0
      move 0,(-ymax%+ffy%) : draw 0,(ymax%+ffy%)
      endproc
      rem ===============================================================
      
      rem Traps Computational Range Errors ONLY
  870 if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then
        on Ctr& goto 2630,3230,3900
      else restore error
      endif
      on Ctr& goto 2630,3230,3900
      rem ===============================================================
      
      rem Main Calculation
      def fn_func(R$,a,t)
      
      rem R$ MUST be expressed in terms of t and a !!!
      =eval(R$)
      rem ===============================================================
      
      rem Menu of Preset Curves
      def proc_menu
      
      colour 2: print tab(5,10);"Pick one of the curves from the list below, by pressing the relevant number."
      
      colour 3 : print tab(10,15);"a ";:colour 7:print;"CIRCLE "
      colour 3 : print tab(10,17);"b ";:colour 7:print;"PARABOLA"
      colour 3 : print tab(10,19);"c ";:colour 7:print;"CARDIOID"
      colour 3 : print tab(10,21);"d ";:colour 7:print;"DOUBLE FOLIUM"
      colour 3 : print tab(10,23);"e ";:colour 7:print;"COCHLEOID"
      colour 3 : print tab(10,25);"f ";:colour 7:print;"LEMNISCATE of Bernoulli"
      colour 3 : print tab(10,27);"g ";:colour 7:print;"LIMACON of Pascal"
      colour 3 : print tab(10,29);"h ";:colour 7:print;"LITUUS"
      colour 3 : print tab(10,31);"i ";:colour 7:print;"RIGHT STROPHOID"
      colour 3 : print tab(10,33);"j ";:colour 7:print;"EXPERIMENTAL"
      
      repeat
        as$=get$
        as&=asc(as$)-96
      until (as&>0 and as&<11)
      cls
      case as& of
        when 1 : Curv$=" The CIRCLE" : proc_Msg(Curv$)
          r$="a"
        when 2 : Curv$="The PARABOLA" : proc_Msg(Curv$)
          r$="(a*COS(t))/(SIN(t)*SIN(t))"
        when 3 : Curv$="The CARDIOID" : proc_Msg(Curv$)
          r$="a*(1+COS(t))/2"
        when 4 : Curv$="The DOUBLE FOLIUM" : proc_Msg(Curv$)
          r$="3*a*COS(t)*SIN(t)*SIN(t)"
        when 5 : Curv$="The COCHLEOID" : proc_Msg(Curv$)
          r$="(a*SIN(t))/t"
        when 6 : Curv$="The LEMNISCATE of Bernoulli" : proc_Msg(Curv$)
          r$="a*SQR(COS(2*t))"
        when 7 : Curv$="The LIMACON of Pascal" : proc_Msg(Curv$)
          r$="(1+(2*a*COS(t)))/2"
        when 8 : Curv$="The LITUUS" : proc_Msg(Curv$)
          r$="a/(SQR(t))"
        when 9 : Curv$="The RIGHT STROPHOID" : proc_Msg(Curv$)
          r$="a*COS(2*t)/COS(t)"
        when 10 : Curv$="EXPERIMENTAL" : proc_Msg(Curv$)
          r$="a*(1-SIN(2*t))/COS(t)"
          
      endcase
      proc_defpars
      endproc
      rem ==================================================================================================
      
      rem Display the chosen curve
      def proc_Msg(M$)
      
      colour 7 : print tab(5,10);"You have chosen:   ";:colour 2:print;M$
      endproc
      rem ==================================================================================================
      
      rem Default to self-entry of Expression and Parameters
      def proc_default
      
      print " Input the Expression for R in terms of the angle, t, and a parameter, a"
      print:input " Enter the Expression here:  ",r$
      
      print:input " Enter value for the parameter, a ",Pe
      if Pe=0 then Pe=1
      
      print:print " What end-points (in degrees) do you want for the Angle, Theta ? "
      print:input " Theta Start (0)",Ths," Theta End (360)",The
      if The=0 then The=360
      
      print:input " Axis Scale (1) ",hnv
      if hnv=0 then hnv=1
      print:input " How many plot points (2000)",Pts%
      if Pts%=0 then Pts%=2000
      
      print:input " Enter a NAME for the curve, if required, otherwise just press -Enter-.",Curv$
      endproc
      rem ==================================================================================================
      
      rem Default values of parameters
      def proc_defpars
      
      Pe=1
      Ths=0 : The=360
      rem 3 times round if the LITUUS was chosen; more interesting
      if as&=8 then The=1080
      Pts%=5000
      colour 7:print tab(10,17);" Axis Scale (1) ";:colour 1:input;hnv
      if hnv=0 then hnv=1
      
      endproc
      rem ==================================================================================================
      
      rem Defines ASCII character 129 as a symbol for Theta
      def proc_theta
      local r1&,r2&,r3&,r4&,r5&,r6&,r7&,r8&
      
      r1&= %00000000
      r2&= %00111100
      r3&= %01100110
      r4&= %01111110
      r5&= %01100110
      r6&= %00111100
      r7&= %00000000
      r8&= %00000000
      
      vdu 23,129,r1&,r2&,r3&,r4&,r5&,r6&,r7&,r8&
      endproc
      rem ==================================================================================================
      
      rem Identify Preset Curve details
      def proc_iden
      
      *FONT Desdemona,16,B
      colour 12: print tab(2,1);Curv$
      *FONT Arial,10,B
      l%=len(r$)
      pr$=""
      colour 3
      for z%=1 to l%
        tr$=mid$(r$,z%,1)
        case tr$ of
          when "t" : pr$=pr$+chr$(129)
          when "*" : pr$=pr$+"."
          otherwise
            pr$=pr$+tr$
        endcase
      next z%
      
      print tab(2,4);"Polar Form: r=";pr$
      endproc
      rem ==================================================================================================
      
      rem Select Light Source with the mouse
      def proc_light
      
      mouse on 3 : gcol 7
      red=hnv*rmax%
      repeat
        mouse px%,py%,m%
        px=px% : py=py%
        Px=px/red : Py=py/red
      until m%=4
      mouse off
      endproc
      rem ==================================================================================================
      
      rem Displays the Chosen curve
      def proc_curv
      local t%
      
      vdu 19,2,16,0,255,0 : gcol 2
      Ctr&=1
      for t%=0 to Pts%
        Th=rad(Ths + ((The-Ths)*t%/Pts%))
        
        rem Computational Error Trapping
        on error local goto 870
        
        rem Main calculation
        r=(hnv*rmax%)*fn_func(r$,Pe,Th)
        x=r*cos(Th) : y=r*sin(Th)
        
        rem Plot formatting section
        sy%=sgn(y) : sx%=sgn(x)
        x%=sx%*int(abs(x))
        y%=sy%*int(abs(y))
        
        if (abs(y%)<=1600) or (abs(x%)<=2100) then plot 69,x%,y%
        
 2630 next t%
      endproc
      rem ==================================================================================================
      
      rem Draws the rays of light
      def proc_rays
      
      vdu 19,4,16,100,0,0, : gcol 4
      Ctr&=2 : mxo%=10000 : myo%=10000
      for s%=0 to Ray%
        
        Th=rad(Ths + ((The-Ths)*s%/Ray%))
        rem Computational Error Trapping
        on error local goto 870
        
        r=(hnv*rmax%)*fn_func(r$,Pe,Th)
        x=r*cos(Th) : y=r*sin(Th)
        
        rem Plot formatting section
        sy%=sgn(y) : sx%=sgn(x)
        x%=sx%*int(abs(x))
        y%=sy%*int(abs(y))
        
        rem Only plot points if they are actually within the screen area
        if (abs(y%)<=3200) or (abs(x%)<=4200) then
          
          if Ray%>300 then
            move x%,y%
          else
            move px%,py% : draw x%,y%
          endif
          
          rem Calculate estimate for dr/dtheta
          Thdel=Th+0.01
          r1=(hnv*rmax%)*fn_func(r$,Pe,Thdel)
          drdt=(r1-r)/0.01
          
          rem Calculate outward normal to the curve
          den=sqr(drdt*drdt + r*r) : alp=drdt/den : bet=r/den
          nx=bet*cos(Th) + alp*sin(Th) : ny=bet*sin(Th) - alp*cos(Th)
          qdotn=((x-px)*nx) + ((y-py)*ny)
          
          rem Calculate vector for reflected ray of light
          vx=(x-px) - (2*qdotn*nx) : vy=(y-py) - (2*qdotn*ny)
          
          sy%=sgn(vy) : sx%=sgn(vx)
          vx%=sx%*int(abs(vx))
          vy%=sy%*int(abs(vy))
          if (abs(vy%)<=3200) or (abs(vx%)<=4200) then
            
            rem Do not draw individual rays for 1000 rays or more; no detail visible anyway
            if Ray%<1000 then draw (x%+5*vx%),(y%+5*vy%)
            
          endif
          
          rem Draw the highlighted Catacaustic in WHITE if required
          if (cat$="y" or cat$="Y") then
            proc_hicaust(s%)
          endif
          
        endif
 3230 next s%
      endproc
      rem ==================================================================================================
      
      rem Initial Choices on Running
      def proc_init
      
      colour 7:print tab(5,10);" Press -";:colour 1:print;"m";:colour 7:\
      \print;"- to pick a curve from a pre-set list, or any other key to enter your own."
      m$=get$
      case m$ of
        when "m", "M" : cls:proc_menu
        otherwise
          cls:proc_default
      endcase
      
      rem Choose number of Rays
      colour 2
      print tab(5,22);"Behaviour changes as you choose more rays !!!! "
      print tab(10,24);"* Less than 300, all rays are plotted."
      print tab(10,26);"* Less than 1000, only reflected rays are plotted."
      print tab(10,28);"* More than 1000, ONLY the Catacaustic itself is shown."
      colour 3
      print tab(5,20);" How many light rays (150)",;:colour 1:input;Ray%
      if Ray%=0 then Ray%=150
      
      rem Wanna see the Catacaustic picked out?
      if Ray%<1000 then
        colour 3 : print tab(5,30);"Do you want to see a white outline of the Catacaustic (y/n) ";
        repeat
          cat$=get$ : colour 1 : print cat$ : colour 3
        until (cat$="y" or cat$="Y" or cat$="n" or cat$="N")
      else
        cat$="y"
      endif
      endproc
      rem ==================================================================================================
      
      rem Highlight the Catacaustic
      def proc_hicaust(Ss%)
      
      Ctr&=3 : on error local goto 870
      
      if Ss%>0 then
        det=vxo*vy-vx*vyo : rem Determinant for meeting of lines
        if det<>0.0 then
          mu=((x-xo)*vyo-(y-yo)*vxo)/det
          lam=((x-xo)*vy-(y-yo)*vx)/det
          
          rem Lines only meet if both mu and lam are greater than zero
          if (mu>0 and lam>0) then
            mx=x+mu*vx : my=y+mu*vy
            sx%=sgn(mx) : sy%=sgn(my)
            mx%=sx%*int(abs(mx))
            my%=sy%*int(abs(my))
            if (abs(my%)<=1536) or (abs(mx%)<=2048) then
              gcol 15
              if mxo%=10000 and myo%=10000 then
                plot 69,mx%,my%
              else
                if (abs(mx%-mxo%)<100) and (abs(my%-myo%)<100) then
                  move mxo%,myo% : draw mx%,my%
                endif
              endif
              gcol 4
              mxo%=mx% : myo%=my%
            endif
          endif
        endif
      endif
 3900 xo=x : yo=y : vxo=vx : vyo=vy
      Ctr&=2
      endproc
      rem ==================================================================================================