!  There is an apparently unavoidable discontinuity in the variation
!  of the eigenlines: if you sit on the parabola in the defective case
!  and push s across s = 0, the eigendirection vanishes and then
!  reappears in a perpendicular direction. I think this is just the
!  way it is. The vector field itself is quite continuous.

!  Idea for another day: It's always quite annoying that the shapes
!  change so quickly as you near the critical parabola. I wonder if this
!  could be alleviated using the microtuning we did on another tool,
!  using arrow keys on the tr and det sliders. Using the sliders overcomes
!  the problem of sliding up the parabola, so I think it could be done. 

!  - Haynes


!! File: LinearPhaseCursor
!! July 12, 2003 Hubert Hohn

PUBLIC PCFlag,Mac5Flag,M68KFlag,UnixFlag,xmax,ymax
PUBLIC toolLft,toolRgt,toolBas,toolTop,toolhdr,toolHgt,toolWid  ! tool boundaries
PUBLIC winLft,winRgt,winBas,winTop,winHgt,winWid      ! window
PUBLIC workLft,workRgt,workBas,workTop,workMidx  ! work area
PUBLIC backclr,black,drkgry,drkmid,midgry,litmid,litgry,white
PUBLIC red,yellow,green,cyan,blue,magenta,pink,colorscheme
PUBLIC planeclr,gridclr,rimclr,axisclr,axislabelclr,titleclr,rightsclr
PUBLIC numberlineclr,slotdrkclr,slotlgtclr,slideclr
PUBLIC largefonts, title$, SLUmode

LET toolHgt= 560
LET toolWid= 780
LET window$= "The d'Arbeloff Interactive Math Project"
LET colorscheme= 0

LET title$ = "Linear Phase Portraits: Cursor Entry"

SUB ThisProgram
    CALL LinearParameters
    CLEAR
END SUB


!! ---------------------------------------------------------
!! ------ Start TB4 Mac Header and Subs ------

!LET M68KFlag = 1
!LIBRARY "MacTools*", "HHLib.trc"
!ASK PIXELS winWid,winHgt
!LET winLft= 0
!LET winTop= 0
!LET winRgt= winWid-1
!LET winBas= winHgt-1
!SET WINDOW 0,winRgt,winBas,0
!
!CALL Palette
!CLEAR
!
!CALL ToolPanel
!CALL ThisProgram
!
!END
!EXTERNAL
!
!MODULE Mac4Parts
!    SUB SetTextFont(font,size,style$)
!        CALL MacTextFont(font)
!        CALL MacTextSize(size)
!        CALL MacTextFace(style$)
!    END SUB
!
!    SUB StringWidth(sw$,sl)
!        DECLARE DEF MacStringWidth
!        LET sl= MacStringWidth(sw$)
!    END SUB
!
!    SUB SetLineWeight(wgt)
!        CALL MacPenSize(wgt,wgt)
!    END SUB
!
!    SUB BoxDisk(Lft,Rgt,Bas,Top)
!        CALL MacPaintOval(Lft,Rgt,Bas,Top)
!    END SUB
!END MODULE

!! --- End TB4 Mac Header and Subs ---
!!---------------------------------------------------------
!!--- Start TB5 Cross-Platform header and subs ---

LIBRARY "c:\TB Gold 51a\TBLibs\TrueCtrl.trc"     ! windows
LIBRARY "c:\TB Gold 51a\TBLibs\HHLib.trc"

!LIBRARY ":TBLibs:TrueCtrl.trc"    ! macintosh
!LIBRARY "HHLib.trc"

PUBLIC WinID
DECLARE PUBLIC OBJM_SET,OBJM_SYSINFO
LET winHgt= toolHgt
LET winWid= toolWid
DIM values(1)

CALL TC_Init

CALL Object(OBJM_SYSINFO,WinID,"MACHINE",system$,values())
IF system$="MAC" then
   LET Mac5Flag= 1
ELSE IF system$="WIN32" then
   LET PCFlag  = 1
END IF

CALL TC_SetUnitsToPixels          ! 5.1 and up needs this
CALL TC_GetScreenSize(scrnLft,scrnRgt,scrnBas,scrnTop)
LET winLft= int((scrnRgt-scrnLft-winWid)/2)
LET winRgt= winLft+winWid-1
LET winTop= int((scrnBas-scrnTop-winHgt)/2) + 10
LET winBas= winTop+winHgt-1

CALL TC_Win_Create (WinID,"TITLE",winLft,winRgt,winBas,winTop)
LET values(1)= 2
CALL Object(OBJM_SET, WinID, "TYPE", "", values())

IF PCFlag=1 then                  ! kill dithering
   LET  values(1)= 1
   CALL Object(OBJM_SET, WinID, "SOLID MIX", "", values())
END IF

LET  values(1)= 0
CALL TC_SetRect(WinID,winLft,winRgt,winBas,winTop)
CALL TC_Win_SetTitle(WinID,window$)
CALL TC_Show(WinID)

SET MODE "COLORSTANDARD"

ASK PIXELS winWid,winHgt          ! must follow set mode
LET winLft= 0
LET winTop= 0
LET winRgt= winWid-1
LET winBas= winHgt-1
SET WINDOW 0,winRgt,winBas,0

CALL Palette

IF PCFlag=1 then
   LET  values(1)= 0              ! now force solid colors
   CALL Object(OBJM_SET, WinID, "SOLID MIX", "", values())
   CALL TC_Win_RealizePalette(WinID)   ! some PCs need this

   CALL TC_Win_SetFont(WinID,"arial",9,"plain")
   CALL StringWidth("0",sw)
   IF sw>7 then LET largefonts=1 else LET largefonts=0
END IF

CALL TC_Win_Switch(WinID)

CALL ToolPanel
CALL ThisProgram

CALL SetTextFont(1,12,"bold")     ! now shut down and clean up
LET quit$= "click the mouse or press a key to close..."
CALL PlotTextCJ(workmidx,(workbas+worktop)/2,quit$,yellow)
CALL TC_CleanUp

END
EXTERNAL

MODULE TB5Parts
    SUB StringWidth(sw$,sl)
        DECLARE PUBLIC WinID
        LET sl= StrWidth(WinID,sw$)
    END SUB

    SUB SetLineWeight(wgt)
        DECLARE PUBLIC OBJM_SET
        DECLARE PUBLIC WinID
        DIM values(1)
        LET values(1)= wgt
        CALL Object(OBJM_SET,WinID, "WIDTH", "", values())
    END SUB

    SUB SetTextFont(font,size,style$)
        DECLARE PUBLIC WinID,Mac5Flag,PCFlag,largefonts
        IF Mac5Flag=1 then
           SELECT CASE font
           CASE 4
                LET font$= "Courier"
           CASE 16
                LET font$= "Times"
           CASE else
                LET font$= "Geneva"
           END SELECT
        ELSE IF PCFlag=1 then
           IF largefonts=1 then
              IF size<12 then
                 LET size=  6
              ELSE IF size=14 then
                 LET size= 10
              ELSE IF size=18 then
                 LET size= 12
              ELSE IF size=24 then
                 LET size= 14
              ELSE IF size=12 then
                 LET size=  8
              ELSE
                 LET size= round(72/96 * size * .8)
              END IF
           ELSE
              IF size<12 then
                 LET size=  7
              ELSE IF size=14 then
                 LET size= 12
              ELSE IF size=12 then
                 LET size=  9
              ELSE IF size=18 then
                 LET size= 14
              ELSE IF size=24 then
                 LET size= 18
              ELSE
                 LET size= round(72/96 * size)
              END IF
           END IF
           SELECT CASE font
           CASE 4
                LET font$= "Courier New"
           CASE 16
                LET font$= "Times New Roman"
           CASE else
                LET font$= "Verdana"
           END SELECT
        END IF
        IF style$= "normal" then LET style$= "plain"
        CALL TC_Win_SetFont(WinID,font$,size,style$)
    END SUB

    SUB BoxDisk(Lft,Rgt,Bas,Top)
        BOX DISK Lft,Rgt,Bas,Top
    END SUB
END MODULE

