!! File: PhaseLines !! 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$ = "Phase Lines" SUB ThisProgram CALL PhaseLines 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= 7 ! ELSE IF size=12 then ! numeric output - axis labels ! LET CharWidth= 8 ! ELSE IF size=14 then ! rare ! LET CharWidth= 9 ! 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 Unix header -------- ! *** SUB PhaseLines DECLARE PUBLIC backclr,black,drkgry,drkmid,litmid,litgry,white DECLARE PUBLIC red,yellow,green,cyan,blue,magenta,colorscheme DECLARE PUBLIC planeclr,axisclr,axislabelclr,slideclr DECLARE PUBLIC workLft,workRgt,workBas,workTop,workMidx DECLARE PUBLIC true,false DECLARE PUBLIC SLUMode DECLARE DEF InfoWithin, QuitWithin ! ---- Utility Functions ---- DECLARE DEF badarg,roundn,clamp,modulus,rotation ! --- help screen array --- DIM info$(1:16) MAT READ info$ DATA "Phase Lines" DATA "" DATA "This program lets you view the phase line, the solution graphs, the graph of the differential equation, and the bifurcation diagram for families of several first-order differential equations. Stable equilibria are shown in green, unstable equilibria are shown in red, and semistable equilibria are shown in cyan." DATA "" DATA "Click the t,x plane to graph the solution curve for initial values t and x." DATA "" DATA "Click the [v] equation menu button to display the equations. Click on an equation to select it." DATA "" DATA "Click or drag the a slider handle change the value of the parameter and redraw the direction field." DATA "" DATA "Click a [Clear] button to remove solution graphs from the time series plane or the phase line." DATA "Click the [Phase Line] button to see solutions on the phase line." DATA "Click the [Bifurcation Diagram] button display the bifurcation plot." DATA "" DATA "" DATA "Adapted from Differential Equations, by Blanchard, Devaney, and Hall" ! ---------- Utility functions ---------- DECLARE DEF clamp,roundn,e DEF notanum= 987656789 ! --- functions, parameters, equations, constants --- IF M68kFlag=1 then LET deltat= 1/16 ELSE LET deltat= 1/8 END IF ! 1. Logistic ! 2. Saddle Node Bifurcation ! 3. Subcritical Pitchfork Bifurcation ! 4. Transcritical Bifurcation ! 5. Transcritical Bifurcation ! 6. Imperfect Pitchfork ! y' = 1 - a*y -2 < a < 2 ! ! y' = (1/4) - a*y + y^2 -1.5 < a < 1.5 DEF dxdt(eq,a,y) SELECT CASE eq CASE 1 ! logistic LET dy= y*(1-y) - a CASE 2 ! saddle node LET dy= y*y + a CASE 3 ! subcritical pitchfork LET dy= a*y + y^3 CASE 4 ! transcritical LET dy= a*y - y*y CASE 5 ! LET dy= 1 - a*y CASE 6 ! LET dy= .25 - a*y + y^2 ! imperfect pitchfork ! CASE 5 ! ! LET dy= r + a*y - y*y ! CASE 6 ! ! LET dy= r + a*y - y^3 ! imperfect pitchfork END SELECT LET dxdt= dy END DEF DEF df(x) = a - 3*x*x ! --- menu --- ! DIM BUMenu$(0:5) ! DATA "y' = y(1-y) - a" ! DATA "y' = y^[2] + a" ! DATA "y' = ay + y^[3]" ! DATA "y' = ay - y^[2]" ! DATA "y' = r + ay - y^[2]" ! DATA "y' = r + ay - y^[3]" ! MAT READ BUMenu$ ! DIM SLUMenu$(0:5) ! DATA "x' = x(1-x) - a" ! DATA "x' = x^[2] + a" ! DATA "x' = ax + x^[3]" ! DATA "x' = ax - x^[2]" ! DATA "x' = r + ax - x^[2]" ! DATA "x' = r + ax - x^[3]" ! MAT READ SLUMenu$ ! --- Graphing plane parameters and methods --- ! ------------- w1 plane data - t,x -------------- 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 w1xPiAxis= 0 LET w1yPiAxis= 0 LET w1xMult = 1 LET w1yMult = 1 LET whgt = 280 LET w1Lft = workLft+55 ! pixel bounds LET w1Rgt = w1Lft+450 ! for 6 LET w1Top = workTop+30 LET w1Bas = w1Top+whgt LET w1fLft= 0 ! function bounds LET w1fRgt= 8 LET w1fTop= 2 LET w1fBas= -2 LET w1xAx$= "t" ! axis labels LET w1yAx$= "y" LET w1xGridstep= 0 ! horizontal grid intervals LET w1yGridstep= 0 ! vertical grid intervals LET w1xSTik = 0.5 ! horizontal axis Tik marks LET w1xLTik = 1 LET w1xLabel= 1 LET w1xFirst= w1fLft LET w1ySTik = 0.5 ! vertical axis Tik marks LET w1yLTik = 1 LET w1yLabel= 1 LET w1yFirst= w1fBas ! --- w1 Plane methods --- DECLARE DEF w1fncx,w1fncy,w1wndx,w1wndy ! window/function transforms DECLARE DEF w1Within, w1wWithin, w1fWithin CALL w1Variables SUB w1Init 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 CALL PlotTextCJ(w1x0,w1Top-10,w1yAx$,axislabelclr) CALL w1Buttons CALL w1KeepGridLayer END SUB ! ------------------------------------------ ! --- w2 plane data: phase line --- 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 w2xPiAxis= 0 LET w2yPiAxis= 0 LET w2xMult = 1 LET w2yMult = 1 LET w2Lft = w1Rgt + 55 ! pixel bounds LET w2Rgt = w2Lft + 18 LET w2Top = w1Top LET w2Bas = w1Bas LET w2fLft= 0 ! function bounds * pi LET w2fRgt= 1 LET w2fTop= w1fTop LET w2fBas= w1fBas LET w2xAx$= "" ! axis labels LET w2yAx$= "y" LET w2xGridstep= 0 ! horizontal grid intervals LET w2yGridstep= 0 ! vertical grid intervals LET w2xSTik = 0 ! horizontal axis Tik marks LET w2xLTik = 0 LET w2xLabel= 0 LET w2xFirst= w2fLft LET w2ySTik = 0.5 ! vertical axis Tik marks LET w2yLTik = 1 LET w2yLabel= 1 LET w2yFirst= w2fBas ! --- w2 plane methods --- DECLARE DEF w2fncx, w2fncy, w2wndx, w2wndy ! window/function transforms DECLARE DEF w2Within, w2wWithin CALL w2Variables LET w2pntLft= w2Lft+10 LET w2pntRgt= w2Rgt- 7 SUB w2Init CALL w2DrawPlane(0,1,0) ! x, y, zeroaxes CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(w2Rgt+8,w2y0+3,w2xAx$,axislabelclr) ! axis labels CALL PlotTextCJ(w2x0,w2Top-10,w2yAx$,axislabelclr) CALL w2Buttons CALL w2KeepGridLayer END SUB ! ------------------------------- ! --- w3 plane data: DE graph --- 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 w3xPiAxis= 0 LET w3yPiAxis= 0 LET w3xMult = 1 LET w3yMult = 1 LET w3Lft = w1lft ! pixel bounds LET w3Rgt = w3Lft+100 LET w3Top = w1Bas+ 90 LET w3Bas = w3Top+100 LET w3fLft= -2 ! function bounds LET w3fRgt= 2 LET w3fTop= 2 LET w3fBas= -2 LET w3xAx$= "y" ! axis labels LET w3yAx$= "y'" LET w3xGridstep= 0 ! grid line intervals LET w3yGridstep= 0 LET w3xSTik = .5 ! axis Tik marks LET w3xLTik = 1 LET w3xLabel= 1 LET w3xFirst= w3fLft LET w3ySTik = .5 !.1 LET w3yLTik = 1 LET w3yLabel= 1 LET w3yFirst= w3fBas ! --- w3 plane methods --- DECLARE DEF w3fncx,w3fncy,w3wndx,w3wndy ! window/function transforms CALL w3Variables SUB w3Init CALL w3DrawPlane(1,1,1) ! grid, axes, zeroaxes CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(w3Rgt+8,w3y0+3,w3xAx$,axislabelclr) ! axis labels CALL PlotTextCJ(w3x0+2,w3Top-11,w3yAx$,axislabelclr) CALL w3KeepGridLayer END SUB ! --- w4 plane data: bifurcation plane --- DECLARE PUBLIC w4Lft,w4Rgt,w4Bas,w4Top,w4Midx,w4Midy DECLARE PUBLIC w4fLft,w4fRgt,w4fBas,w4fTop,w4x0,w4y0 DECLARE PUBLIC w4xFirst, w4xSTik, w4xLTik, w4xLabel, w4xGridstep DECLARE PUBLIC w4yFirst, w4ySTik, w4yLTik, w4yLabel, w4yGridstep DECLARE PUBLIC w4wWid,w4wHgt,w4fWid,w4fHgt DECLARE PUBLIC w4fxRatio,w4fyRatio,w4wxRatio,w4wyRatio,w4Aspect DECLARE PUBLIC w4xPiFlag, w4xMult, w4yPiFlag, w4yMult LET w4Flag = 1 LET w4xPiAxis= 0 LET w4yPiAxis= 0 LET w4xMult = 1 LET w4yMult = 1 LET w4Lft = workRgt-145 ! pixel bounds LET w4Rgt = w4Lft + 84 LET w4Top = w1Top LET w4Bas = w1Bas LET w4fLft= -2 ! function bounds LET w4fRgt= 2 LET w4fTop= 2 LET w4fBas= -2 LET w4xAx$= "a" ! axis labels LET w4yAx$= "y" LET w4xGridstep= 0 ! grid line intervals LET w4yGridstep= 0 LET w4xSTik = 0.5 ! axis Tik marks LET w4xLTik = 1 LET w4xLabel= 1 LET w4xFirst= w4fLft LET w4ySTik = 0.5 LET w4yLTik = 1 LET w4yLabel= 1 LET w4yFirst= w4fBas ! --- w4 plane methods --- DECLARE DEF w4fncx,w4fncy,w4wndx,w4wndy ! window/function transforms CALL w4Variables SUB w4Init CALL w4DrawPlane(1,1,1) ! grid, axes, zeroaxes CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(w4Rgt+8,w4y0+3,w4xAx$,axislabelclr) ! axis labels CALL PlotTextCJ(w4x0,w4Top-10,w4yAx$,axislabelclr) CALL w4KeepGridLayer END SUB ! ----- menu ----- DECLARE PUBLIC m1Lft, m1Rgt, m1Bas, m1Top, m1Equation DECLARE PUBLIC m1Prefix$, m1Menu1$(),m1tClr DECLARE DEF m1Within MAT redim m1Menu1$(1:6) MAT READ m1Menu1$ DATA "y(1-y) - a" DATA "y^[2] + a" DATA "ay + y^[3]" DATA "ay - y^[2]" DATA "1 - ay" DATA ".25 - ay + y^[2]" !Here they are: ! !y' = 1 - a*y -2 < a < 2 ! !y' = (1/4) - a*y + y^2 -1.5 < a < 1.5 ! DATA "r + ay - y^[2]" ! DATA "r + ay - y^[3]" LET m1Prefix$ = "y' = " LET m1Lft = worklft + 200 LET m1Top = w3Midy - 8 ! workbas - 75 LET m1tClr = yellow LET m1Equation= 1 ! ---------- Slider parameters and methods ---------- ! --------------- horizontal sliders -------------- ! --- h1 slider - a --- DECLARE PUBLIC h1axis,h1wLft,h1wRgt,h1wBas,h1wTop,h1wWid,h1sLft,h1sRgt DECLARE PUBLIC h1fLft,h1fRgt,h1fWid,h1First,h1STik,h1LTik,h1Label DECLARE PUBLIC h1fratio,h1wratio,h1name$,h1form$,h1clr,h1slot$ DECLARE PUBLIC h1PiAxis, h1Mult LET h1Flag = 0 LET h1PiAxis= 0 LET h1Mult = 1 LET h1clr = slideclr LET h1name$ = "a" LET h1form$ = "-%.##" LET h1Places= 2 LET h1Click = 0.5 LET h1axis = m1Top LET h1wLft = w4Lft LET h1wRgt = w4Rgt LET h1fLft = -2 LET h1fRgt = 2 LET h1STik = 0.5 ! short tick marks LET h1LTik = 1 ! long tick marks LET h1Label = 1 ! labels LET h1First = h1fLft ! first tick mark ! --- h1 methods --- CALL h1SliderVariables DECLARE DEF h1Fncx,h1Within !,h1wWithin,h1fWithin ! window/function transforms SUB h1Init CALL h1SliderBounds(m1Equation) CALL h1DrawSlider(h1name$,a) END SUB SUB h1SliderBounds(eq) SELECT CASE eq CASE 1 LET h1fLft = 0 LET h1fRgt = 0.5 LET h1STik = 0.1 ! short tick marks LET h1LTik = 0.1 ! long tick marks LET h1Label= 0.5 ! labels ! CASE 6 ! LET h1fLft = -1.5 ! LET h1fRgt = 1.5 ! ! LET h1STik = 0.5 ! short tick marks ! LET h1LTik = 1.5 ! long tick marks ! LET h1Label= 1.5 ! labels ! CASE 6 ! LET h1fLft = -3 ! LET h1fRgt = 3 ! ! LET h1STik = 0.5 ! short tick marks ! LET h1LTik = 1 ! long tick marks ! LET h1Label= 3 ! labels CASE else LET h1fLft = -2 LET h1fRgt = 2 LET h1STik = 0.5 ! short tick marks LET h1LTik = 1 ! long tick marks LET h1Label= 1 ! labels END SELECT LET w4xSTik = h1Stik ! axis Tik marks LET w4xLTik = h1Ltik LET w4xLabel= h1Label LET w4fLft = h1fLft ! function bounds LET w4fRgt = h1fRgt LET h1First = h1fLft ! first tick mark LET w4xFirst= w4fLft ! first tick mark LET h1Click = h1Stik CALL h1SliderVariables CALL w4Variables END SUB ! --- h2 slider - r --- DECLARE PUBLIC h2axis,h2wLft,h2wRgt,h2wBas,h2wTop,h2wWid,h2sLft,h2sRgt DECLARE PUBLIC h2fLft,h2fRgt,h2fWid,h2First,h2STik,h2LTik,h2Label DECLARE PUBLIC h2fratio,h2wratio,h2name$,h2form$,h2clr,h2slot$ DECLARE PUBLIC h2PiAxis, h2Mult LET h2Flag = 0 LET h2PiAxis= 0 LET h2Mult = 1 LET h2clr = slideclr LET h2name$ = "r" LET h2form$ = "-%.##" LET h2Places= 2 LET h2Click = h2Stik LET h2axis = h1axis + 40 LET h2wLft = w4Lft LET h2wRgt = w4Rgt LET h2fLft = 0 LET h2fRgt = 1 LET h2STik = 0 ! short tick marks LET h2LTik = .5 ! long tick marks LET h2Label = 1 ! labels LET h2First = h2fLft ! first tick mark ! --- h2 methods --- CALL h2SliderVariables DECLARE DEF h2Fncx,h2Within SUB h2SliderBounds(eq) SELECT CASE eq CASE 5 LET h2STik = 0.25 ! short tick marks LET h2LTik = 0.5 ! long tick marks LET h2Label= 1 ! labels LET h2fLft = -1 LET h2fRgt = 1 CASE 6 LET h2STik = 0.5 ! short tick marks LET h2LTik = 1 ! long tick marks LET h2Label= 1 ! labels LET h2fLft = -1 LET h2fRgt = 1 CASE else END SELECT LET h2First = h2fLft ! first tick mark LET h2Click = h2Stik CALL h2SliderVariables END SUB ! ---------- Orbit List ---------- ! not used at MIT ! DIM table(0:16) ! ! LET tblSpc = 15 ! LET tblCnt = 15 ! LET tblLft = w1rgt-155 ! LET tblRgt = tblLft + 155 ! LET tblTop = w1bas + 45 ! LET tblBas = tblTop + tblcnt*tblspc + 4 ! LET tblCol1 = tblLft + 45 ! LET tblCol2 = tblCol1 + 55 ! LET tblCol3 = tblRgt ! LET tblFirst= tbltop+tblSpc-2 ! ! SUB TableClear ! BOX CLEAR tblLft,tblRgt,tblBas,tblTop ! END SUB ! ! SUB DrawTable ! SET COLOR planeclr ! BOX AREA tblLft,tblRgt,tblBas,tblTop ! SET COLOR rimclr ! BOX LINES tblLft,tblRgt,tblBas,tblTop ! PLOT tblCol1,tblTop; tblCol1,tblBas ! PLOT tblCol2,tblTop; tblCol2,tblBas ! CALL SetTextFont(1,12,"bold") ! LET bas= tblTop-3 ! CALL PlotTextCJ((tblLft +tblCol1)/2,bas,"t",trajclr) ! CALL PlotTextCJ((tblCol1+tblCol2)/2,bas,yax$,trajclr) ! CALL PlotTextCJ((tblCol2+ tblRgt)/2,bas,"d"&yax$&"/dt",trajclr) ! END SUB ! ! SUB TableRefresh ! LET tb1Pntr= 0 ! CALL DrawTable ! END SUB ! ! SUB StoreValues(t,x) ! LET tblPntr= round(2*t) ! LET table(tblPntr)= x ! LET t$= using$("%.#",t) ! LET x$= using$("-%.##",x) ! LET ty= tblFirst + 2*t*tblSpc ! CALL PlotTextRJ(tblCol1-10,ty,t$,trajclr) ! CALL PlotTextRJ(tblCol2-10,ty,x$,trajclr) ! LET dx$= using$("-%.##",dxdt(m1Equation,a,x)) ! CALL PlotTextRJ(tblCol3-10,ty,dx$,trajclr) ! END SUB DECLARE PUBLIC tb1Lft,tb1Rgt,tb1Bas,tb1Top,tb1Form$ DECLARE PUBLIC tb1RowCnt,tb1RowStp,tb1ColCnt,tb1ColStp,tb1Pntr DECLARE PUBLIC tb1(,),tb1ColLabels$(),tb1ColColors() LET tb1RowStp= 15 LET tb1RowCnt= 15 ! w1fWid*2 + 1 LET tb1ColStp= 60 LET tb1ColCnt= 2 LET tb1Lft = w1rgt - 121 ! 155 LET tb1Top = w1bas + 53 LET tb1Form$ = "---%.##" CALL tb1Variables DATA "t","y" MAT READ tb1ColLabels$ DATA 8,8 MAT READ tb1ColColors SUB table1Init CALL tb1Init BOX KEEP tb1Lft,tb1Rgt,tb1Bas,tb1Top in tb1$ END SUB SUB tb1Refresh BOX SHOW tb1$ at tb1Lft,tb1Bas LET tb1Pntr= 0 END SUB SUB tb1Store(x,y) LET tb1Pntr= tb1Pntr+1 IF tb1Pntr<=tb1RowCnt then LET tb1(tb1Pntr,1)= x LET tb1(tb1Pntr,2)= y CALL tb1Values(tb1,tb1Pntr) END IF END SUB ! ---------- Text Output Rects ---------- ! ---------- text rectangle 2 --- LET t2BasLn1= w1Bas + 42 LET t2BasLn2= t2BasLn1 + 20 LET t2BasLn3= t2BasLn2 + 20 LET t2Top = t2BasLn1 - 15 LET t2Bas = t2BasLn3 + 5 LET t2Lft = w1Lft + 181 CALL StringWidth("yo = ",sw) LET t2Eqx = t2Lft + sw LET t2Rgt = t2Eqx + 100 SUB t2Label CALL SuperSubScriptRJ(t2Eqx-2,t2BasLn1,w1xax$ & "_[0] = ",trajclr) CALL SuperSubScriptRJ(t2Eqx-2,t2BasLn2,w1yax$ & "_[0] = ",trajclr) CALL PlotTextRJ(t2Eqx,t2BasLn3,w1yax$ & "' = ",trajclr) END SUB SUB t2Set(plane) LOCAL t$,x$,dx$ LET t$ = trim$(using$("-%.##",t)) LET x$ = trim$(using$("-%.##",x)) LET dx$= trim$(using$("-%.##",dx)) CALL t2Clear IF plane=1 then CALL PlotTextLJ(t2Eqx,t2BasLn1, t$,trajclr) CALL PlotTextLJ(t2Eqx,t2BasLn2, x$,trajclr) CALL PlotTextLJ(t2Eqx,t2BasLn3,dx$,trajclr) END SUB SUB t2Clear BOX CLEAR t2Eqx-2,t2Rgt,t2Bas,t2Top END SUB SUB t2Init CALL t2Label ! CALL t2Set END SUB ! --- checkbox - Phase Line --- DECLARE PUBLIC cb1Lft,cb1Rgt,cb1Bas,cb1Top,cb1Txt$,cb1Clr,cb1State DECLARE DEF cb1Within LET cb1Lft = w1Rgt-120 LET cb1Bas = m1Top+ 16 LET cb1Txt$= "Phase Line" LET cb1Clr = litgry CALL cb1Variables ! --- checkbox - Bifurcation Plane --- DECLARE PUBLIC cb2Lft,cb2Rgt,cb2Bas,cb2Top,cb2Txt$,cb2Clr,cb2State DECLARE DEF cb2Within LET cb2Lft = cb1Lft LET cb2Bas = cb1Bas + 20 LET cb2Txt$= "Bifurcation Diagram" LET cb2Clr = litgry CALL cb2Variables ! ----- buttons ----- LET btnhgt= 18 LET clft = w1Lft+1 LET crgt = clft+50 LET ctop = w1Bas+30 LET cbas = ctop+btnhgt LET c2lft = w2Lft-15 LET c2rgt = w2Rgt+15 LET c2top = ctop LET c2bas = cbas SUB w1Buttons CALL SetTextFont(1,12,"bold") CALL DrawButton(clft,crgt,cbas,ctop,5,"Clear") END SUB SUB w2Buttons CALL SetTextFont(1,12,"bold") CALL DrawButton(c2lft,c2rgt,c2bas,c2top,5,"Clear") END SUB ! --- default parameters --- LET trajclr = yellow LET phsclr = trajclr LET fieldclr= drkmid LET w2Flag = 0 LET w3Flag = 1 LET w4Flag = 0 LET t1Flag = 0 LET a,olda = 0 LET r = 0.5 LET h1flag = 1 LET m1equation= 1 LET m1State = 0 LET w2Flag = 0 LET w4Flag = 0 ! --- draw the screen --- ! CALL InitPalette CALL InitScreen CALL SetTimer CALL m1InitMenu1 SUB InitScreen BOX CLEAR worklft,workrgt,workbas,worktop CALL cb1Init CALL cb2Init LET w1Bas= w1Top+whgt CALL w1Variables CALL w1Buttons CALL w1Init LET slopeflag= 1 CALL SlopeField LET w2Bas= w1Bas CALL w2Variables LET w4Bas= w1Bas CALL w4Variables CALL h1Init CALL m1ResetMenu(m1State) CALL t2Init ! ---- CALL EqLines CALL ResetPlanes END SUB SUB ResetPlanes ! phase line IF w2Flag=1 then CALL w2Init CALL w2EqLines CALL PhaseLine(2) END IF ! DEGraph( IF w3Flag=1 then CALL w3Init CALL DEGraph(1) END IF ! parameter plane IF w4Flag=1 then CALL w4Init CALL w4DrawGraph CALL PhaseLine(4) ELSE CALL GetImperfectNose END IF END SUB ! ----------------- Event manager ----------------- DO DO GET MOUSE: mx,my,ms IF w1wWithin(mx,my)=true then IF mx<>oldmx or my<>oldmy then CALL w1RollOver LET oldmx= mx LET oldmy= my END IF ELSE IF w2wWithin(mx,my)=true and w2Flag=1 then IF my<>oldmy then CALL w2RollOver LET oldmx= mx LET oldmy= my END IF ELSE IF clearflag=1 then CALL CleanUp END IF LOOP until ms=2 IF clearflag=1 then CALL CleanUp ! --- IF w1wWithin(mx,my)=true then LET x= w1Fncy(my) LET t= w1Fncx(mx) CALL t2Set(1) CALL w1MouseClick(mx,my) ELSE IF w2wWithin(mx,my)=true then CALL w2MouseClick(mx,my) ELSE IF h1Within(mx,my)=true then ! h1 slider - a value IF mycLft and mxcTop and myc2Lft and mxc2Top and mybLft and mxbTop and my=w1Lft and wt<=w1Rgt and wx>=w1Top and wx<=w1Bas then PLOT mx,my; wt,wx END IF IF w2Flag=1 then CALL w2ShowGraphLayer IF dx<>0 then SET COLOR trajclr DRAW PhaseArrow with scale(1,sgn(dx)) * shift(w2midx+2,my) END IF END IF CALL t2Set(1) LET clearflag= 1 END SUB PICTURE arrow PLOT -2,-2; 2,0; -2,2 END PICTURE ! ----- SUB w2RollOver LET t = 0 LET x = w1fncy(my) LET dx= dxdt(m1Equation,a,x) CALL w1ShowGraphLayer !IF w2Flag=1 then CALL w2ShowGraphLayer IF dx<>0 then SET COLOR trajclr DRAW PhaseArrow with scale(1,sgn(dx)) * shift(w2midx+2,my) END IF !END IF CALL t2Set(2) LET clearflag= 1 END SUB SUB CleanUp CALL w1ShowGraphLayer CALL t2Clear IF w2Flag=1 then CALL w2ShowGraphLayer LET clearflag= 0 END SUB ! --- plane 1 mouse click --- SUB w1MouseClick(mx,my) CALL MouseUp(mx,my,ms) IF mx>w1lft-3 and mxw1top-3 and myw2lft-3 and mxw2top-3 and myolda then CALL h1Mark(a) CALL UpdateGraphs1 LET olda= a END IF END SUB SUB UpdateGraphs1 IF w2Flag=1 then CALL w2ShowGridLayer CALL w2EqLines CALL PhaseLine(2) END IF IF w3Flag=1 then CALL DEGraph(1) IF w4Flag=1 then CALL w4ShowGraphLayer CALL PhaseLine(4) END IF CALL w1ShowGridLayer IF slopeflag=1 then CALL SlopeField CALL EqLines END SUB ! --- h2 slider - r value --- SUB h2MouseClick(mx,my) CALL h2GetClickVal(ms,h2Click,r) CALL h2Action END SUB SUB h2MouseDrag(mx,my) DO CALL h2GetDragVal(ms,h2Places,r) CALL h2Action LOOP until ms=3 END SUB SUB h2Action IF r<>oldr then CALL h2Mark(r) IF w4Flag=0 then CALL GetImperfectNose ! finds nose when 4 is off END IF CALL UpdateGraphs LET oldr= r END IF END SUB ! --- SUB UpdateGraphs CALL w1ShowGridLayer IF slopeflag=1 then CALL SlopeField CALL EqLines IF w2Flag=1 then CALL w2ShowGridLayer CALL w2EqLines CALL PhaseLine(2) END IF IF w3Flag=1 then CALL DEGraph(1) IF w4Flag=1 then CALL w4DrawGraph CALL PhaseLine(4) END IF END SUB ! --- Graph drawing methods --- SUB RungeKutta4(x) LET dx1= dxdt(m1Equation,a,x) LET x1 = x + 0.5*deltat*dx1 LET dx2= dxdt(m1Equation,a,x1) LET x2 = x + 0.5*deltat*dx2 LET dx3= dxdt(m1Equation,a,x2) LET x3 = x + 1.0*deltat*dx3 LET dx4= dxdt(m1Equation,a,x3) LET dx = (dx1 + 2*dx2 + 2*dx3 + dx4) / 6 LET x = x + deltat*dx END SUB SUB Trajectory(x0) CALL w1ShowGraphLayer CALL Orbit LET deltat= -deltat CALL Orbit LET deltat= -deltat ! CALL SetTextFont(1,12,"bold") CALL w1KeepGraphLayer IF w2Flag=1 then CALL w2KeepGraphLayer END IF END SUB SUB Orbit LET x = x0 LET t = t0 LET oldwy= 999 SET COLOR trajclr CALL w1MathToPixels(t,x,wx,wy) CALL CopyCoords(x,t,oldx,oldt) FOR i= 0 to 100000 ! draw trajectory CALL CopyCoords(wx,wy,oldwx,oldwy) CALL w1MathToPixels(t,x,wx,wy) IF w1wWithin(wx,wy)=true then IF w2Flag=1 then PLOT w2pntLft,wy; w2pntRgt,wy IF oldwx>=w1lft and oldwx<=w1rgt and oldwy>=w1top and oldwy<=w1bas then PLOT wx,wy; oldwx,oldwy ELSE CALL CopyCoords(t,x,px,py) CALL CopyCoords(oldt,oldx,qx,qy) CALL ClipLine(w1flft,w1frgt,w1fbas,w1ftop,px,py,qx,qy,drawflag) IF drawflag=1 then PLOT w1Wndx(px),w1Wndy(py); w1Wndx(qx),w1Wndy(qy) END IF ELSE CALL CopyCoords(t,x,px,py) CALL CopyCoords(oldt,oldx,qx,qy) CALL ClipLine(w1flft,w1frgt,w1fbas,w1ftop,px,py,qx,qy,drawflag) IF drawflag=1 then PLOT w1Wndx(px),w1Wndy(py); w1Wndx(qx),w1Wndy(qy) END IF CALL CopyCoords(t,x,oldt,oldx) CALL RungeKutta4(x) LET t= t + deltat ! GET MOUSE: mx,my,ms IF abs(x)>4 or t<= w1fLft or t>=w1fRgt then EXIT FOR CALL Delay(1/100) NEXT i PLOT END SUB ! --- plane 1 subs: t,x plane --- SUB EqLines CALL EqPoints(m1Equation,a) IF eq0<>notanum then CALL eqColor(eq0State,clr) LET wy= w1Wndy(eq0) IF wy>=w1Top and wy<=w1Bas then CALL PlotLine( w1lft+1,wy, w1rgt-1,wy, clr) END IF END IF IF eq1<>notanum then CALL eqColor(eq1State,clr) LET wy= w1Wndy(eq1) IF wy>=w1Top and wy<=w1Bas then CALL PlotLine( w1lft+1,wy, w1rgt-1,wy, clr) END IF END IF IF eq2<>notanum then CALL eqColor(eq2State,clr) LET wy= w1Wndy(eq2) IF wy>=w1Top and wy<=w1Bas then CALL PlotLine( w1lft+1,wy, w1rgt-1,wy, clr) END IF END IF CALL w1KeepGraphLayer END SUB SUB eqcolor(state,clr) SELECT CASE state CASE 0 SET COLOR red ! resistor LET clr= red CASE .5, -.5 SET COLOR cyan ! half attractor LET clr= cyan CASE 1 SET COLOR green ! attractor LET clr= green CASE else !PRINT state END SELECT END SUB SUB SlopeField ! h=450 v=180 CALL w1ShowGridLayer SET COLOR fieldclr ! vector field LET dx= w1Aspect FOR wy= w1top+15 to w1bas-10 step 15 LET y= w1Fncy(wy) IF round(y,2)<>0 then LET dy = dxdt(m1Equation,a,y) LET ang= angle(dx,-dy) LET sa = 6*sin(ang) LET wy1= wy+sa LET wy2= wy-sa LET ca = 6*cos(ang) FOR wx= w1lft+30 to w1rgt-20 step 30 PLOT wx-ca,wy2; wx+ca,wy1 NEXT wx END IF NEXT wy CALL w1KeepGraphLayer END SUB ! ---- plane 2 subs: phase line ---- SUB w2EqLines CALL EqPoints(m1Equation,a) IF eq0<>notanum then CALL eqColor(eq0State,clr) LET wy= w2wndy(eq0) IF wy>=w2Top and wy<=w2Bas then PLOT w2lft+1,wy; w2rgt-1,wy END IF END IF IF eq1<>notanum then CALL eqColor(eq1State,clr) LET wy= w2wndy(eq1) IF wy>=w2Top and wy<=w2Bas then PLOT w2lft+1,wy; w2rgt-1,wy END IF END IF IF eq2<>notanum then CALL eqColor(eq2State,clr) LET wy= w2wndy(eq2) IF wy>=w2Top and wy<=w2Bas then PLOT w2lft+1,wy; w2rgt-1,wy END IF END IF CALL w2KeepGraphLayer END SUB ! --- plane 3 subs: DE graph --- SUB DEGraph(stp) CALL w3ShowGridLayer SET COLOR yellow ! axisclr LET x = w3Fncx(w3Lft) LET y = dxdt(m1Equation,a,x) FOR wx= w3lft to w3rgt step stp ! graph (x,x') CALL CopyCoords(x,y,oldx,oldy) LET x= w3Fncx(wx) LET y= dxdt(m1Equation,a,x) IF y<=w3ftop and y>=w3fbas and oldy<=w3ftop and oldy>=w3fbas then LET wy= w3Wndy(y) PLOT wx,wy; w3Wndx(oldx),w3Wndy(oldy) ELSE CALL CopyCoords(x,y,px,py) CALL CopyCoords(oldx,oldy,qx,qy) CALL ClipLine(lft,rgt,bas,top,px,py,qx,qy,drawflag) IF drawflag=1 then PLOT w3Wndx(px),w3Wndy(py); w3Wndx(qx),w3Wndy(qy) END IF END IF NEXT wx PLOT CALL EqPoints(m1Equation,a) IF eq0<>notanum then CALL eqColor(eq0State,clr) ! SET CURSOR 1,1 ! PRINT eq0state,clr LET eqx= w3Wndx(eq0) LET eqy= w3Wndy(0) IF eqx>=w3Lft and eqx<=w3rgt then CALL EqPointDE(eqx,eqy,eq0State,clr) END IF END IF IF eq1<>notanum then CALL eqColor(eq1State,clr) LET eqx= w3Wndx(eq1) LET eqy= w3Wndy(0) CALL EqPointDE(eqx,eqy,eq1State,clr) END IF IF eq2<>notanum then CALL eqColor(eq2State,clr) LET eqx= w3Wndx(eq2) LET eqy= w3Wndy(0) CALL EqPointDE(eqx,eqy,eq2State,clr) END IF END SUB ! ------ plane 4 subs: bifurcation plane ------ SUB w4FlowField SET COLOR fieldclr LET hstp= round(w4wwid/4) FOR wk= w4Lft+hstp to w4Rgt-10 step hstp ! flow in bifurcation plane LET k = w4fncx(wk) CALL EqPoints(eq,k) LET wy1 = w4Wndy(eq1) ! lower curve LET wy2 = w4Wndy(eq2) ! upper curve LET dist= 6 LET ystp= 12 FOR wy= w4Top+ystp to w4Bas-10 step ystp IF wywy1+dist then ! DRAW FlowPos4 with shift(wk,wy) PLOT wk-4,wy-2; wk,wy+2; wk+4,wy-2 ELSE IF wy>wy2+dist and wy=0 then ! IF r=0 then ! SET COLOR green ! PLOT w4Lft,w4y0; w4x0,w4y0; w4Rgt,w4Wndy(1.5) ! SET COLOR red ! PLOT w4Lft,w4Wndy(-1.5); w4x0,w4y0; w4Rgt,w4y0 ! ELSE ! SET COLOR red ! FOR wk= w4Lft to w4Rgt step 1/8 ! LET k= w4fncx(wk) ! CALL EqPoints(m1Equation,k) ! ! LET wy= w4wndy(eq1) ! IF wy>= w4Top and wy<= w4Bas then ! PLOT wk,wy; ! ELSE ! PLOT ! END IF ! NEXT wk ! PLOT ! ! SET COLOR green ! FOR wk= w4Lft to w4Rgt step 1/8 ! LET k= w4fncx(wk) ! CALL EqPoints(m1Equation,k) ! ! LET wy= w4wndy(eq2) ! IF wy>= w4Top and wy<= w4Bas then ! PLOT wk,wy; ! ELSE ! PLOT ! END IF ! NEXT wk ! PLOT ! END IF ! END IF ! CASE 6 ! LET x1= 0 ! LET x2= 0 ! IF r=0 then ! SET COLOR green ! PLOT w4Lft,w4y0; w4x0,w4y0 ! SET COLOR red ! PLOT w4x0,w4y0; w4Rgt,w4y0 ! ! FOR wx= w4x0 to w4Rgt ! LET r1 = w4Fncx(wx) ! LET oldx1= x1 ! LET oldx2= x2 ! LET x1 = sqr(r1) ! LET x2 = -x1 ! SET COLOR green ! PLOT wx,w4Wndy(x1); wx-1,w4Wndy(oldx1) ! PLOT wx,w4Wndy(x2); wx-1,w4Wndy(oldx2) ! NEXT wx ! ! ELSE ! r <> 0 ! LET stp= 1 ! FOR wy= w4Top to w4Bas step stp ! LET x= w4Fncy(wy) ! LET oldr2= r1 ! IF x=0 then ! IF r<0 then ! LET r1= w4fLft ! ELSE IF r>0 then ! IF oldr2>w4fLft then ! SET COLOR green ! stable ! PLOT w4Lft,wy; w4Wndx(oldr2),wy-stp ! END IF ! LET r1= w4fRgt ! END IF ! END IF ! ! IF x<>0 then ! LET r1= (x*x*x-r)/x ! IF r>0 then ! IF x>0 then ! SET COLOR green ! stable ! IF r1>w4fLft and r1w4fRgt then LET oldr2=w4fRgt ! PLOT w4Wndx(r1),wy; w4Wndx(oldr2),wy-stp ! ELSE IF r1w4fLft then ! PLOT w4Lft,wy; w4Wndx(oldr2),wy-stp ! END IF ! ELSE IF x<0 then ! IF r1>oldr2 then ! IF oldclr<>eq2clr then ! LET nose= r1 ! END IF ! SET COLOR green ! stable ! LET oldclr= eq2clr ! IF r1>w4fLft and r1w4fLft and oldr2w4fLft and r1w4fLft and r10 then ! IF r1w4fLft and r1w4fLft and oldr2eq0clr then LET nose= r1 ! SET COLOR red ! unstable ! LET oldclr= eq0clr ! IF r1>w4fLft and r1=0 then ! half attractor LET wy = w4Midy CALL EqPoint(wa,wy,0,red) LET wy = (w4Top+w4Midy)/2 CALL FlowPos(wa,wy) LET wy = (w4Midy+w4Bas)/2 CALL FlowNeg(wa,wy) END IF CASE 4 ! transcritical LET weq0= w4wndy(eq0) LET weq1= w4wndy(eq1) LET weq2= w4wndy(eq2) IF a<0 then ! SET COLOR phsclr LET wy = (w4Top+weq0)/2 CALL FlowNeg(wa,wy) LET wy = (w4Midy+weq1)/2 CALL FlowPos(wa,wy) LET wy = (w4Bas+weq1)/2 CALL FlowNeg(wa,wy) LET weq0= w4wndy(eq0) CALL EqPoint(wa,weq0,1,green) LET weq1= w4wndy(eq1) CALL EqPoint(wa,weq1,0,red) ELSE IF a>0 then ! half attractor SET COLOR phsclr LET wy = (w4Top+weq1)/2 CALL FlowNeg(wa,wy) LET wy = (w4Midy+weq1)/2 CALL FlowPos(wa,wy) LET wy = (w4Bas+weq0)/2 CALL FlowNeg(wa,wy) LET weq0= w4wndy(eq0) CALL EqPoint(wa,weq0,0,red) LET weq1= w4wndy(eq1) CALL EqPoint(wa,weq1,1,green) ELSE LET wy = (w4Top+w4Midy)/2 CALL FlowNeg(wa,wy) LET wy = (w4Midy+w4Bas)/2 CALL FlowNeg(wa,wy) LET weq1= w4wndy(eq0) CALL EqPoint(wa,weq1,.5,cyan) END IF CASE 5 ! 1 - ay LET weq0= w4wndy(eq0) ! LET weq1= w4wndy(eq1) ! LET weq2= w4wndy(eq2) IF a>-0.5 and a<0.5 then LET wy = (w4Top+w4y0)/2 CALL FlowPos(wa,wy) LET wy = (w4Bas+w4y0)/2 CALL FlowPos(wa,wy) ELSE IF a>0 then SET COLOR phsclr LET wy = (w4Top+weq0)/2 CALL FlowNeg(wa,wy) LET wy = (w4Bas+weq0)/2 CALL FlowPos(wa,wy) CALL EqPoint(wa,weq0,1,green) ELSE IF a<0 then SET COLOR phsclr LET wy = (w4Top+weq0)/2 CALL FlowPos(wa,wy) LET wy = (w4Bas+weq0)/2 CALL FlowNeg(wa,wy) CALL EqPoint(wa,weq0,0,red) END IF CASE 6 IF a>-1 and a<1 then ! left half bifurcation plane SET COLOR trjclr LET beq= w2Wndy(w2fTop/2) CALL FlowPos(wa,beq) LET beq= w2Wndy(w2fBas/2) CALL FlowPos(wa,beq) ELSE IF a= 1 then ! right half LET weq0= w2Wndy(eq0) CALL EqPoint(wa,weq0,-0.5,cyan) LET beq = w2Wndy((eq0+w2fTop)/2) CALL FlowPos(wa,beq) LET beq = w2Wndy((eq0+w2fBas)/2) CALL FlowPos(wa,beq) ELSE IF a=-1 then ! right half LET weq0= w2Wndy(eq0) CALL EqPoint(wa,weq0,-0.5,cyan) LET beq = w2Wndy((eq0+w2fTop)/2) CALL FlowPos(wa,beq) LET beq = w2Wndy((eq0+w2fBas)/2) CALL FlowPos(wa,beq) ELSE ! IF a>1 then ! right half LET weq1= w2Wndy(eq1) LET weq2= w2Wndy(eq2) CALL EqPoint(wa,weq1,1,green) CALL EqPoint(wa,weq2,0,red) LET beq = w2Wndy((eq2+w2fTop)/2) CALL FlowPos(wa,beq) LET beq = w2Wndy((eq2+eq1)/2) CALL FlowNeg(wa,beq) LET beq = w2Wndy((eq1+w2fBas)/2) CALL FlowPos(wa,beq) ! ELSE IF a<-1 then ! right half ! LET weq1= w2Wndy(eq1) ! LET weq2= w2Wndy(eq2) ! CALL EqPoint(wa,weq1,1,green) ! CALL EqPoint(wa,weq2,0,red) ! LET beq = w2Wndy((eq2+w2fTop)/2) ! CALL FlowPos(wa,beq) ! LET beq = w2Wndy((eq2+eq1)/2) ! CALL FlowNeg(wa,beq) ! LET beq = w2Wndy((eq1+w2fBas)/2) ! CALL FlowPos(wa,beq) END IF END SELECT IF plane=2 then CALL w2KeepGraphLayer END SUB SUB EqPoints(eq,a) ! 1 is attractor ! 0 is resistor ! 0.5 is half attractor ! -0.5 is half attractor LET eq0state= 0 LET eq1state= 0 LET eq2state= 0 SELECT CASE eq CASE 1 ! logistic with harvest LET n= 1 - 4*a LET d= -2 IF a<.25 then LET eq0= notanum LET eq1= (-1+sqr(n))/d LET eq2= (-1-sqr(n))/d LET eq1state= 0 LET eq2state= 1 ELSE IF a= 0.25 then LET eq0= 0.5 LET eq0state= 0.5 LET eq1= notanum LET eq2= notanum ELSE LET eq1,eq2,eq0= notanum END IF CASE 2 ! x^2 + a (saddle node) IF a<0 then LET eq0= notanum LET eq1= sqr(-a) LET eq2= -sqr(-a) LET eq1state= 0 LET eq2state= 1 ELSE IF a=0 then LET eq0= 0 LET eq0state= -.5 LET eq1= notanum LET eq2= notanum ELSE LET eq0,eq1,eq2= notanum END IF CASE 3 ! ax + x^3 (subcritical pitchfork) IF a<0 then LET eq1= sqr(-a) LET eq2= -sqr(-a) LET eq0= 0 LET eq0state= 1 LET eq1state= 0 LET eq2state= 0 ELSE IF a>=0 then LET eq0= 0 LET eq1= notanum LET eq2= notanum LET eq0state= 0 LET eq1state= 0 LET eq2state= 0 END IF CASE 4 ! ax - x^2 (transcritical) IF a>0 then LET eq0= 0 LET eq1= a LET eq2= notanum LET eq0state= 0 LET eq1state= 1 ELSE IF a<0 then LET eq0= 0 LET eq1= a LET eq2= notanum LET eq0state= 1 LET eq1state= 0 ELSE LET eq0= 0 LET eq1= notanum LET eq2= notanum LET eq0state= .5 END IF ! CASE 4 ! r + ax + x^3 ! IF a<0 then ! LET eq0= 0 ! LET eq1= sqr(-a) ! LET eq2= -sqr(-a) ! LET eq0state= 1 ! LET eq1state= 0 ! LET eq2state= 0 ! ELSE IF a>=0 then ! LET eq0= r ! LET eq1= notanum ! LET eq2= notanum ! LET eq0state= 0 ! END IF CASE 5 ! 1 + ay IF a=0 then LET eq0= notanum LET eq1= notanum LET eq2= notanum ! LET eq1state= 0 ! LET eq2state= 0 LET eq0state= 0 ELSE IF a>0 then LET eq0= 1/a LET eq1= notanum LET eq2= notanum LET eq0state= 1 ELSE IF a<0 then LET eq0= 1/a LET eq1= notanum LET eq2= notanum LET eq0state= 0 END IF CASE 6 ! x^2 - ax + .25 ! a=-1 b=a c=r ! -a^2/-4 + r ! LET vy= a^2/4 + r ! LET yvertex= -b^2/4*a + c LET vy= sgn(a)*(a^2/4 + 0.25) IF a>-1 and a<1 then LET eq0= notanum LET eq1= notanum LET eq2= notanum ELSE IF a=1 then LET eq0= vy LET eq0state= -0.5 LET eq1= notanum LET eq2= notanum ELSE IF a=-1 then LET eq0= vy LET eq0state= -0.5 LET eq1= notanum LET eq2= notanum ELSE LET n = sqr(a*a - 1) LET d = 2 LET eq0= notanum LET eq1= -(-a + n) / d LET eq2= -(-a - n) / d LET eq1state= 1 LET eq2state= 0 END IF ! CASE 6 ! r + ax - x^3 ! IF r=0 then ! LET eq0= 0 ! IF a>0 then ! LET eq1= +sqr(a) ! LET eq2= -sqr(a) ! LET eq1state= 1 ! LET eq2state= 1 ! LET eq0state= 0 ! ELSE ! LET eq1= notanum ! LET eq2= notanum ! LET eq0state= 1 ! END IF ! ELSE IF r>0 then ! CALL NewtonMethod(w4frgt,eq1) ! newton from frgt ! LET eq1state= 1 ! IF a>=nose then ! CALL NewtonMethod(w4flft,eq2) ! newton from flft ! CALL NewtonMethod(-.01,eq0) ! newton from -.01 ! LET eq0state= 0 ! LET eq2state= 1 ! ELSE ! LET eq2= notanum ! LET eq0= notanum ! END IF ! ELSE IF r<0 then ! CALL NewtonMethod(w4flft,eq2) ! newton from flft ! LET eq2state= 1 ! IF a>=nose then ! CALL NewtonMethod(w4frgt,eq1) ! newton from frgt ! CALL NewtonMethod(.01,eq0) ! newton from +.01 ! LET eq0state= 0 ! LET eq1state= 1 ! ELSE ! LET eq1= notanum ! LET eq0= notanum ! END IF ! END IF CASE else END SELECT END SUB ! SUB EqPoints ! IF r=0 then ! LET eq0= 0 ! IF a>0 then ! LET eq1= +sqr(a) ! LET eq2= -sqr(a) ! END IF ! ELSE IF r>0 then ! CALL NewtonMethod(frgt1,eq1) ! newton from frgt ! IF a>=nose then ! CALL NewtonMethod(flft1,eq2) ! newton from flft ! CALL NewtonMethod(-.01,eq0) ! newton from -.01 ! END IF ! ELSE IF r<0 then ! CALL NewtonMethod(flft1,eq2) ! newton from flft ! IF a>=nose then ! CALL NewtonMethod(frgt1,eq1) ! newton from frgt ! CALL NewtonMethod(.01,eq0) ! newton from +.01 ! END IF ! END IF ! END SUB SUB NewtonMethod(input,root) LET x1= input LET count= 0 DO LET x = x1 LET y = dxdt(m1Equation,a,x) LET d = df(x) IF d<>0 then LET x1= x - y/df(x) ! Newton ELSE LET x= 100000 EXIT DO ! store bignum in root?? END IF LET count= count+1 IF count>30 then EXIT DO LOOP until abs(y)<.01 LET root= x END SUB SUB GetImperfectNose ! needed to find nose value when w4 is turned off IF r<>0 then LET x = w4Fncy(w4Top) LET r1= (x*x*x-r)/x FOR wy= w4Top to w4Bas LET x = w4Fncy(wy) IF x<>0 then LET oldr1= r1 LET r1 = (x*x*x-r)/x IF r>0 and x<0 then ! upper right IF r1>oldr1 then IF oldclr<>eq2clr then LET nose= r1 LET oldclr= eq2clr ELSE LET oldclr= eq0clr END IF ELSE IF r<0 and x>0 then ! lower left IF r1eq0clr then LET nose= r1 LET oldclr= eq0clr END IF END IF END IF NEXT wy END IF END SUB SUB FlowPos(wx,wy) ! up arrow SET COLOR trajclr BOX AREA wx-1,wx+1,wy+6,wy-4 PLOT wx,wy-5 PLOT wx-5,wy; wx-2,wy; wx-2,wy-3; wx-5,wy PLOT wx+5,wy; wx+2,wy; wx+2,wy-3; wx+5,wy PLOT wx-3,wy-1 PLOT wx+3,wy-1 SET COLOR planeclr PLOT wx,wy-6 PLOT wx,wy+7 END SUB SUB FlowNeg(wx,wy) ! down arrow SET COLOR trajclr BOX AREA wx-1,wx+1,wy-6,wy+4 PLOT wx,wy+5 PLOT wx-5,wy; wx-2,wy; wx-2,wy+3; wx-5,wy PLOT wx+5,wy; wx+2,wy; wx+2,wy+3; wx+5,wy PLOT wx-3,wy+1 PLOT wx+3,wy+1 SET COLOR planeclr PLOT wx,wy-7 PLOT wx,wy+6 END SUB SUB EqPoint(eqx,eqy,fill,clr) ! equilibrium points on the bifurcation plane LOCAL sz,r1,r2,r3 LET sz= 4 LET r1= sz-3 LET r2= sz-2 LET r3= sz-1 SET COLOR black CALL BoxDisk(eqx-sz,eqx+sz,eqy+sz,eqy-sz) SET COLOR clr IF fill=1 then ! attractor CALL BoxDisk(eqx-sz,eqx+sz,eqy+sz,eqy-sz) BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz ELSE IF fill= 0 then ! repellor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 ELSE IF fill= 0.5 then ! positive half attractor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 PLOT eqx-r2,eqy-1; eqx+r2,eqy-1 PLOT eqx-r1,eqy-2; eqx+r1,eqy-2 ELSE IF fill= -0.5 then ! negative half attractor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 PLOT eqx-r2,eqy+1; eqx+r2,eqy+1 PLOT eqx-r1,eqy+2; eqx+r1,eqy+2 END IF END SUB SUB EqPointDE(eqx,eqy,fill,clr) ! equilibrium points on the DE plane ! half attractors are sideways for DE graph LOCAL sz,r1,r2,r3 LET sz= 4 LET r1= sz-3 LET r2= sz-2 LET r3= sz-1 SET COLOR black ! erase disk area CALL BoxDisk(eqx-sz,eqx+sz,eqy+sz,eqy-sz) SET COLOR clr IF fill=1 then ! attractor CALL BoxDisk(eqx-sz,eqx+sz,eqy+sz,eqy-sz) BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz ELSE IF fill= 0 then ! repellor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 ELSE IF fill= 0.5 then ! positive half attractor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 PLOT eqx+1,eqy-r2; eqx+1,eqy+r2 PLOT eqx+2,eqy-r1; eqx+2,eqy+r1 ELSE IF fill= -0.5 then ! negative half attractor BOX CIRCLE eqx-sz,eqx+sz,eqy+sz,eqy-sz BOX CIRCLE eqx-r3,eqx+r3,eqy+r3,eqy-r3 PLOT eqx-1,eqy-r2; eqx-1,eqy+r2 PLOT eqx-2,eqy-r1; eqx-2,eqy+r1 END IF END SUB END SUB ! --- end of phase lines tool code -------------------