! In the reverse engineering game, the fact that theta is sometimes ! undefined will cause us a problem. I hadn't thought of that before. ! Here's one solution: if your matrix happens to have a = d and b = -c, ! you are on or above the critical parabola. As you move around above ! that parabola, just keep the same property: in terms of T and D, ! ??? ! this is either ! w = sqrt(D - (T/2)^2) ! A = [ T/2, w ; -w , T/2 ] ! or ! A = [ T/2, -w ; w , T/2 ] ! (and you should keep the same sign) ! If (T/2)^2 >= D define ! w = sqrt((T/2)^2 - D) ! and ! A(s) = [ (T/2)-w s , 0 (T/2)+w ] ! If (T/2) < D define ! w3 = D - (T/2)^2 ! and ! A(s) = [ T/2 s , -w3/s T/2 ] ! ??? ! When you first go south of ! the parabola, we have to artifically choose a value of theta; theta = pi/2 ! is as good as any choice. The unavoidable effect will be that if the ! user picks a matrix below the parabola with theta different from pi/2 ! and moves it around using the (T,D) cursor, if he moves it above the ! parabola then it won't come back to the same matrix it started from ! when he gets back to the original (T,D) value. !! File: LinearPhasePortraits.bu !! September 9, 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$= "D'Arbeloff Interactive Math Project" LET colorscheme= 0 LET title$= "Linear Phase Portraits: Matrix Entry" SUB ThisProgram CALL ParameterPlane 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 ------ !! ----------------------------------------------------------- ! --- End TB5 header and subs --- SUB ParameterPlane DECLARE PUBLIC PCFlag,Mac5Flag,M68KFlag,UnixFlag,xmax,ymax DECLARE PUBLIC backclr,black,drkgry,drkmid,midgry,litmid,litgry,white DECLARE PUBLIC red,yellow,green,cyan,blue,magenta,colorscheme DECLARE PUBLIC planeclr,gridclr,rimclr,axisclr,axislabelclr,slideclr DECLARE PUBLIC workLft,workRgt,workBas,workTop,workMidx DECLARE PUBLIC true,false DECLARE DEF quitWithin, infoWithin ! --- help screen array --- DIM info1$(1:14) MAT READ info1$ DATA "Linear Phase Portraits: Matrix Entry" DATA "" DATA "This illustration supports visual exploration of the relationship between the elements a, b, c and d of matrix A, the trace and determinant of matrix A, the corresponding vector field, and the trajectories of a 2-dimensional homogeneous system:" DATA "-dx/dt = ax + by" DATA "-dy/dt = cx + dy" DATA "or in vector terms" DATA "-x' = Ax" DATA "" DATA "The direction field, eigenlines, and trajectories are displayed on the x,y phase plane to the right." DATA "Click the phase plane to set initial values and graph the corresponding trajectory." DATA "" DATA "The matrix variable, A, designates the square matrix of constant coefficients [a,b, c,d]." DATA "The trace, Tr, of matrix A is a+d, and the determinant, Det, of matrix A is ad - bc." DATA "The trace/determinant plane on the left is divided into regions corresponding to the geometry and stability of the phase portraits of the differential equations. The color of the eigenlines and trajectory match the color of the region of the point on the (tr,det) plane." ! DATA "Click or drag the mouse on the trace/determinant plane to set the trace and determinant values and use to redefine the A matrix constants, to reclassify the geometry and stability of the phase portrait, and to redraw the corresponding direction field and eigenlines in the phase plane." DIM info2$(1:9) MAT READ info2$ DATA "Linear Phase Portraits: Matrix Entry" DATA "" ! DATA "If the companion matrix checkbox is on, the matrix shown is A = [0,1 -Det A,Tr A]." ! DATA "If the companion matrix checkbox is off, the matrix shown is A = [more complicated]." ! DATA "" DATA "Click or drag the A matrix sliders, a, b, c, or d, to set the matrix elements. The trace and determinate of the A matrix are computed as the sliders are moved, and the trace/determinate cursor is shown on the Tr,Det plane." DATA "The trace and determinate values are then used to reclassify the geometry and stability of the phase portrait, and to redraw the corresponding direction field and eigenlines on the phase plane." DATA "Click the Tr,Det 'draw' checkbox to turn the cursor drawing option on or off." DATA "When 'draw' is on, the cursor draws a trail on the Tr,Det plane when the sliders are moved. Click the Tr,Det [clear] button to remove cursor trails from the Tr,Det plane." DATA "" DATA "Click the phase plane [clear] button to erase the phase plane trajectories." DATA "Click the [eigenvalues] button at the lower left to toggle the complex eigenvalue plane on and off. When the complex plane is on, the eigenvalues of the characteristic equation are displayed as the matrix is changed." ! This tool supports the visual exploration of the relationship between the matrix A, its trace and determinant, the corresponding vector field, and the trajectories of solutions to x' = Ax. ! ! The matrix entries are controlled by sliders at lower right and the matrix A is displayed at left. ! Above that is a readout of the trace and determinant of A. ! At upper left the trace/determinant plane is displayed, divided into regions corresponding to geometric type of the phase portrait of the differential equation. ! If tr A and det A are both between -4 and 4, the point is indicated in that window by a "+." ! ! The matrix may also be controlled by a mouse click or drag in that window. ! In this mode, the matrix will either be the companion matrix, with top row [ 0 1 ] , ! or the matrix with the currently active angular and asymetry properties (which are then fixed ! until the matrix entries are changed manually or the companion matrix is selected). ! ! At upper right the direction field is displayed, along with all eigenlines. ! Rolling over that window gives a readout of the coordinates as an initial value, ! and clicking in it creates the trajectory through that point. ! The color of the trajectory and the color of the eigenlines match the color of the region in the (tr,det) plane. ! Below this window the name of the geometric and stability type is displayed. ! The "clear" button clears all trajectories except those along eigenlines. ! The "eigenvalues" button at bottom toggles a window displaying the complex plane with the eigenvalues of A marked, and a readout of the two eigenvalues. ! trace = a+d ! dterm = ad - bc ! characteristic polynomial = lam^2 - trace*lam + dterm ! lam = (trace +- sqr(trace^2 - 4 dterm)) / 2 ! dscrm = trace^2 - 4 dterm ! dterm = trace^2/4 is the parabola on the trace,dterm plane ! IF dterm0 ! outside the parabola - two real lams ! the theta,s plane is continous ! IF dterm=trace^2/4, dscrm=0 ! on the parabola - one real lam = trace/2 ! the theta,s plane is continous - the horizontal axis matters ! IF dterm>trace^2/4, dscrm<0 ! inside the parabola - two complex lams ! the theta,s plane has two active regions ! s^2 must be >= (t2^2 - dterm) ! ----- modify the palette ----- LET b0= 0 ! web safe intervals LET b1= 0.20 LET b2= 0.40 LET b3= 0.60 LET b4= 0.80 LET b5= 1.00 SET COLOR MIX( 7) b5,b3,b1 ! red SET COLOR MIX( 8) b5,b5,b0 ! yellow SET COLOR MIX( 9) b5,b2,b1 ! green SET COLOR MIX(10) b0,b4,b3 ! cyan SET COLOR MIX(11) b2,b3,b5 ! blue SET COLOR MIX(12) b0,b4,b1 ! magenta SET COLOR MIX(13) b5,b1,b5 ! pink FOR clr= 7 to 13 SET COLOR clr NEXT clr LET red = 7 LET yellow = 8 LET green = 9 LET cyan = 10 LET blue = 11 LET magenta= 12 LET pink = 13 ! SET COLOR MIX( 7) b5,b3,b1 ! red ! SET COLOR MIX( 8) b5,b5,b0 ! yellow ! SET COLOR MIX( 9) b5,b2,b1 ! green ! SET COLOR MIX(10) b0,b4,b3 ! cyan ! SET COLOR MIX(11) b2,b3,b5 ! blue ! SET COLOR MIX(12) b0,b4,b1 ! magenta ! SET COLOR MIX(13) b5,b1,b5 ! pink ! FOR clr= 7 to 13 ! SET COLOR clr ! NEXT clr ! ! LET red = 7 ! LET yellow = 8 ! LET green = 9 ! LET cyan = 10 ! LET blue = 11 ! LET magenta= 12 ! LET pink = 13 LET axisclr= litgry ! ---------- Utility functions --- DECLARE DEF clamp,roundn,e DEF notanum= 987656789 ! --- functions --- DEF dxdt(x,y)= ma*x + mb*y DEF dydt(x,y)= mc*x + md*y SUB Complex(trace,dterm,dscrm,cmplx) LET dscrm= trace*trace - 4*dterm ! find discriminant IF dscrm<0 then LET cmplx= 1 else LET cmplx= 0 ! complex? END SUB SUB LinearParams(a,b,c,d,trace,dterm,dscrm,cmplx) LET trace= a+d LET dterm= a*d - b*c LET dscrm= trace*trace - 4*dterm ! find discriminant IF dscrm<0 then LET cmplx= 1 else LET cmplx= 0 ! complex? END SUB DEF fT(ma,mb,mc,md) = ma+md DEF fD(ma,mb,mc,md) = ma*md - mb*mc DEF fDsc(trace,dterm)= trace*trace - 4*dterm SUB MatToTD(ma,mb,mc,md,trace,dterm) LET trace= fT(ma,mb,mc,md) ! find trace LET dterm= fD(ma,mb,mc,md) ! find determinant END SUB DEF omIV(dterm,trace) LET omIV= sqr(abs(dterm - trace^2/4)) END DEF DEF om(dterm,trace) IF dscrm<0 then LET om= sqr(abs(dterm - trace^2/4)) ELSE LET om= 0 END IF END DEF ! ---------- Graphing plane parameters and methods ---------- ! --- w1 plane data - trace,dterminant --- 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 w1xPiFlag, w1xMult, w1yPiFlag, w1yMult DECLARE PUBLIC w1wWid, w1fWid, w1wHgt, w1fHgt, w1Aspect LET w1Flag = 1 ! pixel window visibility LET w1xPiFlag= 0 ! pi switch for x axis LET w1xMult = 1 LET w1yPiFlag= 0 ! pi switch for y axis LET w1yMult = 1 LET w1Lft = workLft+60 ! pixel bounds LET w1Rgt = w1Lft+200 ! for 6 LET w1Top = workTop+45 LET w1Bas = w1Top+200 LET w1size= 8 LET w1fLft= -4 ! function bounds LET w1fRgt= 4 LET w1fBas= -4 LET w1fTop= 4 LET w1Xax$= "Tr A" ! axis labels LET w1Yax$= "Det A" LET w1xGridstep= 0 ! horizontal grid intervals LET w1yGridstep= 0 ! vertical grid intervals LET w1xSTik = 1/2 ! horizontal axis Tik marks LET w1xLTik = 1 LET w1xLabel= 1 LET w1xFirst= w1fLft LET w1ySTik = 1/2 ! vertical axis Tik marks LET w1yLTik = 1 LET w1yLabel= 2 LET w1yFirst= w1fBas ! --- w1 methods --- DECLARE DEF w1fncx,w1fncy,w1wndx,w1wndy ! window/function transforms DECLARE DEF w1Within, w1wWithin CALL w1Variables SUB w1Init LET aspect= (w1wwid/w1fwid)/(w1whgt/w1fhgt) CALL w1DrawPlane(1,1,1) ! grid, axes, zeroaxes CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(w1Rgt+8,w1y0+3,w1Xax$,axislabelclr) ! axis labels CALL PlotTextCJ(w1x0,w1Top-10,w1Yax$,axislabelclr) CALL DrawTDPlane CALL w1KeepGridLayer CALL w1KeepGraphLayer END SUB ! --- w2 data: eigenvalue plane --- DECLARE PUBLIC w2Lft, w2Rgt, w2Bas, w2Top, w2Midx, w2Midy DECLARE PUBLIC w2fLft, w2fRgt, w2fBas, w2fTop, w2x0, w2y0 DECLARE PUBLIC w2xFirst, w2xSTik, w2xLTik, w2xLabel, w2xGridstep DECLARE PUBLIC w2yFirst, w2ySTik, w2yLTik, w2yLabel, w2yGridstep DECLARE PUBLIC w2xPiFlag, w2xMult, w2yPiFlag, w2yMult DECLARE PUBLIC w2wWid, w2fWid, w2wHgt, w2fHgt, w2Aspect LET w2Flag = 1 ! pixel window visibility LET w2xPiFlag= 0 ! pi switch for x axis LET w2xMult = 1 LET w2yPiFlag= 0 ! pi switch for y axis LET w2yMult = 1 LET w2Lft = w1Lft ! pixel bounds LET w2Rgt = w2Lft + 80 LET w2Top = workBas - 110 LET w2Bas = w2Top + 80 LET w2fLft = -4 ! function bounds LET w2fRgt = 4 LET w2fBas = -4 LET w2fTop = 4 LET w2xAx$ = "Re" ! axis labels LET w2yAx$ = "Im" LET w2xGridstep= 0 ! grid line intervals LET w2yGridstep= 0 LET w2xSTik = 1 ! axis Tik marks LET w2xLTik = 2 LET w2xLabel= 2 LET w2xFirst= w2fLft LET w2ySTik = 1 LET w2yLTik = 2 LET w2yLabel= 2 LET w2yFirst= w2fBas ! --- w2 methods --- DECLARE DEF w2fncx,w2fncy,w2wndx,w2wndy ! window/function transforms CALL w2Variables SUB w2Clear BOX CLEAR w2Lft-25,w2Rgt+25,w2Bas+25,w2Top-25 END SUB SUB w2Init CALL w2DrawPlane(1,1,1) ! grid, axes, zeroaxes CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(w2Rgt+8,w2y0+3,w2xAx$,axislabelclr) ! axis labels CALL PlotTextCJ(w2x0,w2Top-10,w2yAx$,axislabelclr) CALL w2KeepGridLayer CALL w2KeepGraphLayer END SUB ! --- w3 data: phase plane --- DECLARE PUBLIC w3Lft,w3Rgt,w3Bas,w3Top,w3Midx,w3Midy DECLARE PUBLIC w3fLft,w3fRgt,w3fBas,w3fTop,w3x0,w3y0 DECLARE PUBLIC w3xFirst, w3xSTik, w3xLTik, w3xLabel, w3xGridstep DECLARE PUBLIC w3yFirst, w3ySTik, w3yLTik, w3yLabel, w3yGridstep DECLARE PUBLIC w3xPiFlag, w3xMult, w3yPiFlag, w3yMult DECLARE PUBLIC w3wWid, w3fWid, w3wHgt, w3fHgt, w3Aspect LET w3Flag = 1 LET w3xPiFlag= 0 LET w3xMult = 1 LET w3yPiFlag= 0 LET w3yMult = 1 LET w3Lft = workRgt-385 ! pixel bounds LET w3Rgt = w3Lft + 300 LET w3Top = w1Top LET w3Bas = w3Top + 300 LET w3Size= 3 LET w3fLft= -w3Size ! function bounds * pi LET w3fRgt= w3Size LET w3fBas= -w3Size LET w3fTop= w3Size LET w3xGridstep= 0 ! horizontal grid intervals LET w3yGridstep= 0 ! vertical grid intervals LET w3xAx$= "x" ! axis labels LET w3yAx$= "y" LET w3xSTik = 1/2 ! horizontal axis Tik marks LET w3xLTik = 1 LET w3xLabel= 1 LET w3xFirst= w3fLft LET w3ySTik = 1/2 ! vertical axis Tik marks LET w3yLTik = 1 LET w3yLabel= 1 LET w3yFirst= w3fBas ! --- w3 methods --- DECLARE DEF w3fncx, w3fncy, w3wndx, w3wndy ! window/function transforms DECLARE DEF w3Within, w3wWithin, w3fWithin CALL w3Variables SUB w3Init CALL w3DrawPlane(1,1,1) ! grid, axes, zeroaxes CALL SetTextFont(1,12,"normal") CALL PlotTextLJ(w3Rgt+8,w3y0+3,w3xAx$,axislabelclr) ! axis labels CALL PlotTextCJ(w3x0,w3Top-11,w3yAx$,axislabelclr) CALL SetMatDiffEq CALL w3KeepGridLayer CALL w3KeepGraphLayer END SUB SUB w3KeepFieldLayer BOX KEEP w3Lft-5,w3Rgt+5,w3Bas+5,w3Top-5 in w3FieldLayer$ END SUB SUB w3ShowFieldLayer BOX SHOW w3FieldLayer$ at w3Lft-5,w3Bas+5 END SUB SUB SetMatDiffEq LET DEBas= w1Top - 28 CALL SetTextFont(1,12,"bold") LET txt$ = "x = Ax" CALL StringWidth(txt$,sw) LET DELft= w3x0 - sw/2 CALL PlotTextLJ(DELft,DEBas,txt$,litgry) CALL PlotTextLJ(DELft+2,DEBas-9,".",white) END SUB ! ------ matrices ------ SUB MatrixBrackets(ml,mr,mb,mt) SET COLOR litmid ! matrix brackets PLOT ml+5,mt ; ml,mt ; ml,mb ; ml+5,mb PLOT ml+5,mt+1; ml+1,mt+1; ml+1,mb-1; ml+5,mb-1 PLOT mr-5,mt ; mr,mt ; mr,mb ; mr-5,mb PLOT mr-5,mt+1; mr-1,mt+1; mr-1,mb-1; mr-5,mb-1 END SUB ! ------ X matrix ------ LET lnspc= 15 LET m0Top = w1Bas + 80 LET m0Bas = m0Top + 40 LET m0Lft = w1Lft + 150 LET m0Rgt = m0Lft + 75 LET m1txt0= m0Top + 14 LET m1txt1= m0Top + 33 LET m1col1= m0Lft + 10 LET m0midx= int((m0Lft+m0Rgt)/2) ! ------ A matrix ------ LET lnspc= 15 LET mTop = w1Bas + 90 LET mBas = mTop + 40 LET mLft = w1Lft LET mRgt = mLft + 115 LET mtxt0= mTop + 14 LET mtxt1= mTop + 33 LET mcol1= mLft + 10 LET mcol2= mcol1 + 50 LET mmidx= int((mLft+mRgt)/2) LET mmidy= int((mTop+mBas)/2) SUB DrawMatrix CALL PlotTextLJ(mLft-20,mmidy+3,"A=",litgry) CALL MatrixBrackets(mLft,mRgt,mBas,mTop) END SUB SUB ClearMatrix BOX CLEAR mLft-30,mRgt+10,mBas+10,mTop-10 END SUB SUB ShowMatrixValues CALL SetTextFont(1,12,"bold") LET matclr= litgry CALL ClearMatrixValues CALL PlotTextLJ(mcol1,mtxt0,using$("--%.##",ma),matclr) CALL PlotTextLJ(mcol2,mtxt0,using$("--%.##",mb),matclr) CALL PlotTextLJ(mcol1,mtxt1,using$("--%.##",mc),matclr) CALL PlotTextLJ(mcol2,mtxt1,using$("--%.##",md),matclr) END SUB SUB ClearMatrixValues BOX CLEAR mLft+3,mRgt-3,mtxt0+3,mtxt0-10 ! upper row BOX CLEAR mLft+3,mRgt-3,mtxt1+3,mtxt1-10 ! lower row END SUB SUB MatAReset(ma) CALL SetTextFont(1,12,"bold") BOX CLEAR mcol1-2,mcol2-5,mtxt0+5,mtxt0-10 CALL PlotTextLJ(mcol1,mtxt0,using$("--%.##",ma),matclr) END SUB SUB MatBReset(mb) CALL SetTextFont(1,12,"bold") BOX CLEAR mcol2-2,mRgt-5,mtxt0+5,mtxt0-10 CALL PlotTextLJ(mcol2,mtxt0,using$("--%.##",mb),matclr) END SUB SUB MatCReset(mc) CALL SetTextFont(1,12,"bold") BOX CLEAR mcol1-2,mcol2-5,mtxt1+5,mtxt1-10 CALL PlotTextLJ(mcol1,mtxt1,using$("--%.##",mc),matclr) END SUB SUB MatDReset(md) CALL SetTextFont(1,12,"bold") BOX CLEAR mcol2-2,mRgt-5,mtxt1+5,mtxt1-10 CALL PlotTextLJ(mcol2,mtxt1,using$("--%.##",md),matclr) END SUB SUB MatResetAll(a,b,c,d) CALL MatAReset(a) CALL MatBReset(b) CALL MatCReset(c) CALL MatDReset(d) END SUB ! --- matrix sliders --- LET matslidebound= 4 LET matslidewid = 120 ! --- h1 slider --- DECLARE PUBLIC h1axis,h1wLft,h1wRgt,h1wBas,h1wTop,h1fLft,h1fRgt DECLARE PUBLIC h1name$,h1form$,h1clr,h1First,h1STik,h1LTik,h1Label DECLARE PUBLIC h1PiAxis,h1Mult,h1fMin,h1fMax DECLARE DEF h1Fncx, h1Within ! window/function transforms LET h1PiAxis= 0 LET h1Mult = 1 LET h1clr = slideclr LET h1name$ = "a" LET h1form$ = "-%.##" LET h1Places= 2 LET h1Click = 1/2 LET h1axis = w3Bas + 115 LET h1wLft = w3Lft LET h1wRgt = h1wLft+matslidewid LET h1fLft = -matslidebound LET h1fRgt = matslidebound LET h1STik = 1/2 ! short tick marks LET h1LTik = 1 ! long tick marks LET h1Label = matslidebound ! labels LET h1First = h1fLft ! first tick mark CALL h1SliderVariables SUB h1Init CALL h1DrawSlider(h1name$,a) END SUB ! --- h2 Slider --- DECLARE PUBLIC h2axis,h2wLft,h2wRgt,h2wBas,h2wTop,h2fLft,h2fRgt DECLARE PUBLIC h2name$,h2form$,h2clr,h2First,h2STik,h2LTik,h2Label DECLARE PUBLIC h2PiAxis,h2Mult,h2fMin,h2fMax DECLARE DEF h2Fncx, h2Within ! window/function transforms LET h2PiAxis= 0 LET h2Mult = 1 LET h2clr = slideclr LET h2name$ = "b" LET h2form$ = "-%.##" LET h2Places= 2 LET h2Click = 1/2 LET h2axis = h1axis LET h2wLft = h1wRgt + 85 LET h2wRgt = h2wLft + matslidewid LET h2fLft = -matslidebound LET h2fRgt = matslidebound LET h2STik = 1/2 LET h2LTik = 1 LET h2Label = matslidebound LET h2First = h2fLft CALL h2SliderVariables SUB h2Init CALL h2DrawSlider(h2name$,b) END SUB ! --- h3 slider --- 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 h3Fncx, h3Within ! window/function transforms LET h3PiAxis= 0 LET h3Mult = 1 LET h3clr = slideclr LET h3name$ = "c" LET h3form$ = "-%.##" LET h3Places= 2 LET h3axis = h1axis + 45 LET h3wLft = h1wLft LET h3wRgt = h1wRgt LET h3fLft = -matslidebound LET h3fRgt = matslidebound LET h3STik = 1/2 LET h3LTik = 1 LET h3Label = matslidebound LET h3First = h3fLft LET h3Click = 1/2 CALL h3SliderVariables SUB h3Init CALL h3DrawSlider(h3name$,c) END SUB ! --- h4 slider --- 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 h4Fncx, h4Within ! window/function transforms LET h4PiAxis= 0 LET h4Mult = 1 LET h4clr = slideclr LET h4name$ = "d" LET h4form$ = "-%.##" LET h4Places= 2 LET h4axis = h3axis LET h4wLft = h2wLft LET h4wRgt = h2wRgt LET h4fLft = -matslidebound LET h4fRgt = matslidebound LET h4STik = 1/2 LET h4LTik = 1 LET h4Label = matslidebound LET h4First = h4fLft LET h4Click = 1/2 CALL h4SliderVariables SUB h4Init CALL h4DrawSlider(h4name$,d) END SUB ! -------- text locations ------ ! ----- t1 - w3 initial values ----- LET t1Lft = w3Midx + 60 LET t1Rgt = w3Rgt LET t1Bas1 = w3Bas + 39 LET t1Bas2 = t1Bas1 + 25 LET t1Top = t1Bas1 - 15 LET t1Bas = t1Bas2 + 5 LET t1Label1$= "x(0) = " LET t1Label2$= "y(0) = " CALL StringWidth(t1Label1$,sw) LET t1Eqx = t1Lft + sw LET t1clr = litgry SUB t1Init CALL SetTextFont(1,12,"normal") CALL PlotTextRJ(t1Eqx,t1Bas1,t1Label1$,t1clr) CALL PlotTextRJ(t1Eqx,t1Bas2,t1Label2$,t1clr) END SUB SUB t1Clear BOX CLEAR t1Lft-2,t1Rgt,t1Bas,t1Top END SUB SUB t1SetValues(x0,y0) CALL t1ClearValues LET x$= using$("-%.##",x0) CALL PlotTextRJ(t1Rgt,t1Bas1,x$,t1clr) LET y$= using$("-%.##",y0) CALL PlotTextRJ(t1Rgt,t1Bas2,y$,t1clr) END SUB SUB t1ClearValues BOX CLEAR t1Eqx,t1Rgt+10,t1Bas,t1Top END SUB ! ----- t2 - classification name ----- LET t2Lft = w3Lft LET t2Rgt = t1Lft-2 LET t2BasLn= t1Bas2 LET t2Top = t2BasLn - 12 LET t2Bas = t2BasLn + 5 SUB t2SetName(name$) CALL t2Clear CALL SuperSubScriptLJ(t2Lft,t2BasLn,name$,tclr) END SUB SUB t2Clear BOX CLEAR t2Lft-2,t2Rgt,t2Bas,t2Top END SUB ! ----- t3 - trace and det ----- LET t3Lft1 = mrgt + 20 LET t3Lft2 = t3Lft1 LET t3Rgt = w3Lft - 20 LET t3Bas1 = mtxt0 LET t3Bas2 = mtxt1 LET t3Top = t3Bas1 - 12 LET t3Bas = t3Bas2 + 5 LET t3eqx = t3Lft1+60 CALL StringWidth(" -44.44",sw) LET t3txtRgt= t3eqx + sw LET t3Clr = litgry SUB t3Init CALL t3Clear CALL SetTextFont( 1,12,"bold") CALL PlotTextRJ(t3eqx,t3Bas1,"Tr A = ",t3Clr) CALL PlotTextRJ(t3eqx,t3Bas2,"Det A = ",t3Clr) CALL t3SetValue END SUB SUB t3Clear BOX CLEAR t3Lft1-2,t3Rgt,t3Bas,t3Top END SUB SUB t3SetValue CALL t3ClearValue CALL SetTextFont( 1,12,"bold") CALL PlotTextRJ(t3txtRgt,t3Bas1,using$("--%.##",trace),t3Clr) CALL PlotTextRJ(t3txtRgt,t3Bas2,using$("--%.##",dterm),t3Clr) END SUB SUB t3ClearValue BOX CLEAR t3eqx-2,t3Rgt,t3Bas,t3Top END SUB ! ----- t4 - eigenvalues ----- LET t4Lft = w2Rgt + 40 LET t4Rgt = t4Lft + 155 LET t4Bas1= w2Midy - 5 LET t4Bas2= t4Bas1 + 18 LET t4Bas = t4Bas2 + 5 LET t4Top = t4Bas1 - 12 LET t4Eqx = t4Lft + 44 LET t4Clr = litgry SUB t4Init CALL DrawLam(t4Lft+14,t4Bas1,t4Clr) CALL DrawLam(t4Lft+14,t4Bas2,t4Clr) CALL SetTextFont(1,9,"bold") CALL PlotTextLJ(t4Lft+23,t4Bas1+4,"1",t4Clr) CALL PlotTextLJ(t4Lft+23,t4Bas2+4,"2",t4Clr) CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(t4Lft+32,t4Bas1,"=",t4Clr) CALL PlotTextLJ(t4Lft+32,t4Bas2,"=",t4Clr) CALL t4SetValues END SUB SUB t4SetValues CALL FindLamdas CALL t4ClearValues CALL SetTextFont(1,12,"bold") CALL PlotTextLJ(t4Eqx,t4Bas1,l1$,t4Clr) CALL PlotTextLJ(t4Eqx,t4Bas2,l2$,t4Clr) END SUB SUB t4Clear BOX CLEAR t4Lft-2,t4Rgt,t4Bas,t4Top END SUB SUB t4ClearValues BOX CLEAR t4Eqx,t4Rgt,t4Bas,t4Top END SUB ! ---- classification name ---- LET lclft= w3lft LET lcbas= t1Bas2 SUB ClassText(className$) CALL ClassTextClear CALL SuperSubScriptLJ(lcLft,lcBas,className$,classClr) !CALL InitValLabels END SUB SUB ClassTextClear BOX CLEAR lcLft,lcLft+170,lcbas+5,lcbas-15 ! className$ END SUB ! --- checkbox - companion matrix --- DECLARE PUBLIC cb1Lft,cb1Rgt,cb1Bas,cb1Top,cb1Txt$,cb1Clr,cb1State DECLARE DEF cb1Within LET cb1Lft = w1Lft LET cb1Bas = w1Bas+40 LET cb1Txt$= "Companion Matrix" LET cb1Clr = litgry CALL cb1Variables ! --- checkbox - drawstate --- DECLARE PUBLIC cb2Lft,cb2Rgt,cb2Bas,cb2Top,cb2Txt$,cb2Clr,cb2State DECLARE DEF cb2Within LET cb2Lft = w1Lft LET cb2Bas = w1Bas+40 LET cb2Txt$= "draw" LET cb2Clr = litgry CALL cb2Variables ! --- buttons -- LET bhgt= 19 LET cTop= w3Bas + 25 LET cBas= cTop + bhgt LET cLft= w3Lft LET cRgt= cLft + 50 LET w1cTop= cb2Bas + 10 LET w1cBas= w1cTop + bhgt LET w1cLft= w1Lft LET w1cRgt= w1cLft + 50 LET evbTop= workBas - 20 LET evbBas= evbTop + bhgt LET evbLft= w2Rgt + 55 LET evbRgt= evbLft + 95 SUB w1ButtonDraw CALL SetTextFont(1,12,"bold") CALL DrawButton(w1cLft,w1cRgt,w1cBas,w1cTop,6,"clear") END SUB SUB w1ButtonClear BOX CLEAR w1cLft,w1cRgt,w1cBas,w1cTop END SUB SUB Buttons CALL SetTextFont(1,12,"bold") CALL DrawButton(cLft,cRgt,cBas,cTop,6,"clear") CALL DrawButton(evbLft,evbRgt,evbBas,evbTop,6,"eigenvalues") END SUB ! ----- default parameters ----- IF M68kFlag=1 then LET tstep= 1/128 ELSE LET tstep= 1/64 END IF LET ma,a,oldma = 2 LET mb,b,oldmb = 1 LET mc,c,oldmc = 1 LET md,d,oldmd = -1 ! ----- set default vars and initialize ----- IF M68kFlag=1 then ! mac TB 4.04 LET tstep= 1/128 ELSE LET tstep= 1/32 END IF LET trace,oldtrace= 0 LET dterm,olddterm= -1 LET s,olds = 0 LET theta,oldtheta= pi/2 LET omega = 0 LET close = 0.1 ! proximity to parabola LET companion = 0 LET w1drawState = 0 ! ----- initialize ----- CALL LinearParams(a,b,c,d,trace,dterm,dscrm,cmplx) CALL InitScreen SUB InitScreen BOX CLEAR workLft,workRgt,workBas,workTop CALL w1Init CALL w3Init IF lamflag=1 then CALL w2Init CALL t4Init END IF CALL DrawMatrix CALL h1Init CALL h2Init CALL h3Init CALL h4Init CALL Buttons CALL t1Init CALL t3Init CALL cb2Init CALL ResetDisplay END SUB ! --- event manager --- DO LET clearflag= 0 DO GET MOUSE: mx,my,ms IF w3wWithin(mx,my)=true then IF mx<>omx or my<>omy then CALL w3Rollover LET omx= mx LET omy= my END IF ELSE IF clearflag2=1 then CALL t1ClearValues CALL w3RolloverClear LET clearflag2= 0 END IF LOOP until ms=2 ! IF w1Within(mx,my)=true then ! trace det plane ! CALL w1wClamp(mx,my) ! CALL w3ShowGridLayer ! LET omx,omy= 999 ! DO ! GET MOUSE: mx,my,ms ! IF w1Within(mx,my)=true then ! CALL w1DragEvent ! END IF ! LOOP until ms=3 IF w3Within(mx,my)=true then ! phase plane ! CALL MouseUp(mx,my,ms) ! LET x0= w3fncx(mx) ! LET y0= w3fncy(my) ! CALL t1SetValues(x0,y0) ! CALL Trajectory(x0,y0) LET omx,omy= -999 DO GET MOUSE: mx,my,ms IF w3Within(mx,my)=true then CALL w3wClamp(mx,my) IF omx<>mx or omy<>my then CALL w3PixelsToMath(mx,my,x0,y0) CALL w3ShowGraphLayer CALL SetVector(x0,y0,x2,y2,.2) CALL SetVars(mx,my,omx,omy) END IF END IF LOOP until ms=3 CALL w3ShowGraphLayer CALL Trajectory(x0,y0) CALL w3KeepGraphLayer ELSE IF h1Within(mx,my)=true then LET olda= -999 IF myevbLft and mxevbTop and myw1cLft and mxw1cTop and mycLft and mxcTop and my0 or dy<>0) then CALL DrawVectorArrow(angle(dx,-dy),6,wx2,wy2,classClr) END IF END SUB ! --- SUB w1DragEvent CALL w1wClamp(mx,my) IF omx<>mx or omy<>my then CALL w1PixelsToMath(mx,my,trace,dterm) LET trace= round(trace,4) LET dterm= round(dterm,4) CALL TDProximity IF abs(dterm)<=0.1 then LET dterm= 0 IF abs(trace)<=0.1 then LET trace= 0 CALL TraceDetUpdate CALL t3SetValue IF lamflag=1 then CALL t4SetValues CALL CopyCoords(mx,my,omx,omy) END IF END SUB SUB TraceDetUpdate CALL TraceDetToMatrix(trace,dterm,theta,s) CALL MarkTDCursor(trace,dterm) CALL ShowMatrixValues CALL SliderResetAll(ma,mb,mc,md) CALL Classification(trace,dterm,dscrm,className$,classClr) CALL DirectionField CALL FindEigenVectors(ma,mb,mc,md) CALL w3KeepGraphLayer END SUB ! SUB w3Rollover ! IF omx<>mx or omy<>my then ! LET x0 = w3fncx(mx) ! LET y0 = w3fncy(my) ! CALL t1SetValues(x0,y0) ! LET clearflag2= 1 ! LET omx= mx ! LET omy= my ! END IF ! END SUB SUB w3Redraw CALL DirectionField CALL FindEigenVectors(ma,mb,mc,md) END SUB SUB DtermProximity LET parab= trace^2/4 ! near parabola? IF abs(dterm-parab)=0 then ! upper half plane? ! check for parabola LET signtr= sgn(trace) LET x = sqr(4*dterm) IF abs(abs(trace)-x)=0 then LET s= omega ELSE LET s= -omega END IF IF s<>olds then LET olds= s END IF END IF ELSE LET omega= 0 END IF LET L = sqr(s^2 + td) LET Lcs= L*cos(2*theta) LET Lsn= L*sin(2*theta) LET ma = t2 - Lsn ! define matrix elements LET mb = Lcs + s LET mc = Lcs - s LET md = t2 + Lsn END IF END SUB ! ---- slider events ---- ! --- h1 slider events --- SUB h1MouseClick CALL h1GetClickVal(ms,h1Click,ma) CALL h1MouseAction END SUB SUB h1MouseDrag DO CALL h1GetDragVal(ms,h1Places,ma) CALL h1MouseAction LOOP until ms=3 END SUB SUB h1MouseAction IF ma<>oldma then LET a,oldma= ma CALL SliderInputAction CALL MatAReset(ma) END IF END SUB ! --- h2 slider events --- SUB h2MouseClick CALL h2GetClickVal(ms,h2Click,mb) CALL h2MouseAction END SUB SUB h2MouseDrag DO CALL h2GetDragVal(ms,h2Places,mb) CALL h2MouseAction LOOP until ms=3 END SUB SUB h2MouseAction IF mb<>oldmb then LET b,oldmb= mb CALL SliderInputAction CALL MatBReset(mb) END IF END SUB ! --- h3 slider events --- SUB h3MouseClick CALL h3GetClickVal(ms,h3Click,mc) CALL h3MouseAction END SUB SUB h3MouseDrag DO CALL h3GetDragVal(ms,h3Places,mc) CALL h3MouseAction LOOP until ms=3 END SUB SUB h3MouseAction IF mc<>oldmc then LET c,oldmc= mc CALL SliderInputAction CALL MatCReset(mc) END IF END SUB ! --- h4 slider events --- SUB h4MouseClick CALL h4GetClickVal(ms,h4Click,md) CALL h4MouseAction END SUB SUB h4MouseDrag DO CALL h4GetDragVal(ms,h4Places,md) CALL h4MouseAction LOOP until ms=3 END SUB SUB h4MouseAction IF md<>oldmd then LET d,oldmd= md CALL SliderInputAction CALL MatDReset(md) END IF END SUB ! --- SUB SliderInputAction CALL LinearParams(ma,mb,mc,md,trace,dterm,dscrm,cmplx) CALL Classification(trace,dterm,dscrm,className$,classClr) ! CALL MatToThetaS CALL t3SetValue CALL MarkTDCursor(trace,dterm) IF lamflag=1 then CALL t4SetValues CALL w3Redraw END SUB SUB SliderResetAll(a,b,c,d) IF companion=0 then CALL h1mark(a) CALL h2mark(b) END IF CALL h3mark(c) CALL h4mark(d) END SUB SUB MatToThetaS LET s = (mb-mc) / 2 IF b+c<0 then ! and pi/40 then LET ang= atn((d-a)/(b+c)) / 2 IF a-d>=0 then LET theta= ang ELSE IF a-d<0 then LET theta= pi + ang END IF ELSE IF b+c=0 then IF a-d>0 then LET theta= pi/4 ELSE IF a-d<0 then LET theta= 3*pi/4 ELSE IF a-d=0 then LET mc= -mb LET md= a END IF END IF LET ma,a,oldma= round(ma,2) LET mb,b,oldmb= round(mb,2) LET mc,c,oldmc= round(mc,2) LET md,d,oldmd= round(md,2) END SUB ! SUB TraceDetToMat ! IF companion=true then ! LET ma,a,oldma= 0 ! LET mb,b,oldmb= 1 ! LET mc,c,oldmc= -dterm ! LET md,d,oldmd= trace ! ! ! CALL Classification(trace,dterm,dscrm,className$,classClr) ! ELSE ! ! Create the A matrix from trace,dterm and theta,s ! ! LET dscrm= round(trace^2 - 4*dterm,8) ! LET t2 = trace/2 ! LET td = trace^2/4 - dterm ! ! IF dscrm<0 then ! in parabola ! LET omega,s = sqr(abs(td))+.00001 ! ELSE ! LET omega,s = 0 ! END IF ! LET theta= 0 ! ! LET L = sqr(s^2 + td) ! LET Lcs= L*cos(2*theta) ! LET Lsn= L*sin(2*theta) ! ! LET ma = t2 - Lsn ! define matrix elements ! LET mb = Lcs + s ! LET mc = Lcs - s ! LET md = t2 + Lsn ! ! ! CALL TestResults ! ! LET cs= round(cos(theta),5) ! LET sn= round(sin(theta),5) ! LET omegaIV= sqr(abs(dterm - trace^2/4)) ! ! !IF dscrm=0 then CALL Classification(trace,dterm,dscrm,className$,classClr) ! ! LET ma,a,oldma= round(ma,2) ! LET mb,b,oldmb= round(mb,2) ! LET mc,c,oldmc= round(mc,2) ! LET md,d,oldmd= round(md,2) ! END IF ! ! CALL SliderResetAll(ma,mb,mc,md) ! CALL ShowMatrixValues ! END SUB ! ---- end of slider routines ---- SUB ResetDisplay CALL Classification(trace,dterm,dscrm,className$,classClr) CALL ShowMatrixValues CALL MarkTDCursor(trace,dterm) CALL t3SetValue IF lamflag=1 then CALL t4SetValues CALL w3Redraw END SUB ! --- eigen directions --- SUB FindLamdas ! eigen values LET dscrm= fDsc(trace,dterm) ! find discriminant IF dscrm>=0 then LET sqrtDscrm= sqr(dscrm) SELECT CASE dscrm CASE 0 ! one real root ! Case 0: dscrm=0 on parabola - one root - lam1=lam2 LET lam1,lam2= trace/2 CALL SetStrings(using$("-%.##",lam1),using$("-%.##",lam2),l1$,l2$) CASE is > 0 ! two real roots ! Case 1: dscrm>0 outside parabola - two roots - lam1=w2Lft and wx<=w2Rgt then CALL PlotDiamondRimClr(wx,w2y0,classClr) END IF IF lam2<>lam1 then LET wx= w2Wndx(lam2) IF wx>=w2Lft and wx<=w2Rgt then CALL PlotDiamondRimClr(wx,w2y0,classClr) END IF END IF ELSE LET wx = w2Wndx( lam1) LET wy1= w2Wndy( lam2) LET wy2= w2Wndy(-lam2) IF wy1>=w2Top and wy2<=w2Bas and wx>=w2Lft and wx<=w2Rgt then CALL PlotDiamondRimClr(wx,wy1,classClr) CALL PlotDiamondRimClr(wx,wy2,classClr) ! CALL DrawDiamond5(wx,wy1) ! CALL DrawDiamond5(wx,wy2) END IF END IF END SUB SUB FindEigenVectors(ma,mb,mc,md) LET dscrm= trace*trace - 4*dterm ! find discriminant IF dscrm>=0 then LET sqrtDscrm= sqr(dscrm) LET rma= ma ! round the matrix LET rmb= mb LET rmc= mc LET rmd= md ! LET places= 5 ! LET rma= round(ma,places) ! round the matrix ! LET rmb= round(mb,places) ! LET rmc= round(mc,places) ! LET rmd= round(md,places) IF dscrm=0 then ! on parabola - one real root IF rma=rmd then IF rmc=0 and rmb=0 then ! star ! no eigenlines ELSE IF rmb<>0 and rmc=0 then ! degenerate node CALL CheckEigenPoint(1,0) ! horizontal ELSE IF rmb=0 and rmc<>0 then ! degenerate node CALL CheckEigenPoint(0,1) ! vertical END IF ELSE IF rmc=0 and rmd=0 then ! first row points needed LET x= 2*rmb LET y= rmd-rma-sqrtDscrm CALL CheckEigenPoint(x,y) ELSE IF rmc=0 and rma=1 and rmb=1 and rmd=1 then CALL CheckEigenPoint(1,0) ! horizontal ELSE ! otherwise second row points LET x= rma-rmd-sqrtDscrm LET y= 2*rmc CALL CheckEigenPoint(x,y) END IF ELSE IF dscrm>0 then ! two real roots IF rmc=0 and rmd=0 then ! first row points needed? LET x= 2*rmb LET y= rmd-rma-sqrtDscrm CALL CheckEigenPoint(x,y) LET y= rmd-rma+sqrtDscrm CALL CheckEigenPoint(x,y) ELSE ! otherwise second row points IF rmc<>0 then LET y= 2*rmc LET x= rma-rmd-sqrtDscrm CALL CheckEigenPoint(x,y) LET x= rma-rmd+sqrtDscrm CALL CheckEigenPoint(x,y) ELSE CALL CheckEigenPoint(1,0) ! horizontal IF rmb=0 then CALL CheckEigenPoint(0,1) ! vertical ELSE LET x= 1 LET y= (rmd-rma)/rmb CALL CheckEigenPoint(x,y) END IF END IF END IF END IF CALL w3KeepGraphLayer END SUB SUB CheckEigenPoint(x,y) CALL SetVars(round(x,8),round(y,8),ex,ey) IF ex=0 and ey=1 then ! vertical CALL PlotLine(w3x0,w3Top+1,w3x0,w3Bas-1,classClr) ELSE IF ex=1 and ey=0 then ! horizontal CALL PlotLine(w3Lft+1,w3y0,w3Rgt-1,w3y0,classClr) ELSE IF abs(ey)>=abs(ex) then ! tall - scale y to hgt IF ey<>0 then LET scl= w3Size/ey CALL DrawEDirection(ex,ey,scl) ELSE CALL PlotLine(w3x0,w3Top+1,w3x0,w3Bas-1,classClr) END IF ELSE IF abs(ex)>abs(ey) then ! wide - scale x to wid IF ex<>0 then LET scl= w3Size/ex CALL DrawEDirection(ex,ey,scl) ELSE CALL PlotLine(w3Lft+1,w3y0,w3Rgt-1,w3y0,classClr) END IF END IF END SUB SUB DrawEDirection(ex,ey,scl) ! print ex,ey CALL SetVars(ex*scl,ey*scl,bx,by) CALL w3MathToPixels( bx, by,wx1,wy1) CALL w3MathToPixels(-bx,-by,wx2,wy2) ! PRINT wx1,wy1; wx2,wy2 CALL PlotLine(wx1,wy1,wx2,wy2,classClr) END SUB ! --- end of slider matrix routines --- ! (left branch) -- maybe "repeated eigenvalues (source)". same thing for ! "zero eigenvalue (source)" (right branch), (sink) for left branch. SUB Classification(trace,dterm,dscrm,className$,classClr) LET dscrm= trace*trace - 4*dterm ! find discriminant LET oldClassName$= className$ IF dterm=0 then ! horizontal axis of TD plane LET classClr = pink LET className$= "degenerate" ! trace= 0: points ! "origin" , "isolated fixed points" , "zero eigenvalues" ! trace<>0: parallel lines ! "nonisolated fixed points" , "multiple equilibrium points" , "zero eigenvalue" ELSE IF dterm<0 then ! saddle LET classClr = blue LET className$= "saddle" ELSE IF dterm>0 then ! upper half plane IF dscrm=0 then ! discriminate=0... on the parabola LET classClr= yellow IF ma=md and mb=0 and mc=0 then ! star IF trace>0 then LET className$= "star source" ELSE IF trace<0 then LET className$= "star sink" END IF ELSE ! node ! LET className$= "improper node" LET className$= "defective node" END IF ! "tr^[2] = 4 det" , "repeated eigenvalue" ELSE IF trace<0 then ! left half - stable IF dscrm>0 then ! complex - within parab LET classClr = green LET className$= "nodal sink" ! "stable node" ELSE LET classClr = red LET className$= "spiral sink" ! "stable spiral" END IF ELSE IF trace>0 then ! right half - unstable IF dscrm>0 then ! complex LET classClr = magenta LET className$= "nodal source" ! "unstable node" ELSE LET classClr = cyan LET className$= "spiral source" ! "unstable spiral" END IF ELSE IF trace=0 then ! neutral center LET classClr = white LET className$= "center" END IF END IF IF oldClassName$<>className$ then CALL ClassText(className$) END SUB ! ----- SUB MarkTDCursor(trace,dterm) CALL w1MathToPixels(trace,dterm,wt,wd) IF w1drawState=0 then CALL w1ShowGridLayer ELSE CALL w1ShowGraphLayer END IF IF w1drawState=1 then IF w1wWithin(wt,wd)=true then SET COLOR black PLOT wt-2,wd PLOT wt-1,wd-1; wt-1,wd+1 PLOT wt,wd-2; wt,wd+2 PLOT wt+1,wd-1; wt+1,wd+1 PLOT wt+2,wd ! PLOT wt,wd ! PLOT wt-1,wd; wt+1,wd ! PLOT wt,wd-1; wt,wd+1 CALL w1KeepGraphLayer END IF END IF IF w1wWithin(wt,wd)=true then SET COLOR backclr BOX AREA wt-3,wt+3,wd+1,wd-1 BOX AREA wt-1,wt+1,wd+3,wd-3 SET COLOR white PLOT wt-2,wd; wt+2,wd PLOT wt,wd-2; wt,wd+2 END IF END SUB SUB DrawDiamond5(wdx,wdy) PLOT wdx-2,wdy ; wdx+2,wdy PLOT wdx ,wdy-2; wdx ,wdy+2 BOX LINES wdx-1,wdx+1,wdy+1,wdy-1 END SUB ! SUB FindVector(x,y) ! LET wx= w3Wndx(x) ! LET wy= w3Wndy(y) ! LET dx= dxdt(x,y) ! LET dy= dydt(x,y) ! ! CALL RungeKutta4(x,y,dx,dy,tstep) ! IF dx<>0 or dy<>0 then ! LET ang= angle(dx,-dy) ! DRAW BigVector with rotate(ang) * shift(wx,wy) ! END IF ! END SUB ! PICTURE BigVector ! PLOT -6, 0; 6,0 ! PLOT 3,-3; 6,0; 3, 3 ! END PICTURE ! PICTURE ICPointer ! PLOT -3, 3; 0,0; -3,-3 ! END PICTURE ! ! PICTURE Vector ! PLOT -4, 0; 4,0 ! PLOT 2,-2; 4,0; 2, 2 ! END PICTURE ! ! PICTURE Pointer ! PLOT 0,0; -4,-4 ! PLOT 0,0; -4, 4 ! END PICTURE ! ---- Vector Field routines ---- ! --------- direction and vector fields ------------ SUB DirectionField CALL w3ShowGridLayer SET COLOR drkmid LET vstp= 30 LET hstp= vstp/2 FOR wx= w3Lft+hstp to w3Rgt-hstp step vstp ! draw vector field FOR wy= w3Top+hstp to w3Bas-hstp step vstp LET x = w3fncx(wx) LET y = w3fncy(wy) LET dx= dxdt(x,y) LET dy= dydt(x,y) IF dx<>0 or dy<>0 then LET ang= angle(dx,-dy) LET sa = 5*sin(ang) LET ca = 5*cos(ang) LET wx0= wx+ca LET wy0= wy+sa LET wx2= wx-ca LET wy2= wy-sa PLOT wx2,wy2; wx0,wy0 CALL VectorHead(ang,wx2,wy2,wx0,wy0) ELSE PLOT wx,wy END IF NEXT wy NEXT wx CALL w3KeepFieldLayer END SUB SUB VectorHead(ang,wx,wy,wx0,wy0) CALL SetMat(ang) CALL Rotate(6, 2,x1,y1) CALL Rotate(6,-2,x2,y2) LET wx1= wx+x1 LET wy1= wy+y1 LET wx2= wx+x2 LET wy2= wy+y2 PLOT wx1,wy1; wx0,wy0; wx2,wy2 END SUB SUB Rotate(xin,yin,xout,yout) LET xout= xin*vma + yin*vmc LET yout= xin*vmb + yin*vmd END SUB SUB SetMat(ang) LET vma= cos(ang) LET vmb= sin(ang) LET vmc= -vmb LET vmd= vma END SUB ! SUB VectorField ! SET COLOR drkmid ! LET vstp= 30 ! LET hstp= vstp/2 ! FOR wx= w3Lft+hstp to w3Rgt-hstp step vstp ! draw vector field ! FOR wy= w3Top+hstp to w3Bas-hstp step vstp ! LET x = w3fncx(wx) ! LET y = w3fncy(wy) ! LET dx= dxdt(x,y) ! LET dy= dydt(x,y) ! IF dx<>0 or dy<>0 then ! LET ang= angle(dx,-dy) ! LET sa = 10*sin(ang) ! LET ca = 10*cos(ang) ! LET wx0= wx+ca ! LET wy0= wy+sa ! PLOT wx,wy; wx0,wy0 ! CALL VectorHead(ang,wx,wy,wx0,wy0) ! ELSE ! PLOT wx,wy ! END IF ! NEXT wy ! NEXT wx ! END SUB ! ---- end of vector field routines --- ! ---- trajectory and RK4 ----- SUB Trajectory(x0,y0) SET COLOR classClr PLOT LET oldstep= tstep IF dscrm<0 then LET tstep= 1/256 FOR n= 1 to -1 step -2 LET oldFlag= 1 LET stp = n*tstep CALL SetVars(x0,y0,x,y) CALL w3MathToPixels(x,y,wx,wy) FOR i= 1 to 40/tstep ! draw trajectory CALL SetVars(wx,wy,oldwx,oldwy) CALL w3MathToPixels(x,y,wx,wy) IF w3wWithin(wx,wy)=true then IF oldFlag=1 then PLOT wx,wy; oldwx,oldwy LET oldFlag= 1 ELSE LET oldFlag= 0 END IF CALL RungeKutta4(x,y,dx,dy,stp) ! equilibrium or off plane? IF (abs(dx)<.02 and abs(dy)<.02) or abs(x)>16 or abs(y)>16 then EXIT FOR END IF IF trace=0 and dterm>0 and i>100 and abs(x-x0)<.05 and abs(y-y0)<.05 then PLOT wx,wy; w3Wndx(x0),w3Wndy(y0) EXIT FOR END IF NEXT i PLOT IF trace=0 and dterm>0 then EXIT FOR NEXT n LET tstep= oldstep END SUB ! ---- Runge Kutta ---- SUB RungeKutta4(x,y,dx,dy,tstp) LET hstp= .5*tstp LET dx1= dxdt(x,y) LET dy1= dydt(x,y) LET x1 = x + dx1*hstp LET y1 = y + dy1*hstp LET t1 = t + hstp LET dx2= dxdt(x1,y1) LET dy2= dydt(x1,y1) LET x2 = x1 + dx2*hstp LET y2 = y1 + dy2*hstp LET t2 = t + hstp LET dx3= dxdt(x2,y2) LET dy3= dydt(x2,y2) LET x3 = x2 + dx3*tstp LET y3 = y2 + dy3*tstp LET dx4= dxdt(x3,y3) LET dy4= dydt(x3,y3) LET dx = (dx1 + 2*dx2 + 2*dx3 + dx4) / 6 LET dy = (dy1 + 2*dy2 + 2*dy3 + dy4) / 6 LET x = x + tstp*dx LET y = y + tstp*dy END SUB SUB DEcoords(inx,iny,tmstp,dx,dy,outx,outy ) CALL DEchange(inx,iny,dx,dy) LET outx= inx + dx*tmstp LET outy= iny + dy*tmstp END SUB SUB DEchange(inx,iny,dx,dy) LET dx= dxdt(inx,iny) LET dy= dydt(inx,iny) END SUB ! --- TD Plane --- ! ---- Draw w1 Trace Determinant plane ---- SUB DrawTDPlane CALL BoxArea(w1Lft,w1Rgt,w1Bas-1,w1y0,blue) CALL BoxArea(w1Lft,w1x0-1,w1y0-1,w1Top,green) CALL BoxArea(w1x0+1,w1Rgt,w1y0,w1Top,magenta) CALL PlotLine(w1Lft,w1y0 ,w1Rgt,w1y0 ,pink) CALL PlotLine(w1Lft,w1y0+1,w1Rgt,w1y0+1,pink) FOR wx= w1Lft to w1Rgt ! fill inside parabola LET fx= w1Fncx(wx) LET fy= 0.25*fx*fx LET wy= w1Wndy(fy) IF fx<0 then CALL PlotLine(wx,wy, wx,w1Top, red) ELSE CALL PlotLine(wx,wy, wx,w1Top, cyan) END IF NEXT wx CALL SetLineWeight(2) SET COLOR yellow ! outline parabola IF m68kFlag=1 or unixFlag=1 then FOR wx= w1Lft to w1Rgt LET fx= w1Fncx(wx) LET fy= 0.25*fx*fx LET wy= w1Wndy(fy) IF round(wy)=w1y0 then LET wy= wy-1 IF fx<0 then PLOT wx,wy; ELSE ! shift for symmetry PLOT wx-1,wy; END IF NEXT wx ELSE FOR wx= w1Lft to w1Rgt LET fx= w1Fncx(wx) LET fy= 0.25*fx*fx LET wy= w1Wndy(fy)+1 PLOT wx+1,wy; NEXT wx END IF PLOT CALL SetLineWeight(1) CALL PlotLine(w1x0,w1Top,w1x0,w1y0,white) END SUB END SUB