!  --- End TB5 Cross-platform header and subs ---
!! ---------------------------------------------------------

!! --- Start Unix Header and Subs ---

!library "HHLib.unix"
!LET UnixFlag= 1
!ASK PIXELS winWid,winHgt
!LET winLft= 0
!LET winTop= 0
!LET winRgt= winWid-1
!LET winBas= winHgt-1
!SET WINDOW 0,winRgt,winBas,0
!CALL Palette
!
!CALL ToolPanel
!CALL ThisProgram
!
!END
!EXTERNAL
!
!MODULE UnixParts
!    SHARE CharWidth
!
!    SUB SetTextFont(font,size,style$)
!        LET font$= "-adobe-courier-"
!        IF style$= "normal" then
!           LET style$= "medium-r-normal--"
!        ELSE
!           LET style$= "bold-r-normal--"
!        END IF
!        IF size=9 then
!           LET size$= str$(10)
!        ELSE
!           LET size$= str$(size)
!        END IF
!        LET test= SetFont(font$&style$&size$&"*")
!
!        IF size=9 then
!           LET CharWidth=  6
!        ELSE IF size=12 then  ! numeric output - axis labels
!           LET CharWidth=  7
!        ELSE IF size=14 then  ! rare
!           LET CharWidth=  8
!        ELSE IF size=18 then  ! rare
!           LET CharWidth= 10
!        END IF
!    END SUB
!
!    SUB StringWidth(sw$,sl)                      ! string width in pixels
!        ! LET sl= StrWidth(sw$)
!        LET chars= len(sw$)
!        LET sl   = chars*CharWidth
!    END SUB
!
!    SUB SetLineWeight(wgt)
!        ! CALL PenSize(wgt,wgt)
!    END SUB
!
!    SUB BoxDisk(Lft,Rgt,Bas,Top)
!        CALL Fill_Circle(Lft,Rgt,Bas,Top)
!    END SUB
!END MODULE

!! ------ End of TB Unix Header and Subs ------
! -----------------------------------------------------------

! ***

! One could display trajectories through the following points.
! I'll take  theta = 0;  for other theta, apply  R(theta)  to what I say.
!
! In real eigenvalue case, show the two eigenlines which go through
!    a1 = (1/2,0)  and  a2 = (s/2, omega)
! and the trajectories though the four vectors
!    a1+a2,  a1-a2,  -a1+a2,  -a1-a2.
!
! For the complex eigenvalue case use 
!    (1/2,0), (-1/2,0)
! and
!    +-((s-1)/2,0)  if  s>0
!    +-((s+1)/2,0)  if  s<0.
!
! One could compute one of each opposite pair using RK.

! ***

SUB LinearParameters
    DECLARE PUBLIC PCFlag,Mac5Flag,M68KFlag,UnixFlag,xmax,ymax
    DECLARE PUBLIC black,drkgry,drkmid,midgry,litmid,litgry,white
    DECLARE PUBLIC red,yellow,green,cyan,blue,magenta,colorscheme
    DECLARE PUBLIC axislabelclr,slideclr,true,false
    DECLARE PUBLIC workLft,workRgt,workBas,workTop,workMidx

    DECLARE DEF InfoWithin,QuitWithin

    ! --- help screen array ---

    DIM  info$(1:11)
    MAT READ info$
    DATA "Linear Phase Portraits: Cursor Entry"
    DATA ""
    DATA "Linear Phase Portraits supports exploration of the phase portraits of linear homogeneous systems with constant coefficients."
    DATA ""
    DATA "You can choose the linear system using the sliders, or by selecting a point on the trace,determinant plane or on the theta,s plane."
    DATA ""
    DATA "Click a slider among the tick marks to set a value rounded to the nearest .25."
    DATA "Click or drag the mouse in the trace,determinant plane to create a matrix with the trace and determinant at the cursor, and display the corresponding vector field and eigenlines."
    DATA "Click or drag the mouse in the theta,s plane to create matrices with the active trace and determinant."
    DATA ""
    DATA "Click a point in the x,y plane to set initial conditions and see the corresponding solution curve."

    !   trace = a+d
    !   dterm = ad - bc
    !   dscrm = trace^2 - 4 dterm
    !   dterm = trace^2/4 is the parabola on the trace,dterm plane

    !   characteristic polynomial = lam^2 - trace*lam + dterm
    !     lam = (trace +- sqr(trace^2 - 4 dterm)) / 2

    !   IF dterm<trace^2/4, dscrm>0
    !      outside the parabola - two real lams
    !      the theta,s plane is continous

    !   IF dterm=trace^2/4, dscrm=0
    !      on the parabola - one real lam = trace/2
    !      the theta,s plane is continous - the horizontal axis matters

    !   IF dterm>trace^2/4, dscrm<0
    !      inside the parabola - two complex lams
    !      the theta,s plane has two active regions
    !      s^2 must be >= (t2^2 - dterm)


    ! ----- modify the palette -----

    LET b0= 0                                    ! web safe intervals
    LET b1= 0.20
    LET b2= 0.40
    LET b3= 0.60
    LET b4= 0.80
    LET b5= 1.00

    SET COLOR MIX( 7)  b5,b3,b1                  ! red
    SET COLOR MIX( 8)  b5,b5,b0                  ! yellow
    SET COLOR MIX( 9)  b5,b2,b1                  ! green
    SET COLOR MIX(10)  b0,b4,b3                  ! cyan
    SET COLOR MIX(11)  b2,b3,b5                  ! blue
    SET COLOR MIX(12)  b0,b4,b1                  ! magenta
    SET COLOR MIX(13)  b5,b1,b5                  ! pink
    FOR clr= 7 to 13
        SET COLOR clr
    NEXT clr

    LET red    =  7
    LET yellow =  8
    LET green  =  9
    LET cyan   = 10
    LET blue   = 11
    LET magenta= 12
    LET pink   = 13

    LET axisclr= drkmid

    ! ---------- Utility functions ---

    DECLARE DEF clamp, roundn, e
    DEF notanum= 987656789


    ! --- functions ---

    DEF dxdt(x,y)= ma*x + mb*y
    DEF dydt(x,y)= mc*x + md*y

    DEF fT(ma,mb,mc,md)  = ma+md
    DEF fD(ma,mb,mc,md)  = ma*md - mb*mc
    DEF fDsc(trace,dterm)= round(trace*trace - 4*dterm,8)

    SUB MatToTD(ma,mb,mc,md,trace,dterm)
        LET trace= fT(ma,mb,mc,md)               ! find trace
        LET dterm= fD(ma,mb,mc,md)               ! find determinant
    END SUB

    DEF omIV(dterm,trace)
        LET omIV= sqr(abs(dterm - trace^2/4))
    END DEF

    DEF om(dterm,trace)
        IF dscrm<0 then
           LET om= sqr(abs(dterm - trace^2/4))
        ELSE
           LET om= 0
        END IF
    END DEF


    ! ---------- Graphing plane parameters and methods ----------

    ! --- w1 plane trace determinate ---

    DECLARE PUBLIC w1Lft,w1Rgt,w1Bas,w1Top,w1Midx,w1Midy
    DECLARE PUBLIC w1fLft,w1fRgt,w1fBas,w1fTop,w1x0,w1y0
    DECLARE PUBLIC w1xFirst, w1xSTik, w1xLTik, w1xLabel, w1xGridstep
    DECLARE PUBLIC w1yFirst, w1ySTik, w1yLTik, w1yLabel, w1yGridstep
    DECLARE PUBLIC w1wWid,w1wHgt,w1fWid,w1fHgt
    DECLARE PUBLIC w1fxRatio,w1fyRatio,w1wxRatio,w1wyRatio,w1Aspect
    DECLARE PUBLIC w1xPiFlag, w1xMult, w1yPiFlag, w1yMult

    LET w1Flag   = 1
    LET w1xPiFlag= 0
    LET w1xMult  = 1
    LET w1yPiFlag= 0
    LET w1yMult  = 1

    LET w1Lft = workLft +  80                    ! pixel bounds
    LET w1Rgt = w1Lft  +  160                    ! for 6
    LET w1Top = workTop + 300
    LET w1Bas = w1Top  +  160

    LET w1size=  4
    LET w1fLft= -w1size                          ! function bounds
    LET w1fRgt=  w1size
    LET w1fTop=  w1size
    LET w1fBas= -w1size

    LET w1Xax$= "Tr A"                           ! axis labels
    LET w1Yax$= "Det A"

    LET w1xGridstep= 0                           ! horizontal grid intervals
    LET w1yGridstep= 0                           ! vertical grid intervals

    LET w1xSTik = 0                              ! horizontal axis Tik marks
    LET w1xLTik = 1
    LET w1xLabel= 2
    LET w1xFirst= w1fLft

    LET w1ySTik = 0                              ! vertical axis Tik marks
    LET w1yLTik = 1
    LET w1yLabel= 2
    LET w1yFirst= w1fBas

    ! --- Plane 1 methods ---

    DECLARE DEF w1Fncx,w1Fncy,w1Wndx,w1Wndy      ! window/function transforms
    DECLARE DEF w1wWithin, w1Within, w1fWithin

    CALL w1Variables

    SUB w1Init
        CALL w1DrawPlane(0,0,1)                  ! grid, axes, zeroaxes
        CALL DrawTDPlane
        CALL PlotTextCJ(w1x0,w1Top-10,w1Yax$,white)
        CALL PlotTextLJ(w1Rgt+8,w1y0+3,w1Xax$,white)
        CALL w1KeepGridLayer
        CALL w1KeepGraphLayer
    END SUB


    ! ------------------------------------------
    ! --- w3 plane - phase plane ---

    DECLARE PUBLIC w3Lft,w3Rgt,w3Bas,w3Top,w3Midx,w3Midy
    DECLARE PUBLIC w3fLft,w3fRgt,w3fBas,w3fTop,w3x0,w3y0
    DECLARE PUBLIC w3xFirst, w3xSTik, w3xLTik, w3xLabel, w3xGridstep
    DECLARE PUBLIC w3yFirst, w3ySTik, w3yLTik, w3yLabel, w3yGridstep
    DECLARE PUBLIC w3wWid,w3wHgt,w3fWid,w3fHgt
    DECLARE PUBLIC w3fxRatio,w3fyRatio,w3wxRatio,w3wyRatio,w3Aspect
    DECLARE PUBLIC w3xPiFlag, w3xMult, w3yPiFlag, w3yMult

    LET w3Flag   = 1
    LET w3xPiFlag= 0
    LET w3xMult  = 1
    LET w3yPiFlag= 0
    LET w3yMult  = 1

    LET w3Lft = workrgt- 400                     ! pixel bounds
    LET w3Rgt = w3Lft +  360
    LET w3Top = worktop + 45
    LET w3Bas = w3Top +  360

    LET w3size=  4
    LET w3fLft= -w3size                          ! function bounds * pi
    LET w3fRgt=  w3size
    LET w3fTop=  w3size
    LET w3fBas= -w3size

    LET w3xGridstep= 0                           ! horizontal grid intervals
    LET w3yGridstep= 0                           ! vertical grid intervals

    LET w3xAx$= "x"                              ! axis labels
    LET w3yAx$= "y"

    LET w3xSTik = 1/2                            ! horizontal axis Tik marks
    LET w3xLTik = 1
    LET w3xLabel= 1
    LET w3xFirst= w3fLft

    LET w3ySTik = 1/2                            ! vertical axis Tik marks
    LET w3yLTik = 1
    LET w3yLabel= 1
    LET w3yFirst= w3fBas

    ! --- plane 2 methods ---

    DECLARE DEF w3Fncx, w3Fncy, w3Wndx, w3Wndy   ! window/function transforms
    DECLARE DEF w3wWithin, w3Within, w3fWithin

    CALL w3Variables

    SUB w3Init
        CALL w3DrawPlane(0,0,1)                  ! grid, axes, zeroaxes
        CALL SetTextFont(1,12,"normal")
        CALL PlotTextLJ(w3Rgt+8,w3y0+3,w3xAx$,white)  ! axis labels
        CALL PlotTextCJ(w3x0,w3Top-10,w3yAx$,white)
        CALL SetMatDiffEq

        CALL w3KeepGridLayer
        CALL w3KeepGraphLayer
        CALL w3KeepFieldLayer
    END SUB

    SUB w3Refresh
        CALL w3ShowGridLayer
        CALL w3KeepGraphLayer
        CALL w3KeepFieldLayer
    END SUB

    SUB w3ShowFieldLayer
        BOX SHOW w3FieldLayer$ at w3Lft-5,w3Bas+5
        CALL w3KeepGraphLayer
    END SUB

    SUB w3KeepFieldLayer
        BOX KEEP w3Lft-5,w3Rgt+5,w3Bas+5,w3Top-5 in w3FieldLayer$
    END SUB

    SUB SetMatDiffEq
        LET DEBas= w3Top - 28
        CALL SetTextFont(1,12,"bold")
        LET txt$ = "x = Ax"
        CALL StringWidth(txt$,sw)
        LET DELft= w3x0 - sw/2
        CALL PlotTextLJ(DELft,DEBas,txt$,white)
        CALL PlotTextLJ(DELft+2,DEBas-9,".",white)
    END SUB

    ! -----------------------------------
    ! --- w2 plane - Theta,S plane ---

    DECLARE PUBLIC w2Lft,w2Rgt,w2Bas,w2Top,w2Midx,w2Midy
    DECLARE PUBLIC w2fLft,w2fRgt,w2fBas,w2fTop,w2x0,w2y0
    DECLARE PUBLIC w2xFirst, w2xSTik, w2xLTik, w2xLabel, w2xGridstep
    DECLARE PUBLIC w2yFirst, w2ySTik, w2yLTik, w2yLabel, w2yGridstep
    DECLARE PUBLIC w2wWid,w2wHgt,w2fWid,w2fHgt
    DECLARE PUBLIC w2fxRatio,w2fyRatio,w2wxRatio,w2wyRatio,w2Aspect
    DECLARE PUBLIC w2xPiFlag, w2xMult, w2yPiFlag, w2yMult

    LET w2Flag   = 1
    LET w2xPiFlag= 1
    LET w2xMult  = pi
    LET w2yPiFlag= 0
    LET w2yMult  = 1

    LET w2Lft = w1Lft                            ! pixel bounds
    LET w2Rgt = w1Rgt
    LET w2Top = w3Top
    LET w2Bas = w2Top + 160

    LET w2fLft=  0                               ! function bounds
    LET w2fRgt=  1
    LET w2fTop=  4
    LET w2fBas= -4

    LET w2xAx$= "theta"                          ! axis labels
    LET w2yAx$= "s"

    LET w2xGridstep= 0                           ! grid line intervals
    LET w2yGridstep= 0

    LET w2xSTik = 0                              ! axis Tik marks
    LET w2xLTik = 1
    LET w2xLabel= 1
    LET w2xFirst= 0

    LET w2ySTik = 1
    LET w2yLTik = 2
    LET w2yLabel= 2
    LET w2yFirst= w2fBas

    ! --- plane 3 methods ---

    DECLARE DEF w2Fncx,w2Fncy,w2Wndx,w2Wndy      ! window/function transforms
    DECLARE DEF w2wWithin,w2Within,w2fWithin

    CALL w2Variables

    SUB w2Init
        CALL w2DrawPlane(0,0,0)                  ! grid, axes, zeroaxes
        SET COLOR litmid
        PLOT w2Lft,w2y0; w2Rgt,w2y0
        CALL w2KeepGridLayer
        CALL w2KeepGraphLayer
    END SUB

    ! --- plane 4 data: bifurcation plane ---


    ! ------ matrix ------

    LET lnspc= 15

    LET mtop = w3Bas + 73
    LET mbas = mtop +  38
    LET mlft = w3Lft + 30
    LET mrgt = mlft + 115
    LET matclr= white

    LET mtxt0= mtop+14
    LET mtxt1= mtop+33
    LET mcol1= mlft+10
    LET mcol2= mcol1+50

    LET mmidx= int((mlft+mrgt)/2)

    ! ---- matrices ----

    DIM Rm(1:2,1:2)
    DIM Rp(1:2,1:2)
    DIM Ht(1:2,1:2)
    DIM As(1:2,1:2)
    DIM A1(1:2,1:2)
    DIM  A(1:2,1:2)

    SUB ClearMatrixValues
        BOX CLEAR mlft+3,mrgt-3,mtxt0+3,mtxt0-10      ! upper row
        BOX CLEAR mlft+3,mrgt-3,mtxt1+3,mtxt1-10      ! lower row
    END SUB

    SUB DrawMatrix
        CALL SetTextFont(1,12,"bold")
        LET midm= int((mtop+mbas)/2)
        CALL PlotTextRJ(mlft-5,midm+3,"A =",white)

        SET COLOR litmid                         ! matrix brackets
        CALL SetLineWeight(2)
        PLOT mlft+4,mtop; mlft,mtop; mlft,mbas; mlft+4,mbas
        PLOT mrgt-4,mtop; mrgt,mtop; mrgt,mbas; mrgt-4,mbas
        CALL SetLineWeight(1)
    END SUB

    SUB ShowMatrixValues
        CALL SetTextFont(1,12,"bold")
        CALL ClearMatrixValues

        LET places= 5
        LET rma= round(ma,places)
        LET rmb= round(mb,places)
        LET rmc= round(mc,places)
        LET rmd= round(md,places)

        CALL PlotTextLJ(mcol1,mtxt0,using$("--%.##",rma),matclr)
        CALL PlotTextLJ(mcol2,mtxt0,using$("--%.##",rmb),matclr)
        CALL PlotTextLJ(mcol1,mtxt1,using$("--%.##",rmc),matclr)
        CALL PlotTextLJ(mcol2,mtxt1,using$("--%.##",rmd),matclr)
        CALL LamValues
    END SUB

    ! ----------- vertical sliders ------------
    ! ----------- vertical slider 1 determinant ---

    DECLARE PUBLIC v1axis,v1wLft,v1wRgt,v1wBas,v1wTop,v1sBas,v1sTop
    DECLARE PUBLIC v1fBas,v1fTop,v1First,v1STik,v1LTik,v1Label
    DECLARE PUBLIC v1name$,v1form$,v1clr,v1PiAxis,v1Mult

    LET v1PiAxis= 0
    LET v1Mult  = 1

    LET v1Clr   = white
    LET v1name$ = ""
    LET v1form$ = "-%.##"
    LET v1Places= 2
    LET v1Click = 1/4

    LET v1axis  = w1Lft - 30
    LET v1wBas  = w1Bas
    LET v1wTop  = w1Top
    LET v1fBas  = -4
    LET v1fTop  =  4

    LET v1STik  =  1/4
    LET v1LTik  =  1/2
    LET v1Label =  1
    LET v1First = v1fBas

    DECLARE DEF v1Fncy,v1Wndy                    ! window/function transforms
    DECLARE DEF v1Within

    CALL v1SliderVariables

    SUB v1Init
        CALL v1DrawSlider(v1name$,dterm)
    END SUB

    ! --- vertical slider 2 ---

    DECLARE PUBLIC v2axis,v2wLft,v2wRgt,v2wBas,v2wTop,v2sBas,v2sTop
    DECLARE PUBLIC v2fBas,v2fTop,v2First,v2STik,v2LTik,v2Label
    DECLARE PUBLIC v2name$,v2form$,v2clr,v2PiAxis,v2Mult

    LET v2PiAxis= 0
    LET v2Mult  = 1

    LET v2Clr   = white
    LET v2name$ = "s"
    LET v2form$ = "-%.##"
    LET v2Places= 2
    LET v2Click = 1/4

    LET v2axis  = w2Lft - 30
    LET v2wBas  = w2Bas
    LET v2wTop  = w2Top
    LET v2fTop  =  4
    LET v2fBas  = -4

    LET v2STik  =  1/4
    LET v2LTik  =  1/2
    LET v2Label =  1
    LET v2First = v2fBas

    DECLARE DEF v2Fncy,v2Wndy                    ! window/function transforms
    DECLARE DEF v2Within

    CALL v2SliderVariables

    SUB v2Init
        CALL v2DrawSlider(v2name$,s)
    END SUB

    ! ---- s and theta sliders ----
    ! ---- slider 5 ---

    DECLARE PUBLIC h1axis,h1wLft,h1wRgt,h1wBas,h1wTop,h1fLft,h1fRgt
    DECLARE PUBLIC h1name$,h1form$,h1clr,h1First,h1STik,h1LTik,h1Label
    DECLARE PUBLIC h1PiAxis,h1Mult,h1fMin,h1fMax

    LET h1PiAxis= 0
    LET h1Mult  = 1

    LET h1clr   = white
    LET h1name$ = ""
    LET h1form$ = "--%.##"
    LET h1Places= 2
    LET h1Click = 1/4

    LET h1axis  = w1Bas + 30
    LET h1wLft  = w1Lft
    LET h1wRgt  = w1Rgt
    LET h1fLft  = -4
    LET h1fRgt  =  4

    LET h1STik  = 1/4
    LET h1LTik  = 1/2
    LET h1Label = 1
    LET h1First = h1fLft

    DECLARE DEF h1Within,h1Fncx

    CALL h1SliderVariables

    SUB h1Init
        CALL h1DrawSlider(h1name$,trace)
    END SUB

    ! --- h2 theta slider ---

    DECLARE PUBLIC h2axis,h2wLft,h2wRgt,h2wBas,h2wTop,h2fLft,h2fRgt
    DECLARE PUBLIC h2name$,h2form$,h2clr,h2First,h2STik,h2LTik,h2Label
    DECLARE PUBLIC h2PiAxis,h2Mult,h2fMin,h2fMax

    LET h2PiAxis= 1
    LET h2Mult  = pi

    LET h2clr   = white
    LET h2name$ = ""
    LET h2form$ = "--%.##"
    LET h2Places= 2
    LET h2Click = pi/4

    LET h2axis  = w2Bas + 30
    LET h2wLft  = w2Lft
    LET h2wRgt  = w2Rgt

    LET h2fLft  = 0
    LET h2fRgt  = 1

    LET h2STik  = 0
    LET h2LTik  = 1/4
    LET h2Label = 1/4
    LET h2First = 0

    DECLARE DEF h2Within,h2Fncx

    CALL h2SliderVariables

    SUB h2Init
        CALL h2DrawSlider(h2name$,theta)
        CALL DrawTheta12(h2wlft-16,h2wbas-1,white)    ! slider variable label
    END SUB

    ! ---- end of slider parameters ----


    ! -------- text locations --------
    ! ---- initial values x and y ----

    LET ivlft = w3Rgt - 100
    LET ivrgt = w3Rgt
    LET ivbas1= w3Bas +  24
    LET ivbas2= ivbas1+  20
    LET ivbas = ivbas2 +  5
    LET ivtop = ivbas1 - 15
    LET iveqx = ivlft +  50

    SUB InitValLabels
        BOX CLEAR ivlft,ivrgt,ivbas,ivtop
        CALL SetTextFont(1,12,"normal")
        LET x$= "x(0) = "
        CALL PlotTextRJ(iveqx,ivbas1,x$,classClr)
        LET y$= "y(0) = "
        CALL PlotTextRJ(iveqx,ivbas2,y$,classClr)
    END SUB

    SUB InitValValues(x0,y0)
        CALL SetTextFont(1,12,"bold")
        CALL InitValClear
        LET x$= using$("-%.##",x0)
        CALL PlotTextLJ(iveqx,ivbas1,x$,classClr)
        LET y$= using$("-%.##",y0)
        CALL PlotTextLJ(iveqx,ivbas2,y$,classClr)
    END SUB

    SUB InitValClear
        BOX CLEAR iveqx-2,ivrgt,ivbas,ivtop
    END SUB

    ! ---- classification name ----

    LET lclft= w3lft
    LET lcbas= ivbas2

    SUB ClassText(className$)
        CALL ClassTextClear
        CALL SuperSubScriptLJ(lcLft,lcBas,className$,classClr)
        CALL InitValLabels
    END SUB

    SUB ClassTextClear
        BOX CLEAR lcLft,lcLft+170,lcbas+5,lcbas-15    ! className$
    END SUB


    ! ---- lam text ---

    LET lamlft = mrgt  + 40
    LET lamrgt = lamlft+120
    LET lambas1= mtxt0
    LET lambas2= mtxt1
    LET lameqx = lamlft + 48

    SUB LamLabels
        CALL LamClear

        CALL DrawLam(lameqx-30,lambas1,white)
        CALL DrawLam(lameqx-30,lambas2,white)

        CALL SetTextFont(1,9,"bold")
        CALL PlotTextLJ(lameqx-20,lambas1+4,"1",white)
        CALL PlotTextLJ(lameqx-20,lambas2+4,"2",white)

        CALL SetTextFont( 1,12,"bold")
        CALL PlotTextRJ(lameqx,lambas1,"=",white)
        CALL PlotTextRJ(lameqx,lambas2,"=",white)

        CALL LamValues
    END SUB

    SUB LamValues
        CALL FindLamdas
        CALL LamClear

        CALL SetTextFont(1,12,"bold")
        CALL PlotTextLJ(lameqx+4,lambas1,l1$,white)
        CALL PlotTextLJ(lameqx+4,lambas2,l2$,white)
    END SUB

    SUB LamClear
        BOX CLEAR lameqx-2,lamlft+155,lambas2+5,lambas1-15
    END SUB

    ! ----- buttons -----

    LET bhgt= 18
    LET ctop= w3Bas+10
    LET cbas= ctop + bhgt
    LET clft= w3lft
    LET crgt= clft+47

    SUB Buttons
        CALL SetTextFont(1,12,"bold")
        CALL DrawButton(clft,crgt,cbas,ctop,5,"clear")
    END SUB

    ! ---- end of text locations ----


    ! ----- set default vars and initialize -----

    IF M68kFlag=1 then                           ! mac TB 4.04
       LET tstep= 1/128
    ELSE
       LET tstep= 1/64
    END IF

    LET trace,oldtrace=  0
    LET dterm,olddterm= -1
    LET s,olds,ss     =  0
    LET theta,oldtheta=  0
    LET omega=  0
    LET close=  0.05                             ! close to parabola?

    ! IF UnixFlag=1 then CALL SetTimer

    CALL InitScreen

    SUB InitScreen
        BOX CLEAR worklft,workrgt,workbas,worktop

        CALL w1Init
        CALL w2Init
        CALL w3Init

        CALL TraceDetUpdate

        CALL v1Init
        CALL h1Init

        CALL h2Init
        CALL v2Init

        CALL MarkTSCursor(theta,s)

        CALL Buttons
        CALL DrawMatrix

        CALL InitValLabels
        CALL LamLabels
    END SUB


    ! --- event manager ---

    DO
       LET clearFlag2=  0
       LET omx,omy   = -999
       DO
          GET MOUSE: mx,my,ms
          IF w3Within(mx,my)=true then
             CALL w3wClamp(mx,my)
             IF omx<>mx or omy<>my then
                CALL w3RollOver(mx,my)
                CALL SetVars(mx,my,omx,omy)
             END IF
          ELSE IF clearFlag2=true then
             CALL w3RollOverClear
          END IF
       LOOP until ms=2                           ! mouse down?

       IF w3Within(mx,my)=true then              ! x,y Phase Plane
          LET omx,omy= -999
          DO
             GET MOUSE: mx,my,ms
             IF w3Within(mx,my)=true then
                CALL w3wClamp(mx,my)
                IF omx<>mx or omy<>my then
                   CALL w3PixelsToMath(mx,my,x0,y0)
                   CALL w3ShowGraphLayer
                   CALL SetVector(x0,y0,x2,y2,.2)
                   CALL InitValValues(x0,y0)
                   CALL SetVars(mx,my,omx,omy)
                END IF
             END IF
          LOOP until ms=3
          CALL w3ShowGraphLayer
          CALL Trajectory(x0,y0)
          CALL w3KeepGraphLayer

          ! --- trace/determinant window ---

       ELSE IF w1Within(mx,my)=true then         ! T,D Parameter Plane
          LET oldw1x,oldw1y= -999
          DO
             GET MOUSE: mx,my,ms
             IF ms=3 then EXIT DO
             IF w1Within(mx,my)=true then
                CALL w1wClamp(mx,my)
                IF oldw1x<>mx or oldw1y<>my then
                   CALL w1PixelsToMath(mx,my,trace,dterm)
                   CALL TDProximity
                   CALL TraceDetUpdate
                   CALL h1Mark(trace)
                   CALL v1Mark(dterm)
                   CALL SetVars(trace,dterm,oldtrace,olddterm)
                   CALL SetVars(mx,my,oldw1x,oldw1y)
                END IF
             END IF
          LOOP until ms=3

          ! --- Theta,s window ---

       ELSE IF w2Within(mx,my)=true then         ! Theta,s Parameter Plane
          LET oldw3x,oldw3y= -999999
          DO
             GET MOUSE: mx,my,ms
             IF ms=3 then EXIT DO
             IF w2Within(mx,my)=true then
                CALL w2wClamp(mx,my)
                IF mx<>oldw3x or my<>oldw3y then
                   CALL w2PixelsToMath(mx,my,theta,s)
                   CALL v2ProximityCheck
                   CALL SThetaUpdate
                   CALL v2Mark(s)
                   CALL h2Mark(theta)
                   CALL SetVars(s,theta,olds,oldtheta)
                   CALL SetVars(mx,my,oldw3x,oldw3y)
                   IF s<>ss then LET ss= s
                END IF
             END IF
          LOOP until ms=3

          ! --- trace/determinant sliders ---

       ELSE IF v1Within(mx,my)=true then         ! determinant
          IF mx<v1axis then
             CALL v1MouseClick
          ELSE
             CALL v1MouseDrag
          END IF
       ELSE IF h1Within(mx,my)=true then         ! trace
          IF my<h1axis then
             CALL h1MouseClick
          ELSE
             CALL h1MouseDrag
          END IF


          ! --- Theta,s sliders ---

       ELSE IF v2Within(mx,my)=true then         ! s param
          IF mx<v2axis then
             CALL v2MouseClick
          ELSE
             CALL v2MouseDrag
          END IF

       ELSE IF h2Within(mx,my)=true then         ! theta param
          IF my<h2axis then
             CALL h2MouseClick
          ELSE
             CALL h2MouseDrag
          END IF

       ELSE IF mx>clft and mx<crgt and my>ctop and my<cbas then      ! clear
          CALL MouseButtonUp(cLft,cRgt,cBas,cTop,ms)
          CALL w3ShowFieldLayer

       ELSE IF infoWithin(mx,my,ms)=true then
          CALL infoButtonUp(ms)
          CALL InfoPage(Info$)
          CALL InitScreen
       ELSE IF quitWithin(mx,my,ms)=true then
          CALL quitButtonUp(ms)
          EXIT SUB
       ELSE
          CALL MouseUp(mx,my,ms)
       END IF
    LOOP

    ! ----- w3 x,y window events -----

    SUB w3RollOver(mx,my)
        CALL w3PixelsToMath(mx,my,x0,y0)
        CALL w3ShowGraphLayer
        CALL SetVector(x0,y0,x2,y2,.2)
        CALL InitValValues(x0,y0)
        LET clearFlag2= 1
    END SUB

    SUB w3RollOverClear
        CALL w3ShowGraphLayer
        CALL InitValClear
        LET clearFlag2= 0
    END SUB

    ! ------ update matrix ------

    SUB SetMatrixValues
        LET L  = sqr(s^2 + td)
        LET Lcs= L*cos(2*theta)
        LET Lsn= L*sin(2*theta)

        LET ma = t2 - Lsn                        ! define matrix elements
        LET mb = Lcs + s
        LET mc = Lcs - s
        LET md = t2 + Lsn
    END SUB

    ! ----- trace and determinant sliders -----

    SUB v1MouseClick
        CALL MouseUp(mx,my,ms)
        LET my   = clamp(my,v1wTop,v1wBas)
        LET dterm= roundn(v1Fncy(my),v1Click)
        CALL v1EventAction
    END SUB
    SUB v1MouseDrag
        DO
           GET MOUSE: mx,my,ms
           LET my   = clamp(my,v1wTop,v1wBas)
           LET dterm= round(v1Fncy(my),v1Places)
           CALL v1EventAction
        LOOP until ms=3
    END SUB
    SUB v1EventAction
        CALL DtermProximity
        IF dterm<>olddterm then
           CALL v1mark(dterm)
           CALL TraceDetUpdate
           LET  olddterm= dterm
        END IF
    END SUB

    SUB DtermProximity
        LET parab= trace^2/4                     ! near parabola?
        IF abs(dterm-parab)<close then
           LET dterm= parab
        END IF
        IF abs(dterm)<close then LET dterm= 0    ! near zero?
    END SUB

    ! -----

    SUB h1MouseClick
        CALL MouseUp(mx,my,ms)
        LET mx   = clamp(mx,h1wLft,h1wRgt)
        LET trace= roundn(h1Fncx(mx),h1Click)
        CALL h1EventAction
    END SUB
    SUB h1MouseDrag
        DO
           GET MOUSE: mx,my,ms
           LET mx   = clamp(mx,h1wLft,h1wRgt)
           LET trace= round(h1Fncx(mx),h1Places)
           CALL h1EventAction
        LOOP until ms=3
    END SUB
    SUB h1EventAction
        CALL TraceProximity
        IF trace<>oldtrace then
           CALL h1mark(trace)
           CALL TraceDetUpdate
           LET oldtrace= trace
        END IF
    END SUB

    SUB TraceProximity
        IF dterm>=0 then                         ! upper half plane?
           ! check for parabola
           LET signtr= sgn(trace)
           LET x     = sqr(4*dterm)
           IF abs(abs(trace)-x)<close then
              LET trace= signtr*x
              LET dterm= trace^2/4
           END IF
        END IF
        IF abs(trace)<close then LET trace= 0
    END SUB

    ! -----

    SUB TDProximity
        CALL TraceProximity
        CALL DtermProximity
    END SUB

    SUB TraceDetUpdate
        CALL TraceDetToMatrix(trace,dterm,theta,s)
        CALL MarkTDCursor(trace,dterm)
        CALL ShowMatrixValues
        CALL Classification(trace,dterm,dscrm,className$,classClr)
        CALL DrawThetaPlane(trace,dterm)

        CALL DirectionField
        CALL FindEigenVectors(ma,mb,mc,md)
        CALL w3KeepGraphLayer
    END SUB

    SUB TraceDetToMatrix(trace,dterm,theta,s)
        ! Create the A matrix from trace,dterm and theta,s

        LET dscrm= fDsc(trace,dterm)             ! trace^2 - 4*dterm
        LET t2   = trace/2
        LET td   = trace^2/4 - dterm

        IF dscrm<0 then                          ! in parabola
           LET omega= sqr(abs(td))+.00001        ! minimum s value
           IF abs(s)<omega then
              IF s>=0 then
                 LET s=  omega
              ELSE
                 LET s= -omega
              END IF
              IF s<>olds then
                 CALL v2Mark(s)
                 LET olds= s
              END IF
           END IF
        ELSE
           LET omega= 0
        END IF

        CALL SetMatrixValues

        IF round(dscrm,8)=0 then CALL Classification(trace,dterm,dscrm,className$,classClr)
    END SUB

    SUB TestResults                              ! debugging routine - don't translate
        LET tr= ma+md
        LET dt= ma*md - mb*mc
        LET ds= trace^2 - 4*dterm
        BOX CLEAR 0,300,60,0
        SET CURSOR 1,1
        PRINT "trace = " & using$("--%.#####",tr)
        PRINT "dterm = " & using$("--%.#####",dt)
        PRINT "dscrm = " & using$("--%.#####",ds)
    END SUB

    ! --- end of trace/determinate action ---


    ! ----- s and theta plane -----

    SUB v2MouseClick
        CALL MouseUp(mx,my,ms)                   ! needed here for gap
        LET my= clamp(my,v2wTop,v2wBas)
        LET s = roundn(v2Fncy(my),v2Click)
        CALL v2EventAction
    END SUB
    SUB v2MouseDrag
        DO
           GET MOUSE: mx,my,ms
           LET my= clamp(my,v2wTop,v2wBas)
           LET s = round(v2Fncy(my),v2Places)
           CALL v2EventAction
        LOOP until ms=3
    END SUB
    SUB v2EventAction
        CALL v2ProximityCheck
        IF s<>olds then
           CALL v2Mark(s)
           CALL SThetaUpdate
           IF round(dscrm,8)= 0 then             ! star or defective node?
              CALL Classification(trace,dterm,dscrm,className$,classClr)
           END IF
           LET olds= s
           IF s<>ss then LET ss= s
        END IF
    END SUB

    SUB v2ProximityCheck
        IF omega>0 then
           IF abs(s)<omega then
              IF s>=0 then
                 LET s=  omega+.00001
              ELSE IF s<0 then
                 LET s= -omega-.00001
              END IF
           END IF
        ELSE
           IF abs(s)<close then LET s= 0
        END IF
    END SUB

    ! ---

    SUB h2MouseClick
        CALL MouseUp(mx,my,ms)
        LET mx   = clamp(mx,h2wLft,h2wRgt)
        LET theta= roundn(h2Fncx(mx),h2Click)
        CALL h2EventAction
    END SUB
    SUB h2MouseDrag
        DO
           GET MOUSE: mx,my,ms
           LET mx   = clamp(mx,h2wLft,h2wRgt)
           LET theta= round(h2Fncx(mx),h2Places)
           CALL h2EventAction
        LOOP until ms=3
    END SUB
    SUB h2EventAction
        IF theta<>oldtheta then
           CALL h2Mark(theta)
           CALL SThetaUpdate
           LET oldtheta= theta
        END IF
    END SUB

    ! ----

    SUB SThetaUpdate
        CALL ThetaSToMatrix(trace,dterm,theta,s)
        CALL MarkTSCursor(theta,s)
        CALL ShowMatrixValues
        CALL DirectionField
        CALL FindEigenVectors(ma,mb,mc,md)
        CALL w3KeepGraphLayer
    END SUB

    SUB ThetaSToMatrix(trace,dterm,theta,s)
        ! Create the A matrix from trace,dterm and theta,s

        LET t2 = trace/2
        LET td = trace^2/4 - dterm

        CALL SetMatrixValues
    END SUB

    ! ------- end of slider routines -------


    ! -------- Draw w2 Theta s Plane ---------

    SUB DrawThetaPlane(trace,dterm)
        LOCAL wop,won
        IF dscrm=0 then                          ! on the parabola
           ! s= 0 is a star
           ! s<>0 is a defective node

           CALL w2ShowGridLayer
           CALL PlotLine(w2Lft,w2y0,w2Rgt,w2y0,yellow)
        ELSE IF dscrm<0 then                     ! inside the parabola
           !  The theta,s plane has two active regions - the region
           !  from -omega to omega on the vertical is inactive

           CALL SetVars(w2Wndy(omega),w2Wndy(-omega),wop,won)
           CALL BoxClear(w2Lft-3,w2Rgt+3,w2Bas+3,w2Top-3)

           CALL BoxArea(w2Lft+1,w2Rgt-1,wop-1,w2Top+1,black)
           CALL BoxLines(w2Lft,w2Rgt,wop,w2Top,drkmid)
           CALL PlotLine(w2Lft,wop,w2Rgt,wop,yellow)

           CALL BoxArea(w2Lft+1,w2Rgt-1,w2Bas-1,won+1,black)
           CALL BoxLines(w2Lft,w2Rgt,w2Bas,won,drkmid)
           CALL PlotLine(w2Lft,won,w2Rgt,won,yellow)
        ELSE IF dscrm>0 then                     ! outside the parabola
           !  The theta,s plane is continous

           CALL w2ShowGridLayer
        END IF

        CALL w2KeepGraphLayer
        CALL MarkTSCursor(theta,s)
    END SUB


    ! ---- Draw w1 Trace Determinant plane ----

    SUB DrawTDPlane
        CALL BoxArea(w1Lft,w1Rgt,w1Bas-1,w1y0,blue)
        CALL BoxArea(w1Lft,w1x0-1,w1y0-1,w1Top,green)
        CALL BoxArea(w1x0+1,w1Rgt,w1y0,w1Top,magenta)
        CALL PlotLine(w1Lft,w1y0  ,w1Rgt,w1y0  ,pink)
        CALL PlotLine(w1Lft,w1y0+1,w1Rgt,w1y0+1,pink)

        FOR wx= w1Lft to w1Rgt                   ! fill inside parabola
            LET fx= w1Fncx(wx)
            LET fy= 0.25*fx*fx
            LET wy= w1Wndy(fy)
            IF fx<0 then
               CALL PlotLine(wx,wy, wx,w1Top, red)
            ELSE
               CALL PlotLine(wx,wy, wx,w1Top, cyan)
            END IF
        NEXT wx

        CALL SetLineWeight(2)
        SET COLOR yellow                         ! outline parabola
        IF m68kFlag=1 or unixFlag=1 then
           FOR wx= w1Lft to w1Rgt
               LET fx= w1Fncx(wx)
               LET fy= 0.25*fx*fx
               LET wy= w1Wndy(fy)
               IF round(wy)=w1y0 then LET wy= wy-1
               IF fx<0 then
                  PLOT wx,wy;
               ELSE                              ! shift for symmetry
                  PLOT wx-1,wy;
               END IF
           NEXT wx
        ELSE
           FOR wx= w1Lft to w1Rgt
               LET fx= w1Fncx(wx)
               LET fy= 0.25*fx*fx
               LET wy= w1Wndy(fy)+1
               PLOT wx+1,wy;
           NEXT wx
        END IF

        PLOT
        CALL SetLineWeight(1)

        CALL PlotLine(w1x0,w1Top,w1x0,w1y0,white)
    END SUB


    ! ----- eigen values and directions ----

    SUB FindLamdas                               ! eigen values
        LET dscrm= fDsc(trace,dterm)             ! find discriminant
        IF dscrm>=0 then LET sqrtDscrm= sqr(dscrm)

        SELECT CASE dscrm
        CASE 0                                   ! one real root
             ! Case  0: dscrm=0  on parabola - one root - lam1=lam2
             LET lam1,lam2= trace/2
             CALL SetStrings(using$("-%.##",lam1),using$("-%.##",lam2),l1$,l2$)
        CASE is > 0                              ! two real roots
             ! Case  1: dscrm>0  outside parabola - two roots - lam1<lam2
             LET lam1= (trace+sqrtDscrm)/2
             LET lam2= (trace-sqrtDscrm)/2
             CALL SetStrings(using$("-%.##",lam1),using$("-%.##",lam2),l1$,l2$)
        CASE is < 0                              ! complex roots
             ! Case -1: dscrm<0  inside parabola - complex
             LET lam1= trace/2
             LET lam2= sqr(-dscrm)/2
             CALL SetStrings(trim$(using$("-%.##",lam1)),trim$(using$("-%.##i",abs(lam2))),r$,i$)
             CALL SetStrings(r$ & " - " & i$, r$ & " + " & i$,l1$,l2$)
        END SELECT
    END SUB

    SUB FindEigenVectors(ma,mb,mc,md)
        LET places= 5
        LET rma= round(ma,places)                ! round the matrix
        LET rmb= round(mb,places)
        LET rmc= round(mc,places)
        LET rmd= round(md,places)
        IF dscrm>=0 then LET sqrtDscrm= sqr(dscrm)

        IF dscrm=0 then                          ! on parabola - one real root
           IF rma=rmd then
              IF rmc=0 and rmb=0 then            ! star
                 ! no eigenlines
              ELSE IF rmb<>0 and rmc=0 then      ! degenerate node
                 CALL CheckEigenPoint(1,0)       ! horizontal
              ELSE IF rmb=0 and rmc<>0 then      ! degenerate node
                 CALL CheckEigenPoint(0,1)       ! vertical
              END IF
           ELSE IF rmc=0 and rmd=0 then          ! first row points needed
              LET x= 2*rmb
              LET y= rmd-rma-sqrtDscrm
              CALL CheckEigenPoint(x,y)
           ELSE IF rmc=0 and rma=1 and rmb=1 and rmd=1 then
              CALL CheckEigenPoint(1,0)          ! horizontal
           ELSE                                  ! otherwise second row points
              LET x= rma-rmd-sqrtDscrm
              LET y= 2*rmc
              CALL CheckEigenPoint(x,y)
           END IF

        ELSE IF dscrm>0 then                     ! two real roots
           IF rmc=0 and rmd=0 then               ! first row points needed?
              LET x= 2*rmb

              LET y= rmd-rma-sqrtDscrm
              CALL CheckEigenPoint(x,y)

              LET y= rmd-rma+sqrtDscrm
              CALL CheckEigenPoint(x,y)
           ELSE                                  ! otherwise second row points
              IF rmc<>0 then
                 LET y= 2*rmc

                 LET x= rma-rmd-sqrtDscrm
                 CALL CheckEigenPoint(x,y)

                 LET x= rma-rmd+sqrtDscrm
                 CALL CheckEigenPoint(x,y)
              ELSE
                 CALL CheckEigenPoint(1,0)       ! horizontal

                 IF rmb=0 then
                    CALL CheckEigenPoint(0,1)    ! vertical
                 ELSE
                    LET x= 1
                    LET y= (rmd-rma)/rmb
                    CALL CheckEigenPoint(x,y)
                 END IF
              END IF
           END IF
        END IF
        CALL w3KeepGraphLayer
    END SUB

    SUB CheckEigenPoint(x,y)
        CALL SetVars(round(x,8),round(y,8),ex,ey)
        IF ex=0 and ey=1 then                    ! vertical
           CALL PlotLine(w3x0,w3Top+1,w3x0,w3Bas-1,classClr)
        ELSE IF ex=1 and ey=0 then               ! horizontal
           CALL PlotLine(w3Lft+1,w3y0,w3Rgt-1,w3y0,classClr)
        ELSE IF abs(ey)>=abs(ex) then            ! tall - scale y to hgt
           IF ey<>0 then
              LET scl= w3Size/ey
              CALL DrawEDirection(ex,ey,scl)
           ELSE
              CALL PlotLine(w3x0,w3Top+1,w3x0,w3Bas-1,classClr)
           END IF
        ELSE IF abs(ex)>abs(ey) then             ! wide - scale x to wid
           IF ex<>0 then
              LET scl= w3Size/ex
              CALL DrawEDirection(ex,ey,scl)
           ELSE
              CALL PlotLine(w3Lft+1,w3y0,w3Rgt-1,w3y0,classClr)
           END IF
        END IF
    END SUB

    SUB DrawEDirection(ex,ey,scl)
        CALL SetVars(ex*scl,ey*scl,bx,by)
        CALL w3MathToPixels(bx,by,wx1,wy1)
        CALL w3MathToPixels(-bx,-by,wx2,wy2)
        CALL PlotLine(wx1,wy1,wx2,wy2,classClr)
    END SUB

    ! --- linear classification ---

    SUB Classification(trace,dterm,dscrm,className$,classClr)
        LET oldClassName$= className$
        IF dterm=0 then                          ! horizontal axis of TD plane
           LET classClr  = pink
           LET className$= "degenerate"
           ! trace= 0: points         ! "origin" , "isolated fixed points" , "zero eigenvalues"
           ! trace<>0: parallel lines ! "nonisolated fixed points" , "multiple equilibrium points" , "zero eigenvalue"

        ELSE IF dterm<0 then                     ! saddle
           LET classClr  =  blue
           LET className$= "saddle"

        ELSE IF dterm>0 then                     ! upper half plane
           IF dscrm=0 then                       ! discriminate=0... on the parabola
              LET classClr= yellow
              IF ma=md and mb=0 and mc=0 then    ! star
                 IF trace>0 then
                    LET className$= "star source"
                 ELSE IF trace<0 then
                    LET className$= "star sink"
                 END IF
              ELSE                               ! node
                 ! LET className$= "improper node"
                 LET className$= "defective node"
              END IF
              ! "tr^[2] = 4 det" , "repeated eigenvalue"

           ELSE IF trace<0 then                  ! left half - stable
              IF dscrm>0 then                    ! complex - within parab
                 LET classClr  =  green
                 LET className$= "nodal sink"
                 ! "stable node"
              ELSE
                 LET classClr  =  red
                 LET className$= "spiral sink"
                 ! "stable spiral"
              END IF

           ELSE IF trace>0 then                  ! right half - unstable
              IF dscrm>0 then                    ! complex
                 LET classClr  =  magenta
                 LET className$= "nodal source"
                 ! "unstable node"
              ELSE
                 LET classClr  =  cyan
                 LET className$= "spiral source"
                 ! "unstable spiral"
              END IF

           ELSE IF trace=0 then                  ! neutral center
              LET classClr  =  white
              LET className$= "center"
           END IF
        END IF

        IF oldClassName$<>className$ then CALL ClassText(className$)
    END SUB


    ! --------- direction and vector fields ------------

    SUB DirectionField
        CALL w3ShowGridLayer
        SET COLOR drkmid
        LET vstp= 30
        LET hstp= vstp/2
        FOR wx= w3Lft+hstp to w3Rgt-hstp step vstp    ! draw vector field
            FOR wy= w3Top+hstp to w3Bas-hstp step vstp
                LET x = w3fncx(wx)
                LET y = w3fncy(wy)
                LET dx= dxdt(x,y)
                LET dy= dydt(x,y)
                IF dx<>0 or dy<>0 then
                   LET ang= angle(dx,-dy)
                   LET sa = 5*sin(ang)
                   LET ca = 5*cos(ang)
                   LET wx0= wx+ca
                   LET wy0= wy+sa
                   LET wx2= wx-ca
                   LET wy2= wy-sa
                   PLOT wx2,wy2; wx0,wy0
                   CALL VectorHead(ang,wx2,wy2,wx0,wy0)
                ELSE
                   PLOT wx,wy
                END IF
            NEXT wy
        NEXT wx
        CALL w3KeepFieldLayer
    END SUB

    SUB VectorHead(ang,wx,wy,wx0,wy0)
        CALL SetMat(ang)
        CALL Rotate(6, 2,x1,y1)
        CALL Rotate(6,-2,x2,y2)
        LET wx1= wx+x1
        LET wy1= wy+y1
        LET wx2= wx+x2
        LET wy2= wy+y2
        PLOT wx1,wy1; wx0,wy0; wx2,wy2
    END SUB

    SUB Rotate(xin,yin,xout,yout)
        LET xout= xin*vma + yin*vmc
        LET yout= xin*vmb + yin*vmd
    END SUB

    SUB SetMat(ang)
        LET vma=  cos(ang)
        LET vmb=  sin(ang)
        LET vmc= -vmb
        LET vmd=  vma
    END SUB

    ! ----- Rollover vector subs -----

    SUB SetVector(x1,y1,x2,y2,tmstp)
        CALL DEcoords(x1,y1,tmstp,dx,dy,x2,y2)
        LET head= 0
        IF w3fWithin(x2,y2)=true then LET head= 1
        CALL Clip(w3fLft,w3fRgt,w3fBas,w3fTop,x1,y1,x2,y2)      ! (x1,y1) is within
        CALL w3MathToPixels(x1,y1,wx1,wy1)
        CALL w3MathToPixels(x2,y2,wx2,wy2)
        CALL PlotLine(wx1,wy1, wx2,wy2, yellow)
        IF head=1 and (dx<>0 or dy<>0) then
           CALL DrawVectorArrow(angle(dx,-dy),6,wx2,wy2,yellow)
        END IF
    END SUB

    ! ---- end of vector field routines ---


    ! ---- trajectory and RK4 -----

    SUB Trajectory(x0,y0)
        SET COLOR classClr
        PLOT
        LET oldstep= tstep

        IF dscrm<0 then LET tstep= 1/256

        FOR n= 1 to -1 step -2
            LET oldFlag= 1
            LET stp    = n*tstep
            CALL SetVars(x0,y0,x,y)
            CALL w3MathToPixels(x,y,wx,wy)

            FOR i= 1 to 40/tstep                 ! draw trajectory
                CALL SetVars(wx,wy,oldwx,oldwy)
                CALL w3MathToPixels(x,y,wx,wy)
                IF w3wWithin(wx,wy)=true then
                   IF oldFlag=1 then PLOT wx,wy; oldwx,oldwy
                   LET oldFlag= 1
                ELSE
                   LET oldFlag= 0
                END IF

                CALL RungeKutta4(x,y,dx,dy,stp)

                ! equilibrium or off plane?
                IF (abs(dx)<.02 and abs(dy)<.02) or abs(x)>16 or abs(y)>16 then
                   EXIT FOR
                END IF

                IF trace=0 and dterm>0 and i>100 and abs(x-x0)<.05 and abs(y-y0)<.05 then
                   PLOT wx,wy; w3Wndx(x0),w3Wndy(y0)
                   EXIT FOR
                END IF
            NEXT i
            PLOT
            IF trace=0 and dterm>0 then EXIT FOR
        NEXT n
        LET tstep= oldstep
    END SUB

    ! ---- Runge Kutta ----

    SUB RungeKutta4(x,y,dx,dy,tstp)
        LET hstp= .5*tstp

        LET dx1= dxdt(x,y)
        LET dy1= dydt(x,y)
        LET x1 = x + dx1*hstp
        LET y1 = y + dy1*hstp

        LET t1 = t + hstp

        LET dx2= dxdt(x1,y1)
        LET dy2= dydt(x1,y1)
        LET x2 = x1 + dx2*hstp
        LET y2 = y1 + dy2*hstp

        LET t2 = t + hstp

        LET dx3= dxdt(x2,y2)
        LET dy3= dydt(x2,y2)
        LET x3 = x2 + dx3*tstp
        LET y3 = y2 + dy3*tstp

        LET dx4= dxdt(x3,y3)
        LET dy4= dydt(x3,y3)

        LET dx = (dx1 + 2*dx2 + 2*dx3 + dx4) / 6
        LET dy = (dy1 + 2*dy2 + 2*dy3 + dy4) / 6

        LET x  = x + tstp*dx
        LET y  = y + tstp*dy
    END SUB

    SUB DEcoords(inx,iny,tmstp,dx,dy,outx,outy )
        CALL DEchange(inx,iny,dx,dy)
        LET outx= inx + dx*tmstp
        LET outy= iny + dy*tmstp
    END SUB

    SUB DEchange(inx,iny,dx,dy)
        LET dx= dxdt(inx,iny)
        LET dy= dydt(inx,iny)
    END SUB


    ! ----- other drawing routines -----

    SUB MarkTDCursor(trace,dterm)
        CALL w1ShowGridLayer
        CALL w1MathToPixels(trace,dterm,wt,wd)   ! mark the spot
        IF w1wWithin(wt,wd)=true then CALL DrawCursor(wt,wd)
    END SUB

    SUB MarkTSCursor(theta,s)
        CALL w2ShowGraphLayer
        CALL w2MathToPixels(theta,s,wt,ws)
        IF w2wWithin(wt,ws)=true then CALL DrawCursor(wt,ws)
    END SUB

    SUB DrawCursor(cx,cy)
        CALL BoxClear(cx-3,cx+3,cy+1,cy-1)
        CALL BoxClear(cx-1,cx+1,cy+3,cy-3)
        CALL PlotLine(cx-3,cy,cx+3,cy,white)
        CALL PlotLine(cx,cy-3,cx,cy+3,white)
    END SUB
END SUB

! ----- end of LinearPhase Cursor code -----