!  Taylor Poly's: There are mistakes in the "terms" here. For example,
!  the terms for the second option, the polynomial f(x), are clearly
!  wrong. There's no hurry on this one, but I will try to look at it over
!  the weekend.


!  The Scottish mathematician James Gregorie, in 1675, found that 
!  a function f that can be expressed as a power series centered at a may be written
!  f(x) = c0 + c1(x-a) + c2(x-a)^2 + c3(x-a)^3 + ...
!  and using term-by-term differentiation the constants cn can be expressed as
!  cn = f^n(a)/n!
!  so the function can be written
!  f(x) = f(a) + f'(a)(x-a) + f''(a)/2! (x-a)^2 + f'''(a)/3! (x-a)^3 + f''''(a)/4! (x-a)^4 + ...
!  Gregorie's method of representing functions was published by Brooks Taylor in 1715, and has come to be called the Taylor Series. When the series is centered at 0 (a=0), it is known as a Maclaurin Series.
!  f(x) = f(0) + f'(0)x + f''(0)/2! x^2 + f'''(0)/3! x^3 + f''''(0)/4! x^4 + ...


!  Taylor Polynomials
!  
!  A function is selected from a menu at lower left. Its graph is drawn
!  in blue on the graphing window. Below it a slider selects a value a
!  of  x.  The value of  a  can be increased or decreased using the 
!  >  and  <  keys, or animated using the  >>>  key.
!  At bottom center a slider selects  n  between  0  and  9.
!  The Taylor polynomial of this order at the point  a  is displaayed
!  in yellow. 
!  
!  A toggle at bottom right allows display of the Taylor polynomial formulae,
!  and a readout of the polynomial at the selected value of a.


!! File: TaylorPolynomials
!! 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$ = "Taylor Polynomials"

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

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


! ***

SUB NthTaylor
    DECLARE PUBLIC worklft,workrgt,workbas,worktop,workmid      ! work area
    DECLARE PUBLIC black,drkgry,drkmid,midgry,litmid,litgry,white
    DECLARE PUBLIC red,yellow,green,cyan,blue,magenta,colorscheme
    DECLARE PUBLIC slideclr,axislabelclr,true,false

    DECLARE DEF InfoWithin, QuitWithin

    ! --- help page text ---

    DIM info$(1:16)
    LET excl$= chr$(33)
    MAT READ info$

    DATA "Taylor Polynomials"
    DATA ""
    DATA "This tool supports graphical exploration of functions that can be expressed as a power series."
    DATA ""
    DATA "Click the [v] popup menu button in the lower left to see the list of equations."
    DATA "Click the mouse on a function to select it and draw its graph in cyan."
    DATA ""
    DATA "Click or drag the n slider to select the order of the Taylor polynomial series."
    DATA "The Taylor polynomial of order  n  centered at the point  a  is displayed in yellow."
    DATA ""
    DATA "Click or drag the  a  slider to define the center of convergence for the Taylor series."
    DATA "Click the [<] or [>] button to step the value of a."
    DATA "Click the [>>>] button to animate the value of a."
    DATA ""
    DATA "Click the [Formulas] button to display the initial terms of the Taylor polynomial."
    DATA "Click the [Terms] button to see the first n terms of the Taylor polynomial with their coefficients."

    !A function is selected from a menu at lower left. Its graph is drawn in blue on the graphing window. 
    !Below it a slider selects a value a of  x.  
    !The value of  a  can be increased or decreased using the >  and  <  keys, or animated using the  >>>  key. 
    !At bottom center a slider selects  n  between  0  and  9. 
    !The Taylor polynomial of this order at the point  a  is displaayed in yellow.
    !
    !A toggle at bottom right allows display of the Taylor polynomial formulae, and a readout of the polynomial at the selected value of a.

    !    LET help$( 2) = ""
    !    LET help$( 3) = "n" & excl$ & " or n factorial"
    !    LET help$( 4) = "By definition"
    !    LET help$( 5) = "n!( n+1)" & excl$ & " = ( n+1)!"
    !    LET help$( 6) = "Then"
    !    LET help$( 7) = "0" & excl$ & "( 0+1)" & excl$ & " = ( 0+1)" & excl$
    !    LET help$( 8) = "0" & excl$ & "1" & excl$ & " = 1" & excl$
    !    LET help$( 9) = "Therefore"
    !    LET help$(10) = "0" & excl$ & " must equal 1"


    ! ------------- Functions and equations -------------

    PUBLIC FncList$(1:4), DrvList$(1:4)          ! eqmarks(1:4)

    DECLARE DEF roundn,clamp,e

    SHARE fact(1:10)
    LET fact[1] = 1
    FOR i= 2 TO 10
        LET fact[i] = i*fact[i-1]                ! fact[n] = n!
    NEXT i

    DEF asin(x)
        IF abs(x)<1 then
           LET as = atn(x/sqr(1-x*x))
        ELSEIF x= 1 then
           LET as = pi/2
        ELSEIF x=-1 then
           LET as = -pi/2
        ELSE
           LET as = badarg
           ! CAUSE EXCEPTION -3000, "Argument not in range"
        END IF
        LET asin= as
    END DEF

    !DEF e= exp(1)

    LET badarg= 987656789                        ! .0123456789
    LET bignum= 10^10
    LET eqcnt = 17


    MAT redim FncList$(1:eqcnt)
    MAT READ FncList$
    DATA "f(x) = 1/sqrt(1+x)"
    DATA "f(x) = 0.3x^[2]-.01x^[5]"
    DATA "f(x) = sin(x)"
    DATA "f(x) = e^[x]"
    DATA "f(x) = 1/(2x)"
    DATA "f(x) = cos(x)"
    DATA "f(x) = sin(x^[2])"
    DATA "f(x) = tan(x)"
    DATA "f(x) = sec(x)"
    DATA "f(x) = -2x/(1+x^[2])^[2]"
    DATA "f(x) = 1/(1+x^[2])"
    DATA "f(x) = tan^[-1](x)"
    DATA "f(x) = x tan^[-1](x)-ln(1+x^[2])/2"
    DATA "f(x) = sin^[-1](x)"
    DATA "f(x) = log(x)"
    DATA "f(x) = cos(x) + cos(9x)/9"
    DATA "f(x) = exp(-1/x^[2])    [x>=0]"

    SUB FixEqnum(eq)
        SELECT CASE eq
        CASE 1
             LET eqnum= 1
        CASE 2
             LET eqnum= 2
        CASE 3
             LET eqnum= 3
        CASE 4
             LET eqnum= 4
        CASE 5
             LET eqnum= 5
        CASE 6
             LET eqnum= 6
        CASE 7
             LET eqnum= 7
        CASE 8
             LET eqnum= 8
        CASE 9
             LET eqnum= 9
        CASE 10
             LET eqnum= 10
        CASE 11
             LET eqnum= 11
        CASE 12
             LET eqnum= 12
        CASE 13
             LET eqnum= 13
        CASE 14
             LET eqnum= 14
        CASE 15
             LET eqnum= 15
        CASE 16
             LET eqnum= 16
        CASE 17
             LET eqnum= 17
        END SELECT
    END SUB

    SUB GetFnc(eqnum,x,y)
        SELECT CASE eqnum
        CASE 1
             IF x>-1 then
                LET y= 1/sqr(1+x)
             ELSE IF x=-1 then
                LET y= bignum
             ELSE
                LET y= badarg
             END IF
        CASE 2
             LET y= 0.3*x*x-.01*x^5
        CASE 3
             LET y= sin(x)
        CASE 4
             LET y= e^(x)
        CASE 5
             IF x<>0 then
                LET y= 1/(2*x)
             ELSE
                LET y= bignum
             END IF
        CASE 6
             LET y = cos(x)
        CASE 7
             LET y = sin(x^2)
        CASE 8
             IF (x>-1.58 and x<-1.56) or (x>4.72 and x<4.74) THEN
                LET y = bignum
             ELSE IF 1.56<x AND x<1.58 THEN
                LET y = bignum
             ELSE
                LET y = tan(x)
             END IF
        CASE 9
             IF (x>-1.58 and x<-1.56) or (x>4.72 and x<4.74) THEN
                LET y = bignum
             ELSE IF 1.56 < x AND x < 1.58 THEN
                LET y = bignum
             ELSE
                LET y = 1/cos(x)
             END IF
        CASE 10
             LET y = -2*x/(1+x^2)^2
        CASE 11
             LET y = 1/(1+x^2)
        CASE 12
             LET y = atn(x)
        CASE 13
             LET y = x*atn(x) - log(1+x^2)/2
        CASE 14
             LET y = asin(x)
        CASE 15
             IF x>0 THEN
                LET y = log(x)
             ELSE
                LET y = badarg
             END IF
        CASE 16
             LET y= cos(x) + cos(9*x)/9
        CASE 17
             IF x>0 then
                LET y= e^(-1/x^2)
             ELSE
                LET y= 0
             END IF
        END SELECT
    END SUB

    SUB GetNine(eqnum,xa,ya,d1,d2,d3,d4,d5,d6,d7,d8,d9)
        LOCAL temp,C,T,Th,F,Fi,S,Se,E,Ni,ten,denom,A,x,y,B,D
        LET x= round(xa,8)
        LET y= ya
        SELECT CASE eqnum
        CASE 1                                   ! y = 1/sqr(1+x)
             IF x>-1 then
                LET denom=   2*(1+x)
                LET d1=     -y/denom
                LET d2=  -3*d1/denom
                LET d3=  -5*d2/denom
                LET d4=  -7*d3/denom
                LET d5=  -9*d4/denom
                LET d6= -11*d5/denom
                LET d7= -13*d6/denom
                LET d8= -15*d7/denom
                LET d9= -17*d8/denom
             ELSE
                LET y= bignum
             END IF
        CASE 2
             !LET y =  0.3*x*x-.01*x^5
             LET d1=  0.6*x - 0.05*x^4
             LET d2=  0.6 - 0.2*x^3
             LET d3= -0.6*x^2
             LET d4= -1.2*x
             LET d5= -1.2
             LET d6=  0
             LET d7=  0
             LET d8=  0
             LET d9=  0
        CASE 3
             !LET y = sin(x)
             LET d1= cos(x)
             LET d2= -y
             LET d3= -d1
             LET d4=  y
             LET d5=  d1
             LET d6= -y
             LET d7= -d1
             LET d8=  y
             LET d9=  d1
        CASE 4
             !LET y = exp(x)
             LET d1= y
             LET d2= y
             LET d3= y
             LET d4= y
             LET d5= y
             LET d6= y
             LET d7= y
             LET d8= y
             LET d9= y
             !        CASE 4
             !             !LET y = exp(-x^2)
             !             LET d1= -y
             !             LET d2=  y
             !             LET d3= -y
             !             LET d4=  y
             !             LET d5= -y
             !             LET d6=  y
             !             LET d7= -y
             !             LET d8=  y
             !             LET d9= -y
             !        CASE 4
             !             !LET y = exp(-x)
             !             LET d1= -y
             !             LET d2=  y
             !             LET d3= -y
             !             LET d4=  y
             !             LET d5= -y
             !             LET d6=  y
             !             LET d7= -y
             !             LET d8=  y
             !             LET d9= -y
        CASE 5
             IF x<>0 then
                !LET y = 1/(2x)
                LET d1= -y/x
                LET d2= -2*d1/x
                LET d3= -3*d2/x
                LET d4= -4*d3/x
                LET d5= -5*d4/x
                LET d6= -6*d5/x
                LET d7= -7*d6/x
                LET d8= -8*d7/x
                LET d9= -9*d8/x
             ELSE
                !LET y= bignum
             END IF
        CASE 6
             !LET y =  cos(x)
             LET d1= -sin(x)
             LET d2= -y
             LET d3= -d1
             LET d4=  y
             LET d5=  d1
             LET d6= -y
             LET d7= -d1
             LET d8=  y
             LET d9=  d1
        CASE 7
             LET T = x*x
             !LET y = sin(T)
             LET C = cos(T)
             LET Th= T*x
             LET F = Th*x
             LET Fi= F*x
             LET S = Fi*x
             LET Se= S*x
             LET E = Se*x
             LET Ni= E*x
             LET d1= 2*C*x
             LET d2=   -4*y*T + 2*C
             LET d3=   -8*C*Th - 12*y*x
             LET d4=   (16*F - 12)*y - 48*C*T
             LET d5=   (32*Fi - 120*x)*C + 160*y*Th
             LET d6=  (-64*S + 720*T)*y + (480*F - 120)*C
             LET d7= (-128*Se + 3360*Th)*C + (-1344*Fi + 1680*x)*y
             LET d8= (256*E -13440*F + 1680)*y + (-3584*S + 13440*T)*C
             LET d9=  (512*Ni - 48384*Fi +30240*x)*C+ (9216*Se  - 80640*Th)*y
        CASE 8
             LET A = sin(x)
             LET C = cos(x)
             LET T = C*C
             !LET y = tan(T)
             LET Th= T*C
             LET F = Th*C
             LET Fi= F*C
             LET S = Fi*C
             LET Se= S*C
             LET E = Se*C
             LET Ni= E*C
             LET ten = Ni*C
             LET d1= 1/T
             LET d2=   2*A/Th
             LET d3=   -2*(2*T - 3)/F
             LET d4=   -8*A*(T - 3)/Fi
             LET d5=   8*(2*F - 15*T + 15)/S
             LET d6=  16*A*(2*F - 30*t + 45)/Se
             LET d7= -16*(4*S -126*F +420*T -315)/E
             LET d8=  -128*A*(S - 63*F + 315*T -315)/Ni
             LET d9=  128*(2*E -255*S + 2205*F -4725*T + 2835)/ten
        CASE 9
             LET A = sin(x)
             LET C = cos(x)
             LET T = C*C
             !LET y = secant(x)
             LET Th= T*C
             LET F = Th*C
             LET Fi= F*C
             LET S = Fi*C
             LET Se= S*C
             LET E = Se*C
             LET Ni= E*C
             LET ten = Ni*C
             LET d1= A/T
             LET d2= (-T + 2)/Th
             LET d3= -A*(T - 6)/F
             LET d4= (F - 20*T + 24)/Fi
             LET d5= A*(F - 60*T + 120)/S
             LET d6= (-S + 182*F - 840*T + 720)/Se
             LET d7= -A*(S - 546*F + 4200*T - 5040)/E
             LET d8= (E - 1640*S + 23184*F - 60480*T + 40320)/Ni
             LET d9= A*(E - 4920*S + 115920*F - 423360*T + 362880)/ten
        CASE 10
             LET T = x*x
             LET F = T*x*x
             LET S = F*x*x
             LET E = S*x*x
             LET ten = E*x*x
             LET temp = 1+T
             LET temp = temp*(1+T)
             !LET y = -2*x/temp
             LET temp = temp*(1+T)
             LET d1 = 2*(3*T - 1)/temp
             LET temp = temp*(1+T)
             LET d2 = -24*x*(T - 1)/temp
             LET temp = temp*(1+T)
             LET d3 = 24*(5*F - 10*T +1)/temp
             LET temp = temp*(1+T)
             LET d4 = -240*x*(3*F - 10*T +3)/temp
             LET temp = temp*(1+T)
             LET d5 = 720*(7*S - 35*F +21*T -1)/temp
             LET temp = temp*(1+T)
             LET d6 = -40320*x*(S -7*F + 7*T -1)/temp
             LET temp = temp*(1+T)
             LET d7 = 40320*(9*E -84*S +126*F -36*T +1)/temp
             LET temp = temp*(1+T)
             LET d8 = -725760*x*(5*E -60*S + 126*F -60*T +5)/temp
             LET temp = temp*(1+T)
             LET d9 = 3628800*(11*ten - 165*E +462*S -330*F +55*T-1)/temp
        CASE 11
             LET T = x*x
             LET F = T*x*x
             LET S = F*x*x
             LET E = S*x*x
             LET temp = 1+T
             !LET y = 1/temp
             LET temp = temp*(1+T)
             LET d1 = -2*x/temp
             LET temp = temp*(1+T)
             LET d2 = 2*(3*T - 1)/temp
             LET temp = temp*(1+T)
             LET d3 = -24*x*(T - 1)/temp
             LET temp = temp*(1+T)
             LET d4 = 24*(5*F - 10*T +1)/temp
             LET temp = temp*(1+T)
             LET d5 = -240*x*(3*F - 10*T +3)/temp
             LET temp = temp*(1+T)
             LET d6 = 720*(7*S - 35*F +21*T -1)/temp
             LET temp = temp*(1+T)
             LET d7 = -40320*x*(S -7*F + 7*T -1)/temp
             LET temp = temp*(1+T)
             LET d8 = 40320*(9*E -84*S +126*F -36*T +1)/temp
             LET temp = temp*(1+T)
             LET d9 = -725760*x*(5*E -60*S + 126*F -60*T +5)/temp
        CASE 12
             LET T = x*x
             LET F = T*x*x
             LET S = F*x*x
             LET E = S*x*x
             !LET y = atn(x)
             LET temp = 1+T
             LET d1 = 1/temp
             LET temp = temp*(1+T)
             LET d2 = -2*x/temp
             LET temp = temp*(1+T)
             LET d3 = 2*(3*T - 1)/temp
             LET temp = temp*(1+T)
             LET d4 = -24*x*(T - 1)/temp
             LET temp = temp*(1+T)
             LET d5 = 24*(5*F - 10*T +1)/temp
             LET temp = temp*(1+T)
             LET d6 = -240*x*(3*F - 10*T +3)/temp
             LET temp = temp*(1+T)
             LET d7 = 720*(7*S - 35*F +21*T -1)/temp
             LET temp = temp*(1+T)
             LET d8 = -40320*x*(S -7*F + 7*T -1)/temp
             LET temp = temp*(1+T)
             LET d9 = 40320*(9*E -84*S +126*F -36*T +1)/temp
        CASE 13
             LET T = x*x
             LET F = T*x*x
             LET S = F*x*x
             LET E = S*x*x
             LET temp = 1+T
             LET A = atn(x)
             !LET y = x*A - log(temp)/2
             LET d1 = A
             LET d2 = 1/temp
             LET temp = temp*(1+T)
             LET d3 = -2*x/temp
             LET temp = temp*(1+T)
             LET d4 = 2*(3*T - 1)/temp
             LET temp = temp*(1+T)
             LET d5 = -24*x*(T - 1)/temp
             LET temp = temp*(1+T)
             LET d6 = 24*(5*F - 10*T +1)/temp
             LET temp = temp*(1+T)
             LET d7 = -240*x*(3*F - 10*T +3)/temp
             LET temp = temp*(1+T)
             LET d8 = 720*(7*S - 35*F +21*T -1)/temp
             LET temp = temp*(1+T)
             LET d9 = -40320*x*(S -7*F + 7*T -1)/temp
        CASE 14
             IF y<>badarg THEN
                LET T = x*x
                LET F = T*x*x
                LET S = F*x*x
                LET E = S*x*x
                !LET y = asin(x)
                IF T>1 then                      ! n<0 to a fractional power?
                   LET y=badarg
                ELSE
                   LET temp = (1-T)^.5
                   IF temp=0 then
                      LET T= .999
                      LET temp= (1-T)^.5
                   END IF
                   LET d1 = 1/temp
                   LET temp = temp*(1-T)
                   LET d2 = x/temp
                   LET temp = temp*(1-T)
                   LET d3 = (2*T + 1)/temp
                   LET temp = temp*(1-T)
                   LET d4 = 3*x*(2*T + 3)/temp
                   LET temp = temp*(1-T)
                   LET d5 = 3*(8*F + 24*T +3)/temp
                   LET temp = temp*(1-T)
                   LET d6 = 15*x*(8*F + 40*T +15)/temp
                   LET temp = temp*(1-T)
                   LET d7 = 45*(16*S +120*F +90*T +5)/temp
                   LET temp = temp*(1-T)
                   LET d8 = 315*x*(16*S +168*F + 210*T +35)/temp
                   LET temp = temp*(1-T)
                   LET d9 = 315*(128*E +1792*S + 3360*F +1120*T +35)/temp
                END IF
             END IF
        CASE 15
             ! LET y = log(x)
             IF x<>0 then
                LET d1=     1/x
                LET d2=   -d1/x
                LET d3= -2*d2/x
                LET d4= -3*d3/x
                LET d5= -4*d4/x
                LET d6= -5*d5/x
                LET d7= -6*d6/x
                LET d8= -7*d7/x
                LET d9= -8*d8/x
             END IF
        CASE 16
             LET T = 9
             LET A = cos(x)
             LET C = cos(T*x)
             LET B = sin(x)
             LET D = sin(T*x)
             LET Th= T*T
             LET F = Th*T
             LET Fi= F*T
             LET S = Fi*T
             LET Se= S*T
             LET E = Se*T
             LET Ni= E*T
             LET d1= -B - D
             LET d2= -A - T*C
             LET d3=  B + Th*D
             LET d4=  A + F*C
             LET d5= -B - Fi*D
             LET d6= -A - S*C
             LET d7=  B + Se*D
             LET d8=  A + E*C
             LET d9= -B - Ni*D
        CASE 17
             IF x>0 THEN
                LET A = x*x
                LET C = exp(-1/A)
                LET T = A*A
                LET Th= T*A
                LET F = Th*A
                LET Fi= F*A
                LET S = Fi*A
                LET Se= S*A
                LET E = Se*A
                LET Ni= E*A
                LET ten = Ni*A
                LET d1= 2*C/(A*x)
                LET d2=  -2*C*(3*A - 2)/Th
                LET d3= 4*C*(6*T - 9*A + 2)/(F*x)
                LET d4= -4*C*(30*Th - 75*T + 36*A - 4)/S
                LET d5= 8*C*(90*F - 330*Th + 255*T - 60*A +4)/(Se*x)
                LET d6= -8*C*(630*Fi-3150*F+3465*Th-1290*T+180*A-8)/Ni
                LET d7= 16*C*(2520*S-16380*Fi+23940*F-12495*Th+2730*T-252*A+8)/(ten*x)
                LET d8= -16*C*(22680*Se-185220*S+343980*Fi-235305*F+71400*Th-10248*T+672*A-16)/(S*S)
                LET d9= 32*C*(113400*E-1134000*Se+2593080*S-2226420*Fi+877905*F-173880*Th+17640*T-864*A+16)/(ten*Th*x)
             ELSE
                LET d1=0
                LET d2=0
                LET d3=0
                LET d4=0
                LET d5=0
                LET d6=0
                LET d7=0
                LET d8=0
                LET d9=0
             END IF
        END SELECT
    END SUB

    DEF Taylor(n,x,xa,ya)
        LOCAL x1,x2,x3,x4,x5,x6,x7,x8,x9,ty
        LET x1= (x-xa)
        LET x2= x1/2
        LET x3= x1/3
        LET x4= x1/4
        LET x5= x1/5
        LET x6= x1/6
        LET x7= x1/7
        LET x8= x1/8
        LET x9= x1/9

        SELECT CASE n
        CASE 0
             LET ty= ya
        CASE 1
             LET ty= ya+x1*d1
        CASE 2
             LET ty= ya+x1*(d1+x2*d2)
        CASE 3
             LET ty= ya+x1*(d1+x2*(d2+x3*d3))
        CASE 4
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*d4)))
        CASE 5
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*(d4+x5*d5))))
        CASE 6
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*(d4+x5*(d5+x6*d6)))))
        CASE 7
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*(d4+x5*(d5+x6*(d6+x7*d7))))))
        CASE 8
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*(d4+x5*(d5+x6*(d6+x7*(d7+x8*d8)))))))
        CASE 9
             LET ty= ya+x1*(d1+x2*(d2+x3*(d3+x4*(d4+x5*(d5+x6*(d6+x7*(d7+x8*(d8+x9*d9))))))))
        CASE else
        END SELECT
        LET Taylor= ty
    END DEF


    ! ------------- Graphing planes  --------------------
    ! ------------- w1 Graphing plane ------------

    ! --- w1 Plane 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= 0
    LET w1yPiFlag= 0
    LET w1xMult  = 1
    LET w1yMult  = 1

    LET w1Top  = worktop + 85                    ! window values
    LET w1Bas  = w1Top + 320
    LET w1Lft  = worklft + 55
    LET w1Rgt  = w1Lft + 480

    LET w1fTop =  2                              ! function values
    LET w1fBas = -2
    LET w1fLft = -2
    LET w1fRgt =  4

    LET w1xAx$ = "x"                             ! 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 = 0
    LET w1xFirst = 0

    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 w1wWithin

    CALL w1Variables

    SUB w1Init
        BOX CLEAR w1Lft-20,w1Rgt+20,w1Top,w1Top-20
        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 w1KeepGridLayer
    END SUB

    SUB w1SetBounds(eqnum)
        LET w1xPiFlag= 0
        LET w1xMult  = 1
        SELECT CASE eqnum
        CASE 1
             LET w1fLft= -1
             LET w1fRgt=  7
             LET w1fBas=  0
             LET w1fTop=  4
        CASE 2,13
             LET w1fLft= -4
             LET w1fRgt=  4
             LET w1fBas= -1
             LET w1fTop=  3
        CASE 11
             LET w1fLft= -4
             LET w1fRgt=  4
             LET w1fBas= -0.5
             LET w1fTop=  1.5
        CASE 10
             LET w1fLft= -4
             LET w1fRgt=  4
             LET w1fBas= -1
             LET w1fTop=  1
        CASE 3,6,16
             LET w1fLft=  -2                     !*pi
             LET w1fRgt=   4                     !*pi
             LET w1fBas=  -2
             LET w1fTop=   2
             LET w1xPiFlag= 1
             LET w1xMult  = pi
        CASE 4
             LET w1fLft= -4
             LET w1fRgt=  2
             LET w1fBas=  0
             LET w1fTop=  4
        CASE 5,7,12
             LET w1fLft= -4
             LET w1fRgt=  4
             LET w1fBas= -2
             LET w1fTop=  2
        CASE 8,9
             LET w1fLft= -3                      !pi/2
             LET w1fRgt=  3                      !3*pi/2
             LET w1fBas= -5
             LET w1fTop=  5
             LET w1xPiFlag= 1
             LET w1xMult  = pi
        CASE 14
             LET w1fLft= -1
             LET w1fRgt=  1
             LET w1fBas= -2
             LET w1fTop=  2
        CASE 15
             LET w1fLft= -1
             LET w1fRgt=  7
             LET w1fBas= -4
             LET w1fTop=  2
        CASE 17
             LET w1fLft= -1
             LET w1fRgt=  3
             LET w1fBas= -0.1
             LET w1fTop=  1
        END SELECT

        SELECT CASE eqnum
        CASE 17
             LET w1xSTik =  0.5                  ! horizontal axis Tik marks
             LET w1xLTik =  1
             LET w1xFirst= w1fLft
             LET w1ySTik =  0.1                  ! vertical axis Tik marks
             LET w1yLTik =  0.5
             LET w1yLabel=  0.5

             LET w1yFirst=  0
        CASE 11
             LET w1xSTik =  0.5                  ! horizontal axis Tik marks
             LET w1xLTik =  1
             LET w1xFirst= w1fLft
             LET w1ySTik =  0.1                  ! vertical axis Tik marks
             LET w1yLTik =  0.5
             LET w1yLabel=  0.5

             LET w1yFirst= -0.5
        CASE else
             LET w1xSTik = 0.5                   ! horizontal axis Tik marks
             LET w1xLTik = 1
             LET w1xFirst= w1fLft
             LET w1ySTik = .5                    ! vertical axis Tik marks
             LET w1yLTik = 1
             LET w1yLabel= 1

             LET w1yFirst= w1fBas
        END SELECT
        LET w1xLabel= 0

        LET w1Rgt= w1Lft+480
        LET w1Bas= w1Top+320
        CALL w1Variables

        LET xa = 0.5
        LET wxa= w1Wndx(xa)
    END SUB

    ! ------------- Sliders --------------
    ! ----- h1 horizontal slider - sets degree n ---

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

    ! --- Slider vars ---

    LET h1PiAxis= 0
    LET h1Mult  = 1

    LET h1clr   = slideclr
    LET h1name$ = "n"
    LET h1form$ = "-%"
    LET h1Places= 0
    LET h1Click = 1

    LET h1axis  = workbas-22
    LET h1wLft  = w1Rgt-108
    LET h1wRgt  = h1wLft+108

    LET h1fLft  = 0
    LET h1fRgt  = 9

    LET h1STik  = 0                              ! short tick marks
    LET h1LTik  = 1                              ! long tick marks
    LET h1Label = 1                              ! labels
    LET h1First = h1fLft                         ! first tick mark

    DECLARE DEF h1Within                         !,h1fWithin  ! window/function transforms

    CALL h1SliderVariables

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

    ! --- h2 horizontal slider - sets t ---

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

    LET h2PiAxis= 1
    LET h2Mult  = pi

    LET h2clr   = slideclr
    LET h2name$ = "a"
    LET h2form$ = "--%.###"
    LET h2Places= 8
    LET h2Click = 0.5

    LET h2axis  = w1Bas + 21
    LET h2wLft  = w1Lft
    LET h2wRgt  = w1Rgt

    LET h2fLft  = -2
    LET h2fRgt  =  4

    LET h2STik  =  0.5                           ! 0.5                             ! short tick marks
    LET h2LTik  =  1                             ! long tick marks
    LET h2Label =  1                             ! labels
    LET h2First = -6                             ! first tick mark

    DECLARE DEF h2Within                         ! window/function transforms
    DECLARE DEF h2AnimWithin, h2LstpWithin, h2RstpWithin, h2AnimStopWithin
    DECLARE DEF h2LmoveWithin, h2RmoveWithin

    CALL h2SliderVariables

    SUB h2Init
        CALL h2SliderVariables
        CALL h2DrawSlider(h2name$,xa)
        CALL h2AnimButtons
        LET h2AnimStep= h2Mult*(h2fRgt-h2fLft) / (h2wRgt-h2wLft)
    END SUB

    SUB h2SetBounds(eqnum)
        SELECT CASE eqnum
        CASE 3,6,8,9,16
             LET h2PiAxis = 1
             LET h2Mult   = pi
        CASE else
             LET h2PiAxis = 0
             LET h2Mult   = 1
        END SELECT

        ! LET h2AnimStep= h2Mult/100

        LET h2fLft    = w1fLft
        LET h2fRgt    = w1fRgt
        LET h2AnimStep= 2 * (h2Mult*(h2fRgt-h2fLft) / (h2wRgt-h2wLft))
        LET h2STik    = 0.5                      ! 0.5                             ! short tick marks
        LET h2LTik    = 1                        ! long tick marks
        LET h2Label   = 1                        ! labels
        LET h2First   = w1xFirst                 ! first tick mark
        LET h2Click   = 0.5*h2Mult
        CALL h2SliderVariables
    END SUB

    ! ---

    PUBLIC errors(0:639)
    MAT redim errors(w1Lft:w1Rgt)

    ! ----------- coefficients --------
    ! ---------- Poly and Polya -------

    SHARE num$,form$,lnspc,topx,topxa,txtlft,txtrgt,bas,sign,lng,clrlft

    ! LET form$ = "------%.########"

    LET form$ = "+++++++%.######"
    LET txtrgt= workrgt- 60
    LET txtlft= workrgt- 55
    LET lnspc = 18
    LET topx  = w1Top + 10 + lnspc
    LET topxa = topx
    LET clrlft= w1Rgt + 20
    LET clrrgt= workRgt
    LET clrbas= topx + 10*lnspc
    LET clrtop= topx - 15
    LET atxclr= yellow

    LET polybtnlft= workrgt-160                  ! on/off btn
    LET polybtnrgt= polybtnlft+90
    LET polybtnbas= w1Bas
    LET polybtntop= w1Bas-17
    LET polystate = 0                            ! reg on/off flag

    SUB ClearTerms
        BOX CLEAR clrLft,clrRgt,clrBas,clrTop
    END SUB

    SUB PolyTextButtons
        CALL DrawButton(polybtnlft,polybtnrgt,polybtnbas,polybtntop,5,"Show Terms")
        !CALL DrawButton(polyabtnlft,polyabtnrgt,polyabtnbas,polyabtntop,0,"")
    END SUB

    SUB PolySwitch
        LET polystate= 1-polystate
        IF polystate=0 then
           CALL ClearPoly
           BOX CLEAR clrLft,clrRgt,topx,topx-30
        ELSE
           CALL PolyInit
        END IF
    END SUB

    SUB PolyInit
        CALL PlotTextLJ(clrlft+20,w1Top+8,"p(x) = ",yellow)
        CALL ShowEquations(xa)
    END SUB

    SUB ShowEquations(xa)
        IF polystate=1 then
!           IF xa=0 THEN
!              CALL ShowTerms(eqnum,n)
!              CALL ShowCoef(eqnum,xa,ya,n)
!           ELSE
              CALL ShowTermsa(eqnum,n)
              CALL ShowCoefa(eqnum,xa,ya,n)
!           END IF
        ELSE
           CALL ClearPoly
        END IF
    END SUB

    SUB ClearPoly
        ! BOX CLEAR clrlft,workrgt,topx+10*lnspc,topx-30
        BOX CLEAR clrLft,clrRgt,clrBas,clrTop
    END SUB

    SUB ShowTerms(eqnum,n)
        CALL ClearPoly
        CALL SetTextFont(1,12,"bold")

        FOR i= 0 to n
            SELECT CASE i
            CASE 0
                 LET txt$= ""
            CASE 1
                 LET txt$= "x"
            CASE else
                 IF i=n then
                    LET txt$= "x^[" & str$(i) & "]"
                 ELSE
                    LET txt$= "x^[" & str$(i) & "]"
                 END IF
            END SELECT
            LET  bas= topx + i*lnspc
            CALL SuperSubScriptLJ(txtlft,bas,txt$,atxclr)
        NEXT i
    END SUB

    SUB ShowCoef(eqnum,xa,ya,n)
        LOCAL c[0:10],a[10]

        BOX CLEAR clrlft,txtlft-3,clrBas,clrTop
        CALL SetTextFont(1,12,"bold")

        LET c[0]= ya
        LET num$= using$(form$,c[0])
        LET txt$= ""
        LET bas = topx
        ! CALL Text(2,0,txtrgt,bas,num$,atxclr)
        CALL PlotTextRJ(txtrgt,bas,num$,atxclr)

        CALL GetNine(eqnum,xa,ya,c[1],c[2],c[3],c[4],c[5],c[6],c[7],c[8],c[9])

        LET a[1] = xa
        FOR I = 2 TO n
            LET a[I] = a[I-1]*xa                 !a[n] = a^n
        NEXT I
        FOR I = 2 TO n
            LET c[I] = c[I]/fact[I]
        NEXT I

        LET c[0]= ya
        LET sign= -1
        FOR I= 1 to n
            LET c[0]= c[0] + sign*a[I]*c[I]
            LET sign= -sign
        NEXT I

        IF n>=1 THEN
           LET sign= -1
           FOR I= 2 to n
               LET c[1]= c[1] + sign*I*a[I-1]*c[I]
               LET sign= -sign
           NEXT I
           LET num$= using$(form$,c[1])
           LET bas = topx + lnspc
           ! CALL Text(2,0,txtrgt,bas,num$,atxclr)
           CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        END IF

        IF n>=2 THEN
           LET sign= -1
           FOR I= 3 to n
               LET c[2]= c[2] + sign*I*(I-1)*a[I-2]*c[I]/2
               LET sign= -sign
           NEXT I
           LET num$= using$(form$,c[2])
           LET bas = topx + 2*lnspc
           ! CALL Text(2,0,txtrgt,bas,num$,atxclr)
           CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        END IF

        IF n>=3 THEN
           LET sign= -1
           FOR I= 4 to n
               LET c[3]= c[3] + sign*I*(I-1)*(I-2)*a[I-3]*c[I]/6
               LET sign= -sign
           NEXT I
           LET num$= using$(form$,c[3])
           LET bas = topx + 3*lnspc
           CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        END IF

        IF n>=4 THEN
           LET sign= -1
           FOR I= 5 to n
               LET c[4]= c[4] + sign*I*(I-1)*(I-2)*(I-3)*a[I-4]*c[I]/24
               LET sign= -sign
           NEXT I
           LET num$= using$(form$,c[4])
           LET bas = topx + 4*lnspc
           CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        END IF

        IF n=5 THEN
           LET num$= using$(form$,c[5])
           LET bas = topx + 5*lnspc
           CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        END IF

        IF n=6 THEN
           LET c[5]= -6*a[1]*c[6]+c[5]

           FOR i= 5 to 6
               LET num$= using$(form$,c[i])
               LET bas = topx + i*lnspc
               CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
           NEXT i
        END IF

        IF n=7 THEN
           LET c[5]= +21*a[2]*c[7]-6*a[1]*c[6]+c[5]
           LET c[6]= -7*a[1]*c[7]+c[6]

           FOR i= 5 to 7
               LET num$= using$(form$,c[i])
               LET bas = topx + i*lnspc
               CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
           NEXT i
        END IF

        IF n=8 THEN
           LET c[5]= -56*a[3]*c[8]+21*a[2]*c[7]-6*a[1]*c[6]+c[5]
           LET c[6]= +28*a[2]*c[8]-7*a[1]*c[7]+c[6]
           LET c[7]= -8*a[1]*c[8]+c[7]

           FOR i= 5 to 8
               LET num$= using$(form$,c[i])
               LET bas = topx + i*lnspc
               CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
           NEXT i
        END IF

        IF n=9 THEN
           LET c[5]= 126*a[4]*c[9]-56*a[3]*c[8]+21*a[2]*c[7]-6*a[1]*c[6]+c[5]
           LET c[6]= -84*a[3]*c[9]+28*a[2]*c[8]-7*a[1]*c[7]+c[6]
           LET c[7]= 36*a[2]*c[9]-8*a[1]*c[8]+c[7]
           LET c[8]= -9*a[1]*c[9]+c[8]

           FOR i= 5 to 9
               LET num$= using$(form$,c[i])
               LET bas = topx + i*lnspc
               CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
           NEXT i
        END IF
    END SUB

    SUB ShowTermsa(eqnum,n)
        CALL ClearPoly
        CALL SetTextFont(1,12,"bold")

        FOR i= 0 to n
            SELECT CASE i
            CASE 0
                 LET txt$= ""
            CASE 1
                 LET txt$= "(x-a)"
            CASE else
                 IF i=n then
                    LET txt$= "(x-a)^[" & str$(i) & "]"
                 ELSE
                    LET txt$= "(x-a)^[" & str$(i) & "]"
                 END IF
            END SELECT
            LET  bas= topxa + i*lnspc
            CALL SuperSubScriptLJ(txtlft,bas,txt$,atxclr)
        NEXT i
    END SUB

    SUB ShowCoefa(eqnum,xa,ya,n)                 ! new
        LOCAL c[0:10]
        CALL  GetNine(eqnum,xa,ya,c[1],c[2],c[3],c[4],c[5],c[6],c[7],c[8],c[9])

        LET c[0]= ya
        FOR I= 2 TO n
            LET c[I] = c[I]/fact[I]
        NEXT I

        BOX CLEAR clrlft,txtlft-3,clrBas,clrTop
        CALL SetTextFont(1,12,"bold")

        FOR i= 0 to n
            LET num$= using$(form$,c[i])
            LET bas = topxa + i*lnspc
            ! CALL Text(2,0,txtrgt,bas,num$,atxclr)
            CALL PlotTextRJ(txtrgt,bas,num$,atxclr)
        NEXT i
    END SUB


    ! ------------- Text output holes -------------

    LET xtxtlft= w1Midx - 18
    LET xtxtrgt= xtxtlft+100
    LET xtxtbas= h2wBas + 20
    LET xtxttop= xtxtbas- 10

    ! ---- formulas text ----

    ! p(x) = a_0 + a_1(x-a) + a_2(x-a)^2 + ...
    ! the coefficient a_k  is f^{(k)}(a)/k!. So either the polynomial
    ! in this form, with numbers, or perhaps a list of values of the a_k's.
    ! LET tp$= "p(x) = f^[0](a)/0! + f^[1](a)/1!(x-a) + f^[2](a)/2!(x-a) + f^[3](a)/3!(x-a) + ... )" 
    ! LET tp$= "f(a)/0" & excl$ & " + (f'(a)/1" & excl$ & ") (x-a) + (f''(a)/2" & excl$ & ") (x-a)^[2] + (f'''(a)/3" & excl$ & ") (x-a)^[3] + (f''''(a)/4" & excl$ & ") (x-a)^[4] + ... )"

    SUB ShowFormulas(eqFlag)
        IF eqFlag=1 then
           LET excl$= chr$(33)
           LET c$ = "k_[0] = f(a)" & "      k_[1] = f'(a)/1" & excl$ & "      k_[2] = f''(a)/2" & excl$  & "      k_[3] = f'''(a)/3" & excl$  & "   ..."  ! ?? notation for f', f'', etc ??
           CALL SuperSubScriptLJ(w1Lft+50,w1Top-65,c$,yellow)
           LET c$ = "k_"
           LET tp$= c$ & "[0] + " & c$ & "[1]" & "(x-a) + " & c$ & "[2]" & "(x-a)^[2] + " & c$ & "[3]" & "(x-a)^[3] + " & c$ & "[4]" & "(x-a)^[4] + ..."
           CALL PlotTextRJ(w1Lft+50,w1Top-40,"p(x) = ",yellow)
           CALL SuperSubScriptLJ(w1Lft+50,w1Top-40,tp$ ,yellow)
        ELSE
           BOX CLEAR w1Lft,workrgt,w1Top-25,worktop
        END IF
    END SUB


    ! --- text rectangle 4 - rollover t,x values ---

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

    LET t4Label1$= "x = "
    LET t4Label2$= "y = "
    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 x = w1Fncx(mx)
        LET y = w1Fncy(my)
        LET x$= trim$(using$("--%.##",x))
        LET y$= trim$(using$("--%.##",y))
        CALL t4Clear
        CALL PlotTextLJ(t4Eqx,t4BasLn1,x$,t4Clr)
        CALL PlotTextLJ(t4Eqx,t4BasLn2,y$,t4Clr)
    END SUB

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

    SUB t4Init
        CALL t4Label
        ! CALL t4SetCoords
    END SUB


    ! ----- menu 1 data -----

    DECLARE PUBLIC m1Lft, m1Rgt, m1Bas, m1Top, m1Equation, m1Prefix$, m1Menu1$()
    DECLARE PUBLIC m1tClr, m2tClr

    DECLARE DEF m1Within

    MAT redim m1Menu1$(1:17)
    MAT READ m1Menu1$
    DATA "1/sqrt(1+x)"
    DATA "0.3x^[2] - .01x^[5]"
    DATA "sin(x)"
    DATA "e^[x]"
    DATA "1/(2x)"
    DATA "cos(x)"
    DATA "sin(x^[2])"
    DATA "tan(x)"
    DATA "sec(x)"
    DATA "-2x/(1+x^[2])^[2]"
    DATA "1/(1+x^[2])"
    DATA "tan^[-1](x)"
    DATA "x tan^[-1](x) - ln(1+x^[2])/2"
    DATA "sin^[-1](x)"
    DATA "log(x)"
    DATA "cos(x) + cos(9x)/9"
    DATA "exp(-1/x^[2])    [x>=0]"

    LET m1Prefix$ = "f(x) = "
    LET m1Lft     = w1Lft
    LET m1Top     = workbas - 30
    LET m1tClr    = white


    ! --- buttons ---

    LET btop= m1Top
    LET bbas= btop + 18
    LET blft= h2wRgt + 70
    LET brgt= blft + 80

    LET polybtnlft= brgt+2
    LET polybtnrgt= polybtnlft + 55
    LET polybtntop= btop
    LET polybtnbas= bbas

    SUB Buttons
        CALL SetTextFont(1,12,"bold")
        CALL DrawButton(blft,brgt,bbas,btop,5,"Formulas")
        CALL DrawButton(polybtnlft,polybtnrgt,polybtnbas,polybtntop,5,"Terms")
    END SUB



    ! --- default parameters ---

    LET n,oldn    = 1
    LET eq        = 3
    CALL FixEqnum(eq)                            ! outputs eqnum

    LET m1Equation= eq
    LET polystate = 0
    ! LET polyastate= 0
    LET xa,oldxa  = 3
    LET eqFlag    = 0
    LET m1State   = 0

    ! --- draw the screen ---

    CALL InitScreen
    CALL m1InitMenu1

    SUB InitScreen
        BOX CLEAR worklft,workrgt,workbas,worktop

        CALL ShowFormulas(eqFlag)

        CALL w1SetBounds(eqnum)
        CALL h2SetBounds(eqnum)
        CALL w1Init
        CALL h2Init
        CALL h1Init
        CALL m1ResetMenu(m1State)
        CALL DrawFunction(1)
        LET  wxa= w1Wndx(xa)
        CALL RollOver1(xa,wxa,0)
        CALL Buttons

        CALL t4Init

        CALL ShowEquations(xa)

        IF polystate=1 then CALL PolyInit
        !    CALL ShowTerms(eqnum,n)
        !    CALL ShowCoef(eqnum,xa,ya,n)
    END SUB

    ! ------------- Main event loop -------------

    DO
       LET clearFlag= 0
       IF ms<>2 then
          BOX KEEP w1Lft-5,w1Rgt+5,w1Bas+5,w1Top-5 in temp$
          DO                                     ! while mouse is up
             GET MOUSE: mx,my,ms

             IF w1wWithin(mx,my)=true then
                IF (mx<>oldmx or my<>oldmy) then
                   ! CALL w1ShowGraphLayer
                   BOX SHOW temp$ at w1Lft-5,w1Bas+5
                   SET COLOR litmid
                   PLOT w1Lft+1,my; w1Rgt-1,my
                   PLOT mx,w1Top+1; mx,w1Bas-1
                   LET  t = w1Fncx(mx)
                   LET  x = w1Fncy(my)
                   CALL t4SetCoords

                   LET w1ClearFlag= 1
                   LET oldmx= mx
                   LET oldmy= my
                END IF
             ELSE IF w1ClearFlag=1 then
                BOX SHOW temp$ at w1Lft-5,w1Bas+5
                CALL t4Clear
                LET w1ClearFlag= 0
             END IF
          LOOP until ms=2                        ! mouse goes down
       END IF

       IF w1ClearFlag=1 then
          BOX SHOW temp$ at w1Lft-5,w1Bas+5
          CALL t4Clear
          LET w1ClearFlag= 0
       END IF

       IF h1Within(mx,my)=true then              ! n (taylor terms)
          IF my<h1axis then
             CALL h1ClickEvent(ms)
          ELSE
             CALL h1DragEvent(ms)
          END IF

       ELSE IF h2Within(mx,my)=true then         ! x slider
          IF my<h2axis then
             CALL h2ClickEvent(ms)
          ELSE
             CALL h2DragEvent(ms)
          END IF

       ELSE IF h2AnimWithin(mx,my)=true then
          CALL h2AnimButtonUp(ms)
          CALL AnimateK
       ELSE IF h2LstpWithin(mx,my)=true then
          CALL h2LstpButtonUp(ms)
          CALL h2LeftStep(xa,h2AnimStep)
          CALL h2Action
       ELSE IF h2RstpWithin(mx,my)=true then
          CALL h2RstpButtonUp(ms)
          CALL h2RightStep(xa,h2AnimStep)
          CALL h2Action

       ELSE IF m1Within(mx,my)=true then         ! menu
          CALL m1ButtonUp(ms)
          CALL m1ReadMenu
          LET eq= m1Equation
          CALL FixEqnum(eq)

          CALL w1SetBounds(eqnum)
          CALL h2SetBounds(eqnum)
          CALL w1Init
          CALL h2Init
          LET h2AnimStep= (h2fRgt-h2fLft)/(h2wRgt-h2wLft)
          CALL DrawFunction(1)

          CALL RollOver1(xa,wxa,0)
          IF polystate=1 then
             CALL ShowTermsa(eqnum,n)
             CALL ShowCoef(eqnum,xa,ya,n)
          END IF

          IF polyastate=1 THEN
             CALL ShowCoefa(eqnum,xa,ya,n)
          END IF

       ELSE IF mx>polybtnlft and mx<polybtnrgt and my>polybtntop and my<polybtnbas then
          CALL MouseButtonUp(polybtnlft,polybtnrgt,polybtnbas,polybtntop,ms)
          CALL polyswitch

       ELSE IF mx>blft and mx<brgt and my>btop and my<bbas then      ! formula button
          CALL MouseButtonUp(blft,brgt,bbas,btop,ms)
          LET eqFlag= 1-eqFlag
          CALL ShowFormulas(eqFlag)

       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

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

    SUB ClearAll
        CALL ClearRollOver1
        IF polystate =1 then CALL ClearPoly
    END SUB

    ! --- mouse events ---

    SUB h1ClickEvent(ms)
        CALL h1GetClickVal(ms,h1Click,n)
        CALL h1Action
        CALL h1MouseUpAction
    END SUB
    SUB h1DragEvent(ms)
        DO
           CALL h1GetDragVal(ms,h1Places,n)
           CALL h1Action
        LOOP until ms=3
        CALL h1MouseUpAction
    END SUB

    SUB h1Action
        IF n<>oldn then
           CALL RollOver1(xa,wxa,0)
           LET oldn= n
        END IF
    END SUB
    SUB h1MouseUpAction
        CALL RollOver1(xa,wxa,0)
        IF polystate=1 THEN
           CALL ShowCoef(eqnum,xa,ya,n)
           CALL ShowEquations(xa)
        END IF
    END SUB

    ! ---

    SUB h2ClickEvent(ms)
        CALL h2GetClickVal(ms,h2Click,xa)
        CALL h2Action
    END SUB
    SUB h2DragEvent(ms)
        DO
           CALL h2GetDragVal(ms,h2Places,xa)
           CALL h2Action
        LOOP until ms=3
    END SUB

    SUB h2Action
        IF xa<>oldxa then
           LET wxa= w1Wndx(xa)
           CALL RollOver1(xa,wxa,0)
           IF polystate=1 THEN
              CALL ShowCoef(eqnum,xa,ya,n)
              !IF (xa=0 and oldxa<>0) or (xa<>0 and oldxa=0) then
                 CALL ShowEquations(xa)
              !END IF
           END IF
           LET oldxa= xa
        END IF
    END SUB

    ! --- drawing methods ---

    SUB RollOver1(xa,wxa,stpflg)
        LOCAL wy,d,cx
        CALL ClearRollOver1
        CALL GetFnc(eqnum,xa,ya)

        ! CALL Convergence

        CALL DrawTaylor(n,xa,ya,stpflg)

        IF eqnum=17 and xa<0 then
           LET wy= w1Wndy(0)
        ELSE
           LET wy= w1Wndy(ya)
        END IF

        IF wy>=w1Top and wy<=w1Bas and ya<>badarg  then
           DRAW diamond5(cyan) with shift(wxa,wy)
        END IF
    END SUB
    SUB ClearRollOver1
        CALL w1ShowGraphLayer
    END SUB

    SUB Convergence
        ! interval of convergence
        SET COLOR 24
        SELECT CASE eqnum
        CASE 1
             IF xa>-1 THEN
                LET cx= w1Wndx(-1)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(2*xa+1)
                IF cx>w1Lft and cx<w1Rgt then
                   PLOT cx,w1Top+1; cx,w1Bas-1
                END IF
             END IF
        CASE 5
             LET cx= w1Wndx(2*xa)
             IF cx>w1Lft and cx<w1Rgt then
                PLOT cx,w1Top+1; cx,w1Bas-1
             END IF
             LET cx= w1Wndx(0)
             PLOT cx,w1Top+1; cx,w1Bas-1
        CASE 8,9
             IF -pi/2<xa and xa <= 0 THEN
                LET cx= w1Wndx(-pi/2) + 1
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(2*xa + pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
             ELSE IF 0<xa and xa < pi/2 - .01 THEN
                LET cx= w1Wndx(2*xa - pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
             ELSE IF pi/2 + .01 < xa and xa <= pi THEN
                LET cx= w1Wndx(pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(2*xa - pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
             ELSE IF  pi < xa and xa < 3*pi/2 THEN
                LET cx= w1Wndx(2*xa - 3*pi/2)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(3*pi/2) - 1
                PLOT cx,w1Top+1; cx,w1Bas-1
             END IF
        CASE 10,11,12,13
             LET cx= w1Wndx(xa-sqr(1+xa*xa))
             IF cx>w1Lft and cx<w1Rgt then
                PLOT cx,w1Top+1; cx,w1Bas-1
             END IF
             LET cx= w1Wndx(xa+sqr(1+xa*xa))
             IF cx>w1Lft and cx<w1Rgt then
                PLOT cx,w1Top+1; cx,w1Bas-1
             END IF
        CASE 14
             IF -1 < xa and xa<0 THEN
                LET cx= w1Wndx(-1)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(2*xa+1)
                PLOT cx,w1Top+1; cx,w1Bas-1
             ELSE IF xa=0 THEN
                LET cx= w1Wndx(-1)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(1)
                PLOT cx,w1Top+1; cx,w1Bas-1
             ELSE IF 0 < xa and xa < 1 THEN
                LET cx= w1Wndx(2*xa-1)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(1)
                PLOT cx,w1Top+1; cx,w1Bas-1
             END IF
        CASE 15
             IF xa>0 THEN
                LET cx= w1Wndx(0)
                PLOT cx,w1Top+1; cx,w1Bas-1
                LET cx= w1Wndx(2*xa)
                IF cx>w1Lft and cx<w1Rgt then
                   PLOT cx,w1Top+1; cx,w1Bas-1
                END IF
             END IF
        CASE 17
             IF xa<>0 THEN
                LET cx= w1Wndx(2*xa)
                IF cx>w1Lft and cx<w1Rgt then
                   PLOT cx,w1Top+1; cx,w1Bas-1
                END IF
                PLOT w1Wndx(0),w1Top+1; w1Wndx(0),w1Bas-1
             END IF
        CASE else
        END SELECT
    END SUB

    ! --- Animate ---

    SUB AnimateK
        DO
           CALL h2RightStep(xa,h2AnimStep)
           IF xa>h2fMax-h2AnimStep then
              EXIT DO
           ELSE
              CALL h2Action
           END IF

           GET MOUSE: mx,my,ms
           IF ms=3 then
              EXIT DO
           END IF
           CALL Delay(1/16)
        LOOP
        CALL h2StopButtonClear
    END SUB

    ! --- graphing routines ---

    SUB DrawTaylor(n,xa,ya,stpflg)
        CALL GetNine(eqnum,xa,ya,d1,d2,d3,d4,d5,d6,d7,d8,d9)
        SET COLOR yellow
        LET stp = 1
        CALL DrawLoopBnd(+stp,wx)
        LET maxx= wx
        CALL DrawLoopBnd(-stp,wx)
        LET minx= wx
    END SUB

    SUB DrawLoopBnd(step,wx)
        ! clipping here
        LOCAL x,yt,wy
        LET wx= wxa
        LET wy= w1Wndy(ya)

        IF eqnum<>17 or wx>=w1x0 then
           DO
              LET wx= wx+step
              IF wx>w1Lft and wx<w1Rgt then
                 LET x = w1Fncx(wx)
                 LET yt= Taylor(n,x,xa,ya)
                 LET wy= w1Wndy(yt)
                 IF wy>=w1Top and wy<=w1Bas then
                    PLOT wx,wy;
                 ELSE
                    PLOT
                 END IF
              ELSE
                 PLOT
                 EXIT DO
              END IF
           LOOP
           PLOT
        END IF
        IF eqnum=17 and wxa<w1x0 then
           PLOT w1Lft,w1y0; w1Rgt,w1y0
        END IF
    END SUB

    !    SUB ClipPlot
    !        LET oldwy= wy
    !        IF y>w1fBas and y<w1fTop and oldy>w1fBas and oldy<w1fTop then
    !           LET wy= w1Wndy(y)
    !           PLOT wx,wy; wx,oldwy
    !        ELSE IF oldy>w1fBas and oldy<w1fTop then
    !           IF y>=w1fTop then
    !              PLOT wx,w1Top; wx,oldwy
    !           ELSE IF y<=w1fBas then
    !              PLOT wx,w1Bas; wx,oldwy
    !           END IF
    !        ELSE IF y>w1fBas and y<w1fTop then
    !           LET wy= w1Wndy(y)
    !           IF oldy>=w1fTop then
    !              PLOT wx,wy; wx,w1Top
    !           ELSE IF oldy<=w1fBas then
    !              PLOT wx,wy; wx,w1Bas
    !           END IF
    !        END IF
    !    END SUB

    SUB DrawLoopErr(stp,maxerr,wx)
        LOCAL x,y,yt,wy,err
        LET wx= wxa
        LET wy= w1Wndy(ya)
        IF wy>=w1Top and wy<w1Bas then
           PLOT wx,wy;
        END IF
        DO
           LET wx= wx+stp
           IF wx>w1Lft and wx<w1Rgt then
              LET x  = w1Fncx(wx)
              CALL GetFnc(eqnum,x,y)
              LET yt = Taylor(n,x,xa,ya)
              LET err= yt-y
              LET errors(wx)= err
              LET wy = w1Wndy(yt)
              SET COLOR yellow
              IF wy>=w1Top and wy<=w1Bas and y<>badarg then
                 PLOT wx,wy;
              ELSE
                 PLOT
              END IF
           ELSE
              PLOT
              EXIT DO
           END IF
        LOOP until abs(err)>maxerr
        PLOT

        IF stp>0 then                            ! new
           LET wx = wx-1
        ELSE
           LET wx = wx+1
        END IF
        PLOT wx,w1Bas; wx,w1Bas-5
        PLOT
    END SUB

    SUB DrawFunction(stp)
        LOCAL x,y,wx,wy
        CALL w1ShowGridLayer
        SET COLOR cyan

        IF eqnum=17 then                         ! zero for negative input
           PLOT w1Lft,w1y0; w1x0,w1y0
           FOR wx= w1x0 to w1Rgt                 ! step stp
               LET x= w1Fncx(wx)
               CALL GetFnc(eqnum,x,y)
               IF y<>badarg and y<>bignum then
                  LET wy= w1Wndy(y)
                  IF wy>=w1Top and wy<=w1Bas then
                     PLOT wx,wy;
                  ELSE
                     PLOT
                  END IF
               ELSE
                  PLOT
               END IF
           NEXT wx
        ELSE
           LET x= w1fLft
           CALL GetFnc(eqnum,x,y)
           FOR wx= w1Lft to w1Rgt                ! step stp
               LET oldx= x
               LET oldy= y
               LET x   = w1Fncx(wx)
               CALL GetFnc(eqnum,x,y)

               IF y<>badarg and y<>bignum then
                  LET oldwy= wy
                  LET wy   = w1Wndy(y)
                  IF wy>w1Top and wy<w1Bas and oldwy>w1Top and oldwy<w1Bas then
                     PLOT wx,wy; wx-1,oldwy
                  ELSE
                     IF (oldy>w1fTop and y>w1fTop) or (oldy<w1fBas and y<w1fBas) then
                     ELSE
                        IF (oldy>w1fTop and y<w1fBas) or (oldy<w1fBas and y>w1fTop) then

                        ELSE
                           LET px= oldx
                           LET py= oldy
                           LET qx= x
                           LET qy= y
                           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
                     END IF
                  END IF
               END IF
           NEXT wx
        END IF
        PLOT
        CALL w1KeepGraphLayer
    END SUB
END SUB

! ----- end of nth Taylor polynomials -----