!! File: Initial Conditions
!! May 26, 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$= "Initial Conditions"

SUB ThisProgram
    CALL InitialConditions
    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 ------
! -----------------------------------------------------------


! ***

SUB InitialConditions
    DECLARE PUBLIC black,drkgry,drkmid,midgry,litmid,litgry,white
    DECLARE PUBLIC red,yellow,green,cyan,blue,magenta
    DECLARE PUBLIC axislabelclr,slideclr
    DECLARE PUBLIC workLft,workRgt,workBas,workTop,workMidx
    DECLARE PUBLIC true,false

    DECLARE DEF quitWithin,infoWithin

    ! --- help screen array ---

    DIM  info$(1:1)
    MAT READ info$

    DATA "Initial Conditions"


    ! --- color definition and adjustment ---

    LET pdclr = green
    LET ivclr = yellow
    LET fldclr= drkmid
    LET posclr= yellow
    LET velclr= cyan
    LET trclr = blue

    SET COLOR MIX(blue) .4,.6,1

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

    DECLARE DEF clamp,roundn,e

    ! ---- Functions ----

    LET dt= 1/64
    DEF dxdt(p,v)=  v
    DEF dvdt(p,v)= -b*v - k*p + cos(omega*t)

    ! I have three signals in mind
    ! 0 ,  cos(omega t) ,  and  sq(omega t)

    ! There are two basic solutions to the "homogeneous equation"
    !   x" + bx' + kx = 0
    !   y(t)  with  y(0) = 1  and  y'(0) = 0
    !   z(t)  with  z(0) = 0  and  z'(0) = 1

    SUB SetW(b,k, w)
        LET w= sqr(abs((b/2)^2 - k))
    END SUB

    SUB SetCritDamp(k,critDamp)
        LET critDamp= sqr(4*k)
    END SUB

    SUB SetNatFreq(k,natFreq)
        LET natFreq= sqr(k)
    END SUB

    SUB SetD(b,k,omega, D)
        LET D= (k-omega^2)^2 + omega^2*b^2
    END SUB

    SUB SetRS(b,w, r,s)
        LET r= -b/2 + w
        LET s= -b/2 - w
    END SUB

    SUB SetVars
        CALL SetW(b,k, w)
        CALL SetD(b,k,omega, D)
        CALL SetRS(b,w, r,s)
        CALL SetCritDamp(k, critDamp)
        LET b2= b*b
        LET k4= k*4
    END SUB

    SUB InitialVals(b,k,omega,D, fxp0,dxp0)
        LET fxp0=  (k-omega^2)/D
        LET dxp0= -(omega^2*b)/D
    END SUB

    SUB InitialDifs
        LET deltax= fx0-fxp0
        LET deltad= dx0-dxp0
    END SUB

    DEF y(t)                                     ! velocity ?
        IF b>critDamp then
           LET y= 1/(s-r) * (s*e^(r*t) - r*e^(s*t))
        ELSE IF b<critDamp then
           LET y= e^(-b*t/2)*(cos(w*t) + (b/(2*w))*sin(w*t))
        ELSE
           LET y= (1 + b*t/2) * e^(-b*t/2)
        END IF
    END DEF

    DEF z(t)                                     ! acceleration ?
        IF b>critDamp then
           LET z= -1/(s-r) * (e^(r*t) - e^(s*t))
        ELSE IF b<critDamp then
           LET z= 1/w * e^(-b*t/2) * sin(w*t)
        ELSE
           LET z= t * e^(-b*t/2)
        END IF
    END DEF

    DEF xp(t)                                    ! periodic solution for f(t)= cos(wt)
        LET xp= 1/D * ( (k-omega^2)*cos(omega*t) - (omega*b)*sin(omega*t) )
    END DEF

    DEF tr(t)                                    ! transient
        LET tr= deltax*y(t) + deltad*z(t)
    END DEF

    DEF fx(t)                                    ! solution
        LET fx= tr(t) + xp(t)
    END DEF


    ! Case:  f(t) = 0 (There's no frequency omega in this case). 
    ! If the initial condition on the phase plane reads (x(0),x'(0)),
    ! then the graph plots:
    !   x(t) = x(0)*y(t) + x'(0)*z(t)


    ! Case:  f(t) = cos(omega t)
    ! Then the periodic solution is given by

    !   D = (k-omega^2)^2 + omega^2*b^2
    !   xp(t) = 1/D * ( (k-omega^2)*cos(omega*t) - (omega*b)*sin(omega*t) )

    ! The initial value is

    !   xp(0)  = (k - omega^2)/D
    !   xp'(0) = -(omega^2*b)/D

    ! and so when this menu choice is opened xp(t) should be drawn 
    ! and this cursor point appears on the phase plane. 
    !
    ! When the user moves to the point (x(0),x'(0)) on the phase plane,
    ! the solution becomes  

    !   x(t) = xp(t) + (x(0)-xp(0))*y(t) + (x'(0)-xp'(0))*z(t)


    SUB RungeKutta4(x,v)
        LET dx1= dxdt(x,v)
        LET dv1= dvdt(x,v)
        LET x1 = x + .5*dx1*dt
        LET v1 = v + .5*dv1*dt
        LET dx2= dxdt(x1,v1)
        LET dv2= dvdt(x1,v1)
        LET x2 = x + .5*dx2*dt
        LET v2 = v + .5*dv2*dt
        LET dx3= dxdt(x2,v2)
        LET dv3= dvdt(x2,v2)
        LET x3 = x + dx3*dt
        LET v3 = v + dv3*dt
        LET dx4= dxdt(x3,v3)
        LET dv4= dvdt(x3,v3)

        LET dv = (dv1 + 2*dv2 + 2*dv3 + dv4) / 6
        LET dx = (dx1 + 2*dx2 + 2*dx3 + dx4) / 6
        LET v  = v + dt*dv
        LET x  = x + dt*dx
    END SUB

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

    LET wsize = 160
    LET fsize =   4
    LET whgt  = 200

    ! --- plane 1 data ---

    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= 1
    LET w1xMult  = pi

    LET w1yPiFlag= 0
    LET w1yMult  = 1

    LET w1Lft= workRgt - 440                     ! pixel bounds
    LET w1Rgt= w1Lft + 400                       !360
    LET w1Top= workTop + 120
    LET w1Bas= w1Top + whgt

    LET w1fLft=   0                              ! function bounds
    LET w1fRgt=  10
    LET w1fTop=  fsize
    LET w1fBas= -fsize

    LET w1xAx$= "t"                              ! axis labels
    LET w1yAx$= "x"

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

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

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

    ! --- Plane 1 methods ---

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

    CALL w1Variables

    SUB w1Init
        BOX CLEAR w1Lft-40,w1Lft-4,w1Bas+8,w1Top-8
        CALL w1DrawPlane(1,1,1)                  ! x axis, y axis, zeroaxes

        CALL SetTextFont(1,12,"bold")
        CALL PlotTextLJ(w1Rgt+8,w1y0+3,w1xAx$,axislabelclr)     ! axis labels

        LET  ybas= w1Top-11
        CALL PlotTextCJ(w1x0,ybas,w1yAx$,yellow)

        CALL w1KeepGridLayer
    END SUB


    ! ------------------------------------------
    ! --- plane 3 data ---

    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 = worklft+50                       ! pixel bounds
    LET w3Rgt = w3Lft+whgt
    LET w3Top = w1Top
    LET w3Bas = w3Top+whgt

    LET w3fLft= -fsize                           ! function bounds
    LET w3fRgt=  fsize
    LET w3fTop=  fsize
    LET w3fBas= -fsize

    LET w3xAx$= "x(0)"                           ! axis labels
    LET w3yAx$= "x(0)"

    LET w3xGridstep= 0                           ! grid line intervals
    LET w3yGridstep= 0

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

    LET w3ySTik = 1
    LET w3yLTik = 2
    LET w3yLabel= 2
    LET w3yFirst= w3fBas


    ! --- plane 3 methods ---

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

    CALL w3Variables

    SUB w3Init
        CALL w3DrawPlane(1,1,1)                  ! grid, axes, zeroaxes
        CALL SetTextFont(1,12,"bold")

        CALL PlotTextLJ(w3Rgt+9,w3y0+3,w3xAx$,cyan)   ! axis labels
        CALL PlotTextCJ(w3Rgt+13,w3y0-7,".",cyan)

        CALL PlotTextCJ(w3x0,w3Top-11,w3yAx$,yellow)

        CALL w3KeepGridLayer
        CALL w3KeepGraphLayer
    END SUB


    ! ---------- Slider parameters and methods ----------

    LET slideWidth= 200

    ! --- Slider 2 ---

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

    DECLARE DEF h2Within                         ! window/function transforms

    LET h2Flag  = 0
    LET h2PiAxis= 0
    LET h2Mult  = 1

    LET h2clr   = slideclr
    LET h2name$ = "b"
    LET h2form$ = "-%.##"
    LET h2Places= 2

    LET h2axis  = w3Bas + 90
    LET h2wLft  = w3Lft
    LET h2wRgt  = h2wLft+slidewidth

    LET h2fLft  = 0
    LET h2fRgt  = 1                              ! needs 4 if critical damping is on

    LET h2STik  = 0.1
    LET h2LTik  = 0.5
    LET h2Label = 1
    LET h2First = h2fLft
    LET h2Click = 0.1

    CALL h2SliderVariables

    DECLARE DEF h2Wndx

    ! --- slider 3 ---

    DECLARE PUBLIC h3axis,h3wLft,h3wRgt,h3wBas,h3wTop,h3fLft,h3fRgt
    DECLARE PUBLIC h3name$,h3form$,h3clr,h3First,h3STik,h3LTik,h3Label
    DECLARE PUBLIC h3PiAxis,h3Mult,h3fMin,h3fMax

    DECLARE DEF h3Within                         ! window/function transforms

    LET h3Flag  = 0
    LET h3PiAxis= 0
    LET h3Mult  = 1

    LET h3clr   = slideclr
    LET h3name$ = "k"
    LET h3form$ = "-%.##"
    LET h3Places= 2

    LET h3axis  = h2axis+60
    LET h3wLft  = h2wLft
    LET h3wRgt  = h3wLft+slideWidth

    LET h3fLft  = 0
    LET h3fRgt  = 2.0

    LET h3STik  = 0.1
    LET h3LTik  = 0.5
    LET h3Label = 1
    LET h3First = h3fLft
    LET h3Click = 0.1

    CALL h3SliderVariables

    ! --- slider 5 ---

    DECLARE PUBLIC h5axis,h5wLft,h5wRgt,h5wBas,h5wTop,h5fLft,h5fRgt
    DECLARE PUBLIC h5name$,h5form$,h5clr,h5First,h5STik,h5LTik,h5Label
    DECLARE PUBLIC h5PiAxis,h5Mult,h5fMin,h5fMax

    DECLARE DEF h5Within                         ! window/function transforms

    LET h5Flag  = 0
    LET h5PiAxis= 0
    LET h5Mult  = 1

    LET h5clr   = slideclr
    LET h5name$ = ""
    LET h5form$ = "-%.##"
    LET h5Places= 2

    LET h5axis  = h3axis
    LET h5wLft  = w1Lft
    LET h5wRgt  = h5wLft+slidewidth

    LET h5fLft  = 0
    LET h5fRgt  = 2

    LET h5STik  = 0.1
    LET h5LTik  = 0.5
    LET h5Label = 1
    LET h5First = h5fLft
    LET h5Click = 0.1

    CALL h5SliderVariables

    SUB h5Init
        CALL h5DrawSlider(h5name$,omega)
        CALL SwapOmega(h5wLft-17,h5wBas-3,"w","w",litgry)
    END SUB

    DECLARE DEF h5Wndx

    ! --- check boxes ---
    ! --- field switch ---

    DECLARE PUBLIC chkSiz

    LET cb1Lft = w3Lft
    LET cb1Rgt = cb1Lft + chkSiz
    LET cb1Bas = w3Bas + 40
    LET cb1Top = cb1Bas - chkSiz
    LET cb1Txt$= "vector field"
    LET cb1Clr = litgry


    ! --- Text Output Rects ----------
    ! --- text rectangle 1 ---

    LET t1BasLn= w1Top - 50
    LET t1txt$ = "x + bx + kx = cos(wt)"
    CALL StringWidth(t1txt$,sw)

    LET t1Lft  = w1Midx - sw/2
    LET t1Rgt  = t1Lft + sw + 20
    LET t1Bas  = t1BasLn +  5
    LET t1Top  = t1BasLn - 15

    LET t1Label$= ""

    SUB t1Label
        CALL SuperSubScriptRJ(t1Lft,t1BasLn,t1Label$,white)
    END SUB

    SUB t1Set
        CALL SetTextFont(1,12,"bold")
        LOCAL r$,i$,z$
        CALL t1Clear
        CALL SwapOmega(t1Lft,t1BasLn,t1txt$,"w",yellow)

        LET dot2= t1Lft + 4
        CALL PlotTextCJ(dot2,t1BasLn-10,"..",yellow)

        CALL StringWidth("x + b",sw)
        LET dot1= t1Lft + sw + 4
        CALL PlotTextCJ(dot1,t1BasLn-10,".",yellow)
    END SUB

    SUB t1Clear
        BOX CLEAR t1Lft-2,t1Rgt,t1Bas,t1Top
    END SUB

    SUB t1Init
        ! CALL t1Label
        CALL t1Set
    END SUB


    ! --- text rectangle 2 ---

    LET t2BasLn1 = w3Top - 70
    LET t2BasLn2 = t2BasLn1 + 20
    LET t2Lft    = w3Lft - 8
    LET t2Rgt    = t2Lft + 99
    LET t2Top    = t2BasLn1 - 15
    LET t2Bas    = t2BasLn2 +  5
    LET t2Label1$= "x(0) = "                     ! x
    LET t2Label2$= "x(0) = "                     ! x'
    CALL StringWidth(t2Label1$,sw)
    LET t2eqx    = t2Lft+sw
    LET t2form$  = "--%.##"
    CALL StringWidth(t2form$,sw)
    LET t2tRgt   = t2eqx + sw

    SUB t2Label
        CALL SetTextFont(1,12,"bold")
        CALL PlotTextRJ(t2eqx,t2BasLn1,t2Label1$,pdclr)
        CALL PlotTextRJ(t2eqx,t2BasLn2,t2Label2$,pdclr)
        CALL PlotTextCJ(t2Lft+4,t2BasLn2-10,".",pdclr)
    END SUB

    SUB t2Set
        CALL SetTextFont(1,12,"bold")
        CALL t2ClearVal
        CALL PlotTextRJ(t2tRgt,t2BasLn1,using$(t2form$,fxp0),pdclr)
        CALL PlotTextRJ(t2tRgt,t2BasLn2,using$(t2form$,dxp0),pdclr)
    END SUB

    SUB t2ClearVal
        BOX CLEAR t2eqx-6,t2Rgt,t2Bas,t2Top
    END SUB
    SUB t2Clear
        BOX CLEAR t2Lft-2,t2Rgt,t2Bas,t2Top
    END SUB

    SUB t2Init
        CALL t2Label
        CALL t2Set
    END SUB

    ! --- text rectangle 3 ---

    LET t3BasLn1 = w3Top - 70
    LET t3BasLn2 = t3BasLn1 + 20
    LET t3Lft    = w3Lft + 115
    LET t3Rgt    = t3Lft + 99
    LET t3Top    = t3BasLn1 - 15
    LET t3Bas    = t3BasLn2 +  5
    LET t3Label1$= "x(0) = "                     ! x
    LET t3Label2$= "x(0) = "                     ! x'
    CALL StringWidth(t3Label1$,sw)
    LET t3eqx    = t3Lft+sw
    LET t3form$  = "--%.##"

    CALL StringWidth(t3form$,sw)
    LET t3tRgt= t3eqx + sw

    SUB t3Label
        CALL SetTextFont(1,12,"bold")
        CALL PlotTextRJ(t3eqx,t3BasLn1,t3Label1$,posclr)
        CALL PlotTextRJ(t3eqx,t3BasLn2,t3Label2$,velclr)
        CALL PlotTextCJ(t3Lft+4,t3BasLn2-10,".",velclr)
    END SUB

    SUB t3Set
        CALL SetTextFont(1,12,"bold")
        CALL t3ClearVal
        CALL PlotTextRJ(t3tRgt,t3BasLn1,using$(t3form$,fx0),posclr)
        CALL PlotTextRJ(t3tRgt,t3BasLn2,using$(t3form$,dx0),velclr)
    END SUB

    SUB t3Clear
        BOX CLEAR t3Lft-2,t3Rgt,t3Bas,t3Top
    END SUB
    SUB t3ClearVal
        BOX CLEAR t3eqx-6,t3Rgt,t3Bas,t3Top
    END SUB

    SUB t3Init
        CALL t3Label
        CALL t3Set
    END SUB

    !    ! --- text rectangle 4 ---
    !
    !    LET t4BasLn1 = w1Bas + 45
    !    LET t4BasLn2 = t4BasLn1 + 20
    !    LET t4BasLn3 = t4BasLn2 + 20
    !    LET t4Label1$= "Steady State"                ! x
    !    LET t4Label2$= "Transient"                   ! x'
    !    LET t4Label3$= "Solution"                    ! x'
    !    CALL StringWidth(t4Label1$,sw)
    !    LET t4Lft= w1Rgt - sw
    !    LET t4Rgt= w1Rgt + 2
    !    LET t4Top= t4BasLn1 - 15
    !    LET t4Bas= t4BasLn3 +  5
    !    LET t4eqx= t4Lft+sw
    !
    !    SUB t4Label
    !        CALL SetTextFont(1,12,"bold")
    !        CALL PlotTextRJ(t4eqx,t4BasLn1,t4Label1$,pdclr)
    !        CALL PlotTextRJ(t4eqx,t4BasLn2,t4Label2$,trclr)
    !        CALL PlotTextRJ(t4eqx,t4BasLn2,t4Label3$,psclr)
    !    END SUB
    !
    !    SUB t4Set
    !        CALL SetTextFont(1,12,"bold")
    !        CALL t4ClearVal
    !        CALL PlotTextLJ(t4eqx,t4BasLn1,using$("--%.##",fx0),posclr)
    !        CALL PlotTextLJ(t4eqx,t4BasLn2,using$("--%.##",dx0),velclr)
    !    END SUB
    !
    !    SUB t4Clear
    !        BOX CLEAR t4Lft-2,t4Rgt,t4Bas,t4Top
    !    END SUB
    !    SUB t4ClearVal
    !        BOX CLEAR t4eqx,t4Rgt,t4Bas,t4Top
    !    END SUB
    !
    !    SUB t4Init
    !        CALL SetTextFont(1,12,"bold")
    !        CALL PlotTextRJ(t4eqx,t4BasLn1,t4Label1$,pdclr)
    !        CALL PlotTextRJ(t4eqx,t4BasLn2,t4Label2$,trclr)
    !        CALL PlotTextRJ(t4eqx,t4BasLn3,t4Label3$,ivclr)
    !    END SUB


    ! --- checkboxes ---

    DECLARE PUBLIC c1cnt, c1stp, c1siz
    DECLARE PUBLIC c1lft, c1rgt, c1bas, c1top
    DECLARE PUBLIC c1NameList$(), c1ColorList(), c1Switch()

    CALL SetTextFont(1,12,"bold")
    CALL StringWidth("Steady State",sw)

    LET c1Lft= w1Rgt - sw - 25
    LET c1Top= w1Bas + 45
    LET c1cnt= 3

    CALL c1Variables

    MAT redim c1NameList$(1:c1cnt)
    MAT READ c1NameList$
    DATA "Steady State","Transient","Solution"

    MAT redim c1ColorList(1:c1cnt)
    MAT READ c1ColorList
    DATA 25,27,24 !9,11,8

    MAT redim c1Switch(1:c1cnt)
    MAT READ c1Switch
    DATA 1,1,1

    DECLARE DEF c1Within

    SUB SetSwitchVars
        LET ssFlag= c1Switch(1)
        LET trFlag= c1Switch(2)
        LET slFlag= c1Switch(3)
    END SUB

    ! --- Buttons ---

    SUB Buttons
        LET btnhgt= 18
        CALL SetTextFont(1,12,"bold")
        LET cllft = w1Lft
        LET clrgt = cllft + 130
        LET cltop = w1Bas +  30
        LET clbas = cltop + btnhgt
        CALL DrawButton(cllft,clrgt,clbas,cltop,5,"Clear Transient")
    END SUB

    ! --- default parameters ---

    LET b,oldb = 0.5
    LET k,oldk = 1.0
    LET dx0,fx0= 0
    LET omega,oldomega= 1
    LET vFieldFlag= 0
    LET critMarks = 0
    CALL SetSwitchVars
    CALL SetVars

    ! --- Draw the screen ---

    CALL InitScreen

    SUB InitScreen
        BOX CLEAR worklft,workrgt,workbas,worktop

        CALL w1Init
        CALL w3Init

        CALL t1Init
        CALL t2Init
        CALL t3Init

        CALL c1DrawCheckBoxes
        CALL c1SetCheckBox

        CALL h2DrawSlider(h2name$,b)
        IF critMarks=1 then
           CALL MarkCritDamp(k)
           CALL MarkNatFreq(k)
        END IF

        CALL h3DrawSlider(h3name$,k)
        CALL h5Init

        IF vFieldFlag=1 then
           CALL CheckBoxDraw(cb1Lft,cb1Bas,cb1Txt$,cb1Clr)
           CALL CheckBoxSet(cb1Lft,cb1Bas,vFieldFlag,cb1Clr)
        END IF

        CALL w1SliderGraphs
    END SUB

    ! ----------------- Event manager -----------------

    DO
       IF ms<>2 then
          CALL MouseDown(mx,my,ms)
       END IF

       ! mouse down event
       ! window

       IF w3Within(mx,my)=true then
          LET oldmx,oldmy= -999
          !          CALL MouseUp(mx,my,ms)
          !          CALL Orbit(fx0,dx0)

          CALL w1xpGraph
          DO
             GET MOUSE: mx,my,ms
             IF mx<>oldmx or my<>oldmy then
                CALL w3wClamp(mx,my)
                CALL w3PixelsToMath(mx,my,dx0,fx0)
                CALL w1PlaneGraphs
                CALL CopyCoords(mx,my,oldmx,oldmy)
             END IF
          LOOP until ms=3

          ! sliders

       ELSE IF h2Within(mx,my)=true then         ! b
          IF my<h2axis then
             CALL h2MouseClick
          ELSE
             CALL h2MouseDrag
          END IF
       ELSE IF h3Within(mx,my)=true then         ! k
          IF my<h3axis then
             CALL h3MouseClick
          ELSE
             CALL h3MouseDrag
          END IF
       ELSE IF h5Within(mx,my)=true then         ! w
          IF my<h5axis then
             CALL h5MouseClick
          ELSE
             CALL h5MouseDrag
          END IF

          ! buttons

       ELSE IF c1Within(mx,my)=true then
          CALL MouseUp(mx,my,ms)
          LET cb= int((my-c1Top)/c1stp) + 1
          LET c1Switch(cb)= 1-c1Switch(cb)
          CALL c1SetCheckBox
          CALL SetSwitchVars
          IF ssFlag=0 then
             CALL ShowGridLayers
             CALL KeepGraphLayers
             CALL t2ClearVal
          END IF
          IF trFlag=0 and slFlag=0 then CALL t3ClearVal
          CALL w1SliderGraphs

       ELSE IF mx>cb1Lft and mx<cb1Rgt and my>cb1Top and my<cb1Bas then   ! check box
          IF vFieldFlag=true then
             CALL CheckBoxMouse(cb1Lft,cb1Bas,vFieldFlag,cb1Clr,ms)
             CALL VectorField
          END IF

       ELSE IF mx>=cllft and mx<=clrgt and my>=cltop and my<=clbas then   ! clear
          CALL MouseButtonUp(clLft,clRgt,clBas,clTop,ms)
          CALL ShowGraphLayers

       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


    ! --- Mouse Event Methods ---

    SUB ShowGridLayers
        CALL w1ShowGridLayer
        CALL w3ShowGridLayer
    END SUB
    SUB KeepGraphLayers
        CALL w1KeepGraphLayer
        CALL w3KeepGraphLayer
    END SUB
    SUB ShowGraphLayers
        CALL w1ShowGraphLayer
        CALL w3ShowGraphLayer
    END SUB

    SUB MarkCritDamp(k)
        CALL SetCritDamp(k, critDamp)
        LET wcd  = h2Wndx(critDamp)
        LET cdTop= h2wBas+2
        CALL ValueMarker(h2wLft,h2wRgt,wcd,cdtop,litgry)
    END SUB

    SUB MarkNatFreq(k)
        CALL SetNatFreq(k, natFreq)
        LET wcd  = h5Wndx(natFreq)
        LET cdTop= h5wBas+2
        CALL ValueMarker(h5wLft,h5wRgt,wcd,cdtop,litgry)
    END SUB

    SUB ValueMarker(wLft,wRgt,wx,wTop,clr)
        SET COLOR litgry
        LET wBas= wTop+6
        BOX CLEAR wLft-5,wRgt+5,wBas,wTop
        PLOT wx,wTop; wx,wBas
        PLOT wx,wTop; wx+2,wTop+2; wx-2,wTop+2; wx,wTop
    END SUB

    ! ----- h2 slider event -----

    SUB h2MouseClick
        CALL h2GetClickVal(ms,h2Click,b)
        CALL H2Action
        CALL KeepGraphLayers
    END SUB
    SUB h2MouseDrag
        DO
           CALL h2GetDragVal(ms,h2Places,b)
           CALL h2Action
        LOOP until ms=3
        CALL KeepGraphLayers
    END SUB

    SUB h2Action
        IF b<>oldb then
           CALL w1SliderGraphs
           LET oldb= b
        END IF
    END SUB

    ! ----- h3 slider event -----

    SUB h3MouseClick
        CALL h3GetClickVal(ms,h3Click,k)
        CALL h3Action
        CALL KeepGraphLayers
    END SUB
    SUB h3MouseDrag
        DO
           CALL h3GetDragVal(ms,h3Places,k)
           CALL h3Action
        LOOP until ms=3
        CALL KeepGraphLayers
    END SUB

    SUB h3Action
        IF k<>oldk then
           IF critMarks=1 then
              CALL MarkCritDamp(k)
              CALL MarkNatFreq(k)
           END IF
           CALL w1SliderGraphs
           LET oldk= k
        END IF
    END SUB

    ! ----- h5 slider event -----

    SUB h5MouseClick
        CALL h5GetClickVal(ms,h5Click,omega)
        CALL h5Action
    END SUB
    SUB h5MouseDrag
        DO
           CALL h5GetDragVal(ms,h5Places,omega)
           CALL h5Action
        LOOP until ms=3
    END SUB

    SUB h5Action
        IF omega<>oldomega then
           IF omega<>oldomega then
              CALL w1SliderGraphs
              LET oldomega= omega
           END IF
        END IF
    END SUB


    ! ----- periodic solution -----

    SUB w1xpGraph
        CALL SetVars
        CALL ShowGridLayers
        IF ssFlag=1 then
           IF D<>0 then
              CALL InitialVals(b,k,omega,D, fxp0,dxp0)
              CALL w3MathToPixels(dxp0,fxp0,wx,wy)
              CALL t2Set
              IF w3wWithin(wx,wy)=true then
                 CALL PlotDiamondClr(wx,wy,pdclr)
              END IF

              SET COLOR pdclr
              LET oldwy= w1Wndy(xp(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = xp(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           ELSE
              CALL w1ShowGridLayer
           END IF
        ELSE
           CALL w1ShowGridLayer
        END IF
        CALL KeepGraphLayers
    END SUB

    SUB w1SliderGraphs
        CALL SetVars
        CALL ShowGridLayers

        IF D<>0 then
           CALL InitialVals(b,k,omega,D, fxp0,dxp0)
           CALL w3MathToPixels(dxp0,fxp0,wx,wy)
           CALL InitialDifs
           IF ssFlag=1 then
              CALL t2Set
              IF w3wWithin(wx,wy)=true then
                 CALL PlotDiamondClr(wx,wy,pdclr)
              END IF
           END IF
           IF slFlag=1 or trFlag=1 then
              CALL t3Set
              CALL w3MathToPixels(dx0,fx0,wx,wy)
              IF w3wWithin(wx,wy)=true then
                 CALL PlotDiamondClr(wx,wy,yellow)
              END IF
           END IF

           IF ssFlag=1 then
              SET COLOR pdclr
              LET oldwy= w1Wndy(xp(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = xp(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           END IF

           IF slFlag=1 then
              SET COLOR yellow
              LET oldwy= w1Wndy(fx(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = fx(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           END IF

           IF trFlag=1 then
              SET COLOR trclr
              LET oldwy= w1Wndy(fx(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = tr(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           END IF
        END IF
    END SUB

    SUB w1PlaneGraphs
        CALL SetVars
        CALL ShowGraphLayers

        IF D<>0 then
           IF slFlag=1 or trFlag=1 then
              CALL InitialVals(b,k,omega,D, fxp0,dxp0)
              CALL InitialDifs
              CALL t3Set
              CALL w3MathToPixels(dx0,fx0,wx,wy)
              IF w3wWithin(wx,wy)=true then
                 CALL PlotDiamondClr(wx,wy,yellow)
              END IF
           END IF

           IF slFlag=1 then
              SET COLOR yellow
              LET oldwy= w1Wndy(fx(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = fx(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           END IF

           IF trFlag=1 then
              SET COLOR trclr
              LET oldwy= w1Wndy(fx(0))
              LET oldwx= w1Lft
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = tr(t)
                  LET wy= w1Wndy(x)

                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then     ! both in plane
                     PLOT wx,wy; oldwx,oldwy
                  ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then  ! one in - clip
                     CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                     CALL CopyCoords(wx,wy,wx2,wy2)
                     CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                     PLOT wx1,wy1; wx2,wy2
                  ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then  ! above and below
                     IF wy>oldwy then
                        PLOT oldwx,w1top+1; wx,w1bas-1
                     ELSE
                        PLOT wx,w1top+1; oldwx,w1bas-1
                     END IF
                  END IF
                  CALL CopyCoords(wx,wy,oldwx,oldwy)
              NEXT wx
              PLOT
           END IF
        END IF
    END SUB

    ! ----- transient solution -----

    SUB MatchxpGraph
        !        LET fx0= fxp0
        !        LET dx0= dxp0
        CALL t3Set
        CALL w1xGraph
    END SUB

    SUB w1xGraph
        IF slFlag=1 or trFlag=1 then
           CALL SetVars
           IF D<>0 then
              CALL InitialDifs
              CALL w3MathToPixels(dx0,fx0,wx,wy)

              CALL w3ShowGraphLayer
              IF w3wWithin(wx,wy)=true then
                 CALL PlotDiamondClr(wx,wy,yellow)
              END IF

              SET COLOR yellow
              LET oldwy= w1Wndy(fx(0))
              LET oldwx= w1Lft
              CALL w1ShowGraphLayer
              FOR wx= w1Lft to w1Rgt
                  LET t = w1Fncx(wx)
                  LET x = fx(t)
                  LET wy= w1Wndy(x)

                  IF slFlag=1 then
                     SET COLOR yellow
                     IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then  ! both in plane
                        PLOT wx,wy; oldwx,oldwy
                     ELSE IF (wy>w1Top and wy<w1Bas) or (oldwy>w1Top and oldwy<w1Bas) then    ! one in - clip
                        CALL CopyCoords(oldwx,oldwy,wx1,wy1)
                        CALL CopyCoords(wx,wy,wx2,wy2)
                        CALL ClipVert(w1Bas,w1Top,wx1,wy1,wx2,wy2)
                        PLOT wx1,wy1; wx2,wy2
                     ELSE IF (wy>w1Bas and oldwy<w1Top) or (wy<w1Bas and oldwy>w1Top) then    ! above and below
                        IF wy>oldwy then
                           PLOT oldwx,w1top+1; wx,w1bas-1
                        ELSE
                           PLOT wx,w1top+1; oldwx,w1bas-1
                        END IF
                     END IF
                     CALL CopyCoords(wx,wy,oldwx,oldwy)
                  END IF
                  IF trFlag=1 then
                     SET COLOR trclr
                     LET x = tr(t)               ! transient !!
                     LET wy= w1Wndy(x)
                     IF wy>w1Top and wy<w1Bas then PLOT wx,wy
                  END IF
              NEXT wx
              PLOT
           ELSE
              CALL w1ShowGridLayer
           END IF
        END IF
    END SUB


    ! --- graphing and interface parts ---

    SUB Orbit(p0,v0)
        LET t,erase= 0
        CALL CopyCoords(p0,v0,pos,vel)
        CALL w1MathToPixels(pos,vel,wpy,wvy)
        LET wx     = w1Lft                       ! time series
        LET resFlag= 0

        DO
           CALL RungeKutta4(pos,vel)

           CALL CopyCoords(wpy,wvy,oldpy,oldvy)
           LET wpy= w1Wndy(pos)

           LET wx3= w3Wndx(vel)                  ! phase plane
           IF w3wWithin(wx3,wpy)=true then
              SET COLOR posclr
              PLOT wx3,wpy
           ELSE
              LET resflag= 1
           END IF

           LET oldx= wx
           LET wx  = w1Wndx(t)                   ! time series t
           IF wx>w1Rgt then EXIT DO
           IF wpy>w1Top and wpy<w1Bas and oldpy>w1Top and oldpy<w1Bas then
              CALL PlotLine(wx,wpy,oldx,oldpy,posclr)
           END IF

           GET MOUSE: mx,my,ms
           IF ms=2 then
              CALL MouseUp(mx,my,ms)
              EXIT DO
           END IF

           ! CALL Delay(1/50)
           LET t= t+dt
           CALL PlotPoint(wx3,wpy,phsclr)
        LOOP                                     ! until ms=2
    END SUB


    ! --- vector field routines ---

    SUB VectorField
        CALL w3ShowGridLayer
        IF f=0 and vFieldFlag=1 then
           SET COLOR fldclr
           LET stp = 20
           LET hstp= 10
           FOR wy= w3Top+hstp to w3Bas-hstp step stp
               FOR wx= w3Lft+hstp to w3Rgt-hstp step stp
                   CALL w3PixelsToMath(wx,wy,v,p)
                   LET dy= dxdt(p,v)
                   LET dx= dvdt(p,v)

                   IF dx<>0 or dy<>0 then
                      LET ang= angle(dx,-dy)
                      LET wx1= wx + 8*cos(ang)
                      LET wy1= wy + 8*sin(ang)
                      PLOT wx,wy; wx1,wy1
                      CALL VectorHead(ang,wx,wy,wx1,wy1)
                   END IF
               NEXT wx
           NEXT wy
        END IF
        CALL w3KeepGraphLayer
    END SUB

    SUB VectorHead(ang,wx,wy,wvx0,wvy0)
        CALL SetVMat(ang)
        CALL RotateTranslate(6, 2,wx,wy,wvx1,wvy1)
        CALL RotateTranslate(6,-2,wx,wy,wvx2,wvy2)
        PLOT wvx1,wvy1; wvx0,wvy0; wvx2,wvy2
    END SUB

    SUB Rotate(xin,yin,xout,yout)
        LET xout= xin*vcta + yin*vctb
        LET yout= xin*vctc + yin*vctd
    END SUB
    SUB RotateTranslate(xin,yin,h,v,xout,yout)
        LET xout= xin*vcta + yin*vctb + h
        LET yout= xin*vctc + yin*vctd + v
    END SUB

    SUB SetVMat(ang)
        LET vcta=  cos(ang)
        LET vctc=  sin(ang)
        LET vctb= -vctc
        LET vctd=  vcta
    END SUB
END SUB

! --- end of mass spring code -------------------