!   Exception 2001 caused at line 6868 in -c1mouseup; entering immediate mode
!   
!   There's a lot in this tool. It'll be fun to write problems for.
!   
!   One oddity: when I set k = omega^2 and push b to zero, 
!   the steady state curve vanishes, correctly, if I had it displayed,
!   but the solution appears whether or not I'd selected it.
!   
!   Also, in some circumstances the steady state can be selected
!   but its initial condition numerical readout is not displayed.


!    SUB w2IVpPoint                               ! steady state initial values

!  When k = omega = 0  and b is not zero, then
!  
!  x = (1/b)t + x(0) + ((x'(0)/b) - (1/b^2))*(1 - e^{-bt}).
!  
!  If all three are zero then 
!  
!  x = (t*t/2) + x(0) + x'(0) t .
!  
!  These are all cases of "resonance," which means that we don't want to
!  give meaning to the "steady state" or "transient," and those graphs
!  should vanish (as they do).

! When b = 0  and omega^2 = k ,  all the displays go off.  This is resonance. Let's be clear here.

! When b = 0, it's somewhat of a misnomer to talk about "transients," 
! since the only homogeneous solution which dies off as t grows is the zero solution. 
! Still, it does make sense to talk about "steady state," for the poor reason that 
! there is a uniform formula for the steady state when b>0 which continues to make sense when b = 0: then it is
! cos(omega t) / (omega_n^2 - omega^2).

! This is what you show as long as  omega  differs from  omega_n You are correct, 
! this formula does NOT make sense when the two are equal. 
! So neither steady state nor transient make sense in that situation, and it's correct to make them vanish. 
! You should also make the readout for the green steady state initial conditions vanish (with the green expressions  x(0) = ,  x'(0) = ).

! Nevertheless there IS a solution with the given initial conditions, and it should be drawn. 
! Here's the formula for it. Remember, this is when  b = 0 and k = omega^2:
! x(t) = (x(0) + t/(2 omega)) cos(omega t) + (x'(0)/omega) sin(omega t) .

! remove all numerical solution code?
! RungeKutta4, orbit, vector field, etc.


!! File: ForcedDampedVibes
!! August 25, 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$ = "Forced Damped Vibrations"

SUB ThisProgram
    CALL ForcedDampedVibes
    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\TBLibs\TrueCtrl.trc"         ! windows
LIBRARY "c:\TB Gold\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.trc"
!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 ForcedDampedVibes
    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:13)
    MAT READ info$

    ! DATA "Initial Conditions"

    ! initial value plane
    ! b, k, w sliders
    ! vertical bounds for time series
    ! graphs
    ! time series rollover

    DATA "Forced Damped Vibrations"
    DATA ""
    DATA "This tool supports graphical exploration of the effects of initial values, spring stiffness, damping, and forcing frequency on the transient, steady state, and solution curves for the second order ODE  x'' + bx' + kx = cos(omega t)."
    DATA ""
    DATA "Set the initial conditions, x(0) and x'(0), by clicking or dragging the mouse on the initial value plane, or by clicking or dragging the x(0) and x'(0) slider handles."
    DATA "The initial conditions for the steady state solution are displayed in green on the initial value plane."
    DATA ""
    DATA "Click or drag the slider handles to set the spring constant, k; the damping constant, b; or the forcing frequency, omega."
    DATA "Click a slider among the hash marks for values rounded to the nearest tenth."
    DATA ""
    DATA "The time series plane displays the steady state curve in green, the transient curve in blue, and the solution curve in yellow. Click the checkboxes to toggle the transient, steady state, and solution curves on or off."
    DATA "Click one of the number buttons to reset the vertical scale of the time series plane."
    DATA "Roll the mouse cursor over the time series plane to display crosshairs and the t and x values at the cursor."


    ! Help Text
    ! In this tool the ODE  x" + bx' + kx = cos(omega t)  is studied.
    ! The parameters b,  k,  and omega,  can be set by sliders at the bottom.

    ! The contents of the time series window are controlled by three buttons labeled "Steady State," "Transient," and "Solution." 
    ! The corresponding functions are displayed in green, blue, and yellow, respectively. 
    ! The vertical scale of this window can be set by buttons below it.

    ! A rollover causes crosshairs to appear, and a readout of the value of (t,x).
    ! Initial conditions are set by cursor click on the window at left, or by sliders.
    ! The initial condition of the steady state solution is indicated on that window in green.
    ! The values of these initial conditions are read out above.


    ! --- color definition and adjustment ---

    LET ssClr = green
    LET ivclr = yellow
    LET fldclr= drkmid
    LET posclr= yellow
    LET velclr= yellow
    LET trclr = blue

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

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

    DECLARE DEF clamp,roundn,e

    ! --- arrays ---

    DIM ssList(0:400)
    DIM trList(0:400)
    DIM slList(0:400)

    ! ---- 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

    ! Case: f(t) = 0.  no omega slider
    !  let tr= fx0*y9t) + dx0*z(t)
    !  let fx= tr


    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 SetRS(b,w, r,s)
        LET b2= b/2
        LET r = -b2 + w
        LET s = -b2 - w
        IF s-r<>0 then
           LET sr= 1/(s-r)
           LET c1= s/(s-r)
           LET c2= r/(s-r)
        END IF
        IF w<>0 then
           LET wr= 1/w
        END IF
    END SUB

    SUB SetD(b,k,omega, D)
        LET D= (k-omega^2)^2 + omega^2*b^2
    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 c1= k-omega^2
        LET c2= omega*b
        IF D<>0 then
           LET a1= c1/D
           LET a2= c2/D
        END IF
    END SUB

    SUB InitialVals(b,k,omega,D, fxp0,dxp0)
        IF D<>0 then
           LET fxp0=  (k-omega^2)/D
           LET dxp0=  (omega^2*b)/D
        END IF
    END SUB


    ! ----- these from transients -----

    !    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


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

    DEF z(t)                                     ! acceleration ?
        LET b2=  b/2
        LET w =  sqr(abs(b2^2 - k))
        LET r = -b2 + w
        LET s = -b2 - w
        IF s<>r then
           LET sr= 1/(s-r)
        ELSE
           ! ???
        END IF
        IF w<>0 then
           LET wr= 1/w
        ELSE
           ! ???
        END IF

        IF b>critDamp then
           LET z= -sr * (e^(r*t) - e^(s*t))
        ELSE IF b<critDamp then
           LET z=  wr * e^(-b2*t) * sin(w*t)
        ELSE
           LET z= t*e^(-b2*t)
        END IF
    END DEF

    DEF y(t)
        LET b2=  b/2
        LET w =  sqr(abs(b2^2 - k))
        LET r = -b2 + w
        LET s = -b2 - w
        IF s<>r then
           LET c1=  s/(s-r)
           LET c2=  r/(s-r)
        ELSE
           ! ???
        END IF
        IF w<>0 then
           LET wr= 1/w
        ELSE
           ! ???
        END IF

        IF b>critDamp then
           LET y = c1*e^(r*t) - c2*e^(s*t)
        ELSE IF b<critDamp then
           LET wt= w*t
           LET y = e^(-b2*t)*(cos(wt) + (wr*b2)*sin(wt))
        ELSE
           LET y = (1 + b2*t) * e^(-b2*t)
        END IF
    END DEF

    !    DEF sc1(t)                                    ! special case:  b=0 and k= omega^2
    !        ! x(t) = (x(0) + t/(2 omega)) cos(omega t) + (x'(0)/omega) sin(omega t)
    !        LET o2 = 2*omega
    !
    !        LET ot = omega*t
    !        LET sc = (fx0 + t/(2*omega))*cos(ot) + (dx0/omega)*sin(ot)
    !    END DEF

    DEF sc1(t)                                   ! special case:  b=0 and k= omega^2
        LET sc1 = fx0*cos(omega*t) + (dx0 + t/(2*omega))*sin(omega*t)
    END DEF

    DEF sc2(t)                                   ! special case:  b<>0 k=0 omega=0
        LET sc2 = ( (1/b)*t + fx0 + ((dx0/b) - (1/b^2))*(1 - e^(-b*t)) )
    END DEF

    DEF sc3(t)                                   ! special case:  b=0 k=0 omega=0
        LET sc3 = (t*t/2) + fx0 + dx0*t
    END DEF

    ! LET xp= a1*cos(ot) - a2*sin(ot)

    DEF xp(t)                                    ! periodic steady state solution for f(t)= cos(wt)
        LET o2= omega^2
        LET D = (k-o2)^2 + o2*b^2
        IF D<>0 then
           LET a1= (k-o2)/D
           LET a2= (omega*b)/D
           LET ot= omega*t
           LET xp= a1*cos(ot) + a2*sin(ot)
        ELSE
           ! ???
        END IF
    END DEF

    DEF tr(t)                                    ! transient
        LET tr= (fx0-fxp0)*y(t) + (dx0-dxp0)*z(t)
    END DEF

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

    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 =   2
    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 - 410                    ! pixel bounds
    LET w1Rgt = w1Lft + 360                      !360
    LET w1Top = workTop + 120
    LET w1Bas = w1Top + whgt

    LET w1fLft=  0                               ! function bounds
    LET w1fRgt=  6
    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
    DECLARE DEF w1wWithin,w1Within

    CALL w1Variables

    SUB w1Init
        BOX CLEAR w1Lft-40,w1Lft-4,w1Bas+8,w1Top-8
        CALL w1Reset(r1Num)
        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
        CALL w1KeepGraphLayer
    END SUB

    SUB w1KeepTempLayer
        BOX KEEP w1Lft,w1Rgt,w1Bas,w1Top in temp$
    END SUB

    SUB w1ShowTempLayer
        BOX SHOW temp$ at w1Lft,w1Bas
    END SUB

    SUB w1Reset(zoomState)
        LET w1fLft  = 0                          ! function bounds
        LET w1fRgt  = 6
        LET w1xSTik = 0.5                        ! horizontal axis Tik marks
        LET w1xLTik = 1
        LET w1xLabel= 1
        LET w1xFirst= w1fLft

        SELECT CASE zoomState
        CASE 1
             LET w1fTop  =  10
             LET w1fBas  = -10
             LET w1ySTik =   1                   ! vertical axis Tik marks
             LET w1yLTik =   2
             LET w1yLabel=   2
             LET w1yFirst= -10

             LET w1fLft  =  0                    ! function bounds
             LET w1fRgt  = 18
             LET w1xSTik =  1                    ! horizontal axis Tik marks
             LET w1xLTik =  2
             LET w1xLabel=  2
             LET w1xFirst= w1fLft
        CASE 2
             LET w1fTop  =  2
             LET w1fBas  = -2
             LET w1ySTik =  0.5                  ! vertical axis Tik marks
             LET w1yLTik =  1
             LET w1yLabel=  1
             LET w1yFirst= -2
        CASE 3
             LET w1fTop  =  0.1
             LET w1fBas  = -0.1
             LET w1ySTik =  0.01                 ! vertical axis Tik marks
             LET w1yLTik =  0.05
             LET w1yLabel=  0.05
             LET w1yFirst= -0.1
        CASE 4
             LET w1fTop  =  0.01
             LET w1fBas  = -0.01
             LET w1ySTik =  0.001                ! vertical axis Tik marks
             LET w1yLTik =  0.005
             LET w1yLabel=  0.005
             LET w1yFirst= -0.01
        CASE else

        END SELECT
        CALL w1Variables
    END SUB

    ! ------------------------------------------

    ! --- w2 plane data ---

    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= 0
    LET w2xMult  = 1
    LET w2yPiFlag= 0
    LET w2yMult  = 1

    LET w2Lft = worklft+70                       ! pixel bounds
    LET w2Rgt = w2Lft+whgt
    LET w2Top = w1Top
    LET w2Bas = w2Top+whgt

    LET w2fLft= -fsize                           ! function bounds
    LET w2fRgt=  fsize
    LET w2fTop=  fsize
    LET w2fBas= -fsize

    LET w2xAx$= "x(0)"                           ! axis labels
    LET w2yAx$= "x(0)"

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

    LET w2xSTik = 0.5                            ! axis Tik marks
    LET w2xLTik = 1
    LET w2xLabel= 1
    LET w2xFirst= w2fLft

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

    ! --- w2 plane methods ---

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

    CALL w2Variables

    SUB w2Init
        CALL w2DrawPlane(1,1,1)                  ! grid, axes, zeroaxes
        CALL SetTextFont(1,12,"bold")

        CALL PlotTextLJ(w2Rgt+9,w2y0+3,w2xAx$,yellow)      ! axis labels
        CALL PlotTextCJ(w2Rgt+13,w2y0-7,".",yellow)

        CALL PlotTextCJ(w2x0,w2Top-11,w2yAx$,yellow)

        CALL w2KeepGridLayer
        CALL w2KeepGraphLayer
    END SUB


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

    ! ----------- vertical slider -------

    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   = yellow
    LET v1name$ = ""
    LET v1form$ = "-%.##"
    LET v1Places= 2

    LET v1axis  = w2Lft - 25
    LET v1wBas  = w2Bas
    LET v1wTop  = w2Top
    LET v1fBas  = w2fBas
    LET v1fTop  = w2fTop

    LET v1STik  =  0.5
    LET v1LTik  =  1
    LET v1Label =  1
    LET v1First = -2
    LET v1Click =  0.5

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

    CALL v1SliderVariables

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


    ! --- horizontal sliders ---

    LET slideWidth= 200

    ! --- h1 Slider  ---

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

    DECLARE DEF h1Within                         ! window/function transforms

    LET h1Flag  = 0
    LET h1PiAxis= 0
    LET h1Mult  = 1

    LET h1clr   = velClr
    LET h1name$ = ""
    LET h1form$ = "-%.##"
    LET h1Places= 2

    LET h1axis  = w2Bas + 23
    LET h1wLft  = w2Lft
    LET h1wRgt  = h1wLft + slidewidth

    LET h1fLft  = w2fLft
    LET h1fRgt  = w2fRgt

    LET h1STik  = w2xSTik
    LET h1LTik  = w2xLTik
    LET h1Label = w2xLabel
    LET h1First = w2xFirst
    LET h1Click = 0.5

    CALL h1SliderVariables

    DECLARE DEF h1Wndx


    ! --- Slider h2 ---

    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  = w2Bas + 100
    LET h2wLft  = w2Lft
    LET h2wRgt  = h2wLft+slidewidth

    LET h2fLft  = 0
    LET h2fRgt  = 4

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

    CALL h2SliderVariables

    DECLARE DEF h2Wndx

    ! --- slider h3 ---

    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+65
    LET h3wLft  = h2wLft
    LET h3wRgt  = h3wLft+slideWidth

    LET h3fLft  = 0
    LET h3fRgt  = 4

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

    CALL h3SliderVariables

    ! --- slider h4 ---

    DECLARE PUBLIC h4axis,h4wLft,h4wRgt,h4wBas,h4wTop,h4fLft,h4fRgt
    DECLARE PUBLIC h4name$,h4form$,h4clr,h4First,h4STik,h4LTik,h4Label
    DECLARE PUBLIC h4PiAxis,h4Mult,h4fMin,h4fMax

    DECLARE DEF h4Within                         ! window/function transforms

    LET h4Flag  = 0
    LET h4PiAxis= 0
    LET h4Mult  = 1

    LET h4clr   = slideclr
    LET h4name$ = ""
    LET h4form$ = "-%.##"
    LET h4Places= 2

    LET h4axis  = h3axis
    LET h4wLft  = w1Lft
    LET h4wRgt  = h4wLft+slidewidth

    LET h4fLft  = 0
    LET h4fRgt  = 2

    LET h4STik  = 0.1
    LET h4LTik  = 0.5
    LET h4Label = 1
    LET h4First = h4fLft
    LET h4Click = 0.1

    CALL h4SliderVariables

    SUB h4Init
        CALL h4DrawSlider(h4name$,omega)
        CALL SwapOmega(h4wLft-17,h4wBas-3,"w","w",litgry)
    END SUB

    DECLARE DEF h4Wndx


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

    DECLARE PUBLIC cb1Lft,cb1Rgt,cb1Bas,cb1Top,cb1Txt$,cb1Clr,cb1State

    DECLARE PUBLIC chkSiz

    LET cb1Lft = w2Lft
    LET cb1Bas = w2Bas + 40
    LET cb1Txt$= "Direction Field"
    LET cb1Clr = litgry

    CALL cb1Variables

    DECLARE DEF cb1Within


    ! --- checkboxes for graphs ---

    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 + 40
    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

    MAT redim c1Switch(1:c1cnt)
    MAT READ c1Switch
    DATA 0,0,1                                   ! initial state

    DECLARE DEF c1Within

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


    ! --- 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

    SUB t1Set
        CALL SetTextFont(1,12,"bold")
        LOCAL r$,i$,z$
        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 t1Init
        CALL t1Set
    END SUB


    ! --- text rectangle 2 ---

    LET t2BasLn1 = w2Top - 70
    LET t2BasLn2 = t2BasLn1 + 20
    LET t2Lft    = w2Lft - 10
    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$,ssClr)
        CALL PlotTextRJ(t2eqx,t2BasLn2,t2Label2$,ssClr)
        CALL PlotTextCJ(t2Lft+4,t2BasLn2-10,".",ssClr)
    END SUB

    SUB t2Set
        CALL SetTextFont(1,12,"bold")
        CALL t2ClearVal
        IF b>0 and ssFlag=1 then
           CALL PlotTextRJ(t2tRgt,t2BasLn1,using$(t2form$,fxp0),ssClr)
           CALL PlotTextRJ(t2tRgt,t2BasLn2,using$(t2form$,dxp0),ssClr)
        END IF
    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 = w2Top - 70
    LET t3BasLn2 = t3BasLn1 + 20
    LET t3Lft    = w2Lft + 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 - rollover t,x values ---

    LET t4BasLn1= w1Bas + 50
    LET t4BasLn2= t4BasLn1 + 20
    LET t4Lft   = w1Lft + 100
    LET t4Rgt   = t4Lft + 100
    LET t4Bas   = t4BasLn2 +  5
    LET t4Top   = t4BasLn1 - 15
    LET t4Clr   = litgry

    LET t4Label1$= "t = "
    LET t4Label2$= "x = "
    CALL StringWidth(t4Label2$,sw)
    LET t4Eqx= t4Lft + sw

    SUB t4Label
        CALL SuperSubScriptRJ(t4Eqx,t4BasLn1,t4Label1$,t4Clr)
        CALL SuperSubScriptRJ(t4Eqx,t4BasLn2,t4Label2$,t4Clr)
    END SUB

    SUB t4SetCoords
        LET t = w1Fncx(mx)
        LET x = w1Fncy(my)
        LET t$= trim$(using$("--%.##",t))
        LET x$= trim$(using$("--%.##",x))
        CALL t4Clear
        CALL PlotTextLJ(t4Eqx,t4BasLn1,t$,t4Clr)
        CALL PlotTextLJ(t4Eqx,t4BasLn2,x$,t4Clr)
    END SUB

    SUB t4Clear
        BOX CLEAR t4Eqx-2,t4Rgt,t4Bas,t4Top
    END SUB

    SUB t4Init
        CALL t4Label
        ! CALL t4SetCoords
    END SUB


    ! --- radio buttons ---

    DECLARE PUBLIC r1cnt,r1stp,r1siz
    DECLARE PUBLIC r1lft,r1rgt,r1bas,r1top,r1Name$,r1NameClr
    DECLARE PUBLIC r1NameList$(),r1ColorList()

    MAT redim r1NameList$(1:4)
    MAT READ r1NameList$
    DATA "10","2","0.1","0.01"

    MAT redim r1ColorList(1:4)
    MAT READ r1ColorList
    DATA 22,22,22,22

    LET r1cnt    = 4
    LET r1Lft    = w1Lft
    LET r1Top    = w1Bas+40
    LET r1Name$  = ""
    LET r1NameClr= 22

    DECLARE DEF r1Within

    CALL r1SetVars


    ! --- default parameters ---

    LET b,oldb        = 1
    LET k,oldk        = 1
    LET omega,oldomega= 1

    LET dx0 ,fx0    = 0
    LET dxp0,fxp0   = 0
    LET p,v         = 0
    LET vFieldFlag  = 0
    LET cb1State    = 1
    LET vFieldSwitch= cb1State

    LET r1num       = 2

    LET critMarks   = 0
    CALL SetSwitchVars
    CALL SetVars


    ! --- Draw the screen ---

    CALL InitScreen

    SUB InitScreen
        BOX CLEAR worklft,workrgt,workbas,worktop

        CALL w1Init
        CALL w2Init

        CALL v1Init
        CALL h1DrawSlider(h1name$,vel)

        CALL t1Init
        CALL t2Init
        CALL t3Init
        CALL t4Init

        CALL c1DrawCheckBoxes
        CALL c1SetCheckBox

        CALL h2DrawSlider(h2name$,b)

        CALL r1DrawCheckBoxes
        CALL r1SetCheckBox(r1num)

        IF critMarks=1 then
           CALL MarkCritDamp(k)
           CALL MarkNatFreq(k)
        END IF
        CALL h3DrawSlider(h3name$,k)
        CALL h4Init

        IF vFieldFlag=1 then
           CALL cb1Init
           CALL VectorField
        END IF

        CALL w1SliderGraphs
    END SUB

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

    DO
       LET w1Clear= 0
       CALL w1KeepTempLayer
       DO
          GET MOUSE: mx,my,ms
          IF w1wWithin(mx,my)=true then
             IF oldmx<>mx or oldmy<>my then
                CALL w1RollOver
             END IF
          ELSE IF w1Clear=1 then
             CALL w1RollOverClear
          END IF
       LOOP until ms=2                           ! mouse down?

       ! mouse down event

       IF w1wWithin(mx,my)=true then
          DO
             GET MOUSE: mx,my,ms
             IF w1wWithin(mx,my)=true then
                IF oldmx<>mx or oldmy<>my then
                   CALL w1RollOver
                END IF
             ELSE IF w1Clear=1 then
                CALL w1RollOverClear
             END IF
          LOOP until ms=3
          IF w1Clear=1 then CALL w1RollOverClear

       ELSE IF w2wWithin(mx,my)=true then
          LET oldmx,oldmy= -999
          CALL w1xpGraph
          DO
             GET MOUSE: mx,my,ms
             IF mx<>oldmx or my<>oldmy then
                CALL w2wClamp(mx,my)
                CALL w2PixelsToMath(mx,my,dx0,fx0)
                CALL v1Mark(fx0)
                CALL h1Mark(dx0)
                CALL w1PlaneGraphs
                CALL CopyCoords(mx,my,oldmx,oldmy)
                CALL CopyCoords(fx0,dx0,oldfx0,olddx0)
             END IF
          LOOP until ms=3

          ! sliders

       ELSE IF v1Within(mx,my)=true then         ! b
          IF mx<v1axis then
             CALL v1MouseClick
          ELSE
             CALL v1MouseDrag
          END IF

       ELSE IF h1Within(mx,my)=true then         ! b
          IF my<h1axis then
             CALL h1MouseClick
          ELSE
             CALL h1MouseDrag
          END IF

       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 h4Within(mx,my)=true then         ! w
          IF my<h4axis then
             CALL h4MouseClick
          ELSE
             CALL h4MouseDrag
          END IF

          ! buttons

       ELSE IF c1Within(mx,my)=true then         ! graphs
          CALL c1MouseUp(mx,my,ms, c1Pntr)
          LET c1Pntr = clamp(c1Pntr,1,3)
          CALL SetSwitchVars
          IF ssFlag=0 then
             CALL t2ClearVal
             !             CALL ShowGridLayers
             !             CALL KeepGraphLayers
             !          ELSE
             !             CALL w2IVpPoint

          END IF
          IF trFlag=0 and slFlag=0 then CALL t3ClearVal
          CALL w1SliderGraphs

       ELSE IF r1Within(mx,my)=1 then            ! window bounds
          CALL MouseUp(mx,my,ms)
          LET r1num= int((my-r1Top)/r1Stp) + 1
          LET r1num= clamp(r1num,1,4)
          CALL r1SetCheckBox(r1num)
          CALL w1Init
          CALL w1SliderGraphs

       ELSE IF vFieldFlag=true and cb1Within(mx,my)=true then   ! check box
          CALL cb1MouseUp(ms)
          LET vFieldSwitch= cb1State
          CALL VectorField

       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 w1RollOver                               ! cross hairs and cursor values
        CALL w1ShowTempLayer
        SET COLOR litmid
        PLOT mx,w1Top+1; mx,w1Bas-1
        PLOT w1Lft+1,my; w1Rgt-1,my
        CALL t4SetCoords

        LET w1Clear= 1
        CALL CopyCoords(mx,my,oldmx,oldmy)
    END SUB

    SUB w1RollOverClear
        CALL w1ShowTempLayer
        LET w1Clear= 0
        CALL t4Clear
    END SUB

    ! --- Mouse Event Methods ---

    SUB ShowGridLayers
        CALL w1ShowGridLayer
        CALL w2ShowGridLayer
    END SUB
    SUB KeepGraphLayers
        CALL w1KeepGraphLayer
        CALL w2KeepGraphLayer
    END SUB
    SUB ShowGraphLayers
        CALL w1ShowGraphLayer
        CALL w2ShowGraphLayer
    END SUB

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

        CALL SetTextFont(1,12,"bold")
        CALL PlotTextCJ(wcd,cdTop+16,"b",litgry)
        CALL SetTextFont(1,9,"bold")
        CALL PlotTextLJ(wcd+5,cdTop+19,"c",litgry)
    END SUB

    SUB MarkNatFreq(k)
        CALL SetNatFreq(k, natFreq)
        LET wnf  = h4Wndx(natFreq)
        LET nfTop= h4wBas+2
        CALL ValueMarker(h4wLft,h4wRgt,wnf,nftop,litgry)
        CALL SetTextFont(1,12,"bold")
        CALL SwapOmegaN(wnf-5,nfTop+16,"w","w",litgry)
    END SUB

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

    ! ----- v1 slider event -----

    SUB v1MouseClick
        CALL v1GetClickVal(ms,v1Click,fx0)
        CALL v1Action
    END SUB
    SUB v1MouseDrag
        CALL w1xpGraph
        DO
           CALL v1GetDragVal(ms,v1Places,fx0)
           CALL v1Action
        LOOP until ms=3
    END SUB

    SUB v1Action
        IF fx0<>oldfx0 then
           CALL v1Mark(fx0)
           CALL w1PlaneGraphs
           LET oldfx0= fx0
        END IF
    END SUB

    ! ----- h1 slider event -----

    SUB h1MouseClick
        CALL h1GetClickVal(ms,h1Click,dx0)
        CALL h1Action
    END SUB
    SUB h1MouseDrag
        CALL w1xpGraph
        DO
           CALL h1GetDragVal(ms,h1Places,dx0)
           CALL h1Action
        LOOP until ms=3
    END SUB

    SUB h1Action
        IF dx0<>olddx0 then
           CALL h1Mark(dx0)
           CALL w1PlaneGraphs
           LET olddx0= dx0
        END IF
    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

    ! ----- h4 slider event -----

    SUB h4MouseClick
        CALL h4GetClickVal(ms,h4Click,omega)
        CALL h4Action
    END SUB
    SUB h4MouseDrag
        DO
           CALL h4GetDragVal(ms,h4Places,omega)
           CALL h4Action
        LOOP until ms=3
    END SUB

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


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

    SUB InitGraphs
        CALL SetVars
        CALL InitialVals(b,k,omega,D, fxp0,dxp0)
        CALL InitialDifs
    END SUB

    SUB w1xpGraph
        CALL SetVars
        LET pntr= 0
        FOR wx= w1Lft to w1Rgt
            LET t   = w1Fncx(wx)
            LET x1  = xp(t)
            LET ssList(pntr)= x1
            LET pntr= pntr+1
        NEXT wx

        CALL ShowGridLayers
        IF ssFlag=1 and D<>0 then
           CALL w1SteadyStateGraph
        ELSE
           CALL w1ShowGridLayer
        END IF
        CALL KeepGraphLayers
    END SUB

    SUB w1SliderGraphs
        CALL InitGraphs
        CALL ShowGridLayers
        ! CALL w2IVpPoint
        CALL w2IVPoint

        IF k=0 and omega=0 then
           IF b<>0 then
              IF slFlag=1 then CALL SpecialCase2
           ELSE
              IF slFlag=1 then CALL SpecialCase3
           END IF
        ELSE IF b=0 and k=omega^2 then
           IF slFlag=1 then CALL SpecialCase
        ELSE
           CALL MakeSliderLists
           IF ssFlag=1 then
              CALL w1SteadyStateGraph
           END IF
           IF trFlag=1 then CALL w1TransientGraph
           IF slFlag=1 then CALL w1SolutionGraph
        END IF
        CALL w1KeepGraphLayer
    END SUB

    SUB w1PlaneGraphs
        CALL InitGraphs
        !CALL w2IVpPoint

        IF k=0 and omega=0 then
           CALL ShowGridLayers
           CALL w2IVPoint
           IF b<>0 then
              IF slFlag=1 then CALL SpecialCase2
           ELSE
              IF slFlag=1 then CALL SpecialCase3
           END IF
        ELSE IF b=0 and k=omega^2 then
           CALL ShowGridLayers
           CALL w2IVPoint
           IF slFlag=1 then CALL SpecialCase
        ELSE
           CALL MakePlaneLists
           CALL ShowGraphLayers
           CALL w2IVPoint
           IF trFlag=1 then CALL w1TransientGraph
           IF slFlag=1 then CALL w1SolutionGraph
        END IF
    END SUB

    SUB w2IVpPoint                               ! steady state initial values
        CALL t2Set
        CALL w2MathToPixels(dxp0,fxp0,wx,wy)
        IF w2wWithin(wx,wy)=true and ssFlag=1 then
           CALL PlotDiamondClr(wx,wy,ssClr)
        END IF
    END SUB

    SUB w2IVPoint                                ! initial values
        CALL t3Set
        CALL w2MathToPixels(dx0,fx0,wx,wy)
        IF w2wWithin(wx,wy)=true then
           CALL PlotDiamondClr(wx,wy,yellow)
        END IF
    END SUB

    SUB MakeSliderLists
        LET pntr= 0
        FOR wx= w1Lft to w1Rgt
            LET t   = w1Fncx(wx)
            LET x1  = xp(t)
            LET x2  = tr(t)
            LET x3  = x1+x2
            LET ssList(pntr)= x1
            LET trList(pntr)= x2
            LET slList(pntr)= x3
            LET pntr= pntr+1
        NEXT wx
    END SUB

    SUB MakePlaneLists
        LET pntr= 0
        FOR wx= w1Lft to w1Rgt
            LET t   = w1Fncx(wx)
            LET x2  = tr(t)
            LET x3  = xp(t)+x2
            LET trList(pntr)= x2
            LET slList(pntr)= x3
            LET pntr= pntr+1
        NEXT wx
    END SUB

    SUB w1SteadyStateGraph
        CALL w2IVpPoint

        SET COLOR ssClr
        LET oldwy= w1Wndy(xp(0))
        LET oldwx= w1Lft
        LET n= 0

        FOR wx= w1Lft to w1Rgt
            LET x = ssList(n)
            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)
            LET n= n+1
        NEXT wx
        PLOT
    END SUB

    SUB SpecialCase
        IF omega>0 then
           SET COLOR yellow
           LET oldwy= w1Wndy(fx0)
           LET oldwx= w1Lft
           LET n= 0

           FOR wx= w1Lft to w1Rgt
               LET t = w1Fncx(wx)
               LET x = sc1(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)
               LET n= n+1
           NEXT wx
           PLOT
        END IF
    END SUB
    SUB SpecialCase2
        SET COLOR yellow
        LET oldwy= w1Wndy(fx0)
        LET oldwx= w1Lft
        LET n= 0

        FOR wx= w1Lft to w1Rgt
            LET t = w1Fncx(wx)
            LET x = sc2(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)
            LET n= n+1
        NEXT wx
        PLOT
    END SUB
    SUB SpecialCase3
        SET COLOR yellow
        LET oldwy= w1Wndy(fx0)
        LET oldwx= w1Lft
        LET n= 0

        FOR wx= w1Lft to w1Rgt
            LET t = w1Fncx(wx)
            LET x = sc3(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)
            LET n= n+1
        NEXT wx
        PLOT
    END SUB

    SUB w1SolutionGraph
        SET COLOR yellow
        LET oldwy= w1Wndy(fx(0))
        LET oldwx= w1Lft
        LET n= 0
        FOR wx= w1Lft to w1Rgt
            LET x = slList(n)
            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)
            LET n= n+1
        NEXT wx
        PLOT
    END SUB

    SUB w1TransientGraph
        SET COLOR trclr
        LET oldwy= w1Wndy(tr(0))
        LET oldwx= w1Lft
        LET n= 0
        FOR wx= w1Lft to w1Rgt
            LET x = trList(n)
            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)
            LET n= n+1
        NEXT wx
        PLOT
    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= w2Wndx(vel)                  ! phase plane
           IF w2wWithin(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

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

    ! --- vector field routines ---

    SUB VectorField
        CALL w2ShowGridLayer
        IF f=0 and vFieldFlag=1 and vFieldSwitch=1 then
           SET COLOR fldclr
           LET stp = 20
           LET hstp= 10
           FOR wy= w2Top+hstp to w2Bas-hstp step stp
               FOR wx= w2Lft+hstp to w2Rgt-hstp step stp
                   CALL w2PixelsToMath(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 w2KeepGraphLayer
    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 forced damped vibes code ------------------