axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PLOT.spad

)abbrev domain PLOT Plot
++ Author: Michael Monagan (revised by Clifton J. Williamson)
++ Date Created: Jan 1988
++ Date Last Updated: 30 Nov 1990 by Jonathan Steinbach
++ Description:
++ The Plot domain supports plotting of functions defined over a
++ real number system.  A real number system is a model for the real
++ numbers and as such may be an approximation.  For example
++ floating point numbers and infinite continued fractions.
++ The facilities at this point are limited to 2-dimensional plots
++ or either a single function or a parametric function.

Plot() : SIG == CODE where

  B   ==> Boolean
  F   ==> DoubleFloat
  I   ==> Integer
  L   ==> List
  N   ==> NonNegativeInteger
  OUT ==> OutputForm
  P   ==> Point F
  RN  ==> Fraction Integer
  S   ==> String
  SEG ==> Segment
  R   ==> Segment F
  C   ==> Record(source: F -> P,ranges: L R,knots: L F,points: L P)

  SIG ==> PlottablePlaneCurveCategory with

    --% function plots

    plot : (F -> F,R) -> %
      ++ plot(f,a..b) plots the function \spad{f(x)} 
      ++ on the interval \spad{[a,b]}.
      ++ 
      ++X fp:=(t:DFLOAT):DFLOAT +-> sin(t)
      ++X plot(fp,-1.0..1.0)$PLOT

    plot : (F -> F,R,R) -> %
      ++ plot(f,a..b,c..d) plots the function \spad{f(x)} on the interval
      ++ \spad{[a,b]}; y-range of \spad{[c,d]} is noted in Plot object.

    --% multiple function plots

    plot : (L(F -> F),R) -> %
      ++ plot([f1,...,fm],a..b) plots the functions \spad{y = f1(x)},...,
      ++ \spad{y = fm(x)} on the interval \spad{a..b}.

    plot : (L(F -> F),R,R) -> %
      ++ plot([f1,...,fm],a..b,c..d) plots the functions \spad{y = f1(x)},...,
      ++ \spad{y = fm(x)} on the interval \spad{a..b}; y-range of \spad{[c,d]}
      ++ is noted in Plot object.

    --% parametric plots

    plot : (F -> F,F -> F,R) -> %
      ++ plot(f,g,a..b) plots the parametric curve \spad{x = f(t)},
      ++ \spad{y = g(t)} as t ranges over the interval \spad{[a,b]}.

    plot : (F -> F,F -> F,R,R,R) -> %
      ++ plot(f,g,a..b,c..d,e..f) plots the parametric curve \spad{x = f(t)},
      ++ \spad{y = g(t)} as t ranges over the interval \spad{[a,b]}; x-range
      ++ of \spad{[c,d]} and y-range of \spad{[e,f]} are noted in Plot object.

    --% parametric plots

    pointPlot : (F -> P,R) -> %
      ++ pointPlot(t +-> (f(t),g(t)),a..b) plots the parametric curve
      ++ \spad{x = f(t)}, \spad{y = g(t)} as t ranges over the interval
      ++ \spad{[a,b]}.

    pointPlot : (F -> P,R,R,R) -> %
      ++ pointPlot(t +-> (f(t),g(t)),a..b,c..d,e..f) plots the parametric
      ++ curve \spad{x = f(t)}, \spad{y = g(t)} as t ranges over the interval
      ++ \spad{[a,b]}; x-range of \spad{[c,d]} and y-range of \spad{[e,f]}
      ++ are noted in Plot object.

    --% polar plots

    plotPolar : (F -> F,R) -> %
      ++ plotPolar(f,a..b) plots the polar curve \spad{r = f(theta)} as
      ++ theta ranges over the interval \spad{[a,b]}; this is the same as
      ++ the parametric curve \spad{x = f(t)*cos(t)}, \spad{y = f(t)*sin(t)}.

    plotPolar : (F -> F) -> %
      ++ plotPolar(f) plots the polar curve \spad{r = f(theta)} as theta
      ++ ranges over the interval \spad{[0,2*%pi]}; this is the same as
      ++ the parametric curve \spad{x = f(t)*cos(t)}, \spad{y = f(t)*sin(t)}.

    plot : (%,R) -> %              -- change the range
        ++ plot(x,r) is not documented

    parametric? : % -> B
      ++ parametric? determines whether it is a parametric plot?

    zoom : (%,R) -> %
        ++ zoom(x,r) is not documented

    zoom : (%,R,R) -> %
        ++ zoom(x,r,s) is not documented

    refine : (%,R) -> %
        ++ refine(x,r) is not documented

    refine : % -> %
      ++ refine(p) performs a refinement on the plot p

    tRange : % -> R
      ++ tRange(p) returns the range of the parameter in a parametric plot p

    minPoints : () -> I
      ++ minPoints() returns the minimum number of points in a plot

    setMinPoints : I -> I
      ++ setMinPoints(i) sets the minimum number of points in a plot to i

    maxPoints : () -> I
      ++ maxPoints() returns the maximum number of points in a plot

    setMaxPoints : I -> I
      ++ setMaxPoints(i) sets the maximum number of points in a plot to i

    screenResolution : () -> I
      ++ screenResolution() returns the screen resolution

    setScreenResolution : I -> I
      ++ setScreenResolution(i) sets the screen resolution to i

    adaptive? : () -> B
      ++ adaptive?() determines whether plotting be done adaptively

    setAdaptive : B -> B
      ++ setAdaptive(true) turns adaptive plotting on
      ++ \spad{setAdaptive(false)} turns adaptive plotting off

    numFunEvals : () -> I
      ++ numFunEvals() returns the number of points computed

    debug : B -> B
      ++ debug(true) turns debug mode on
      ++ \spad{debug(false)} turns debug mode off

  CODE ==> add

    import PointPackage(DoubleFloat)

--% local functions

    checkRange     : R -> R
      -- checks that left-hand endpoint is less than right-hand endpoint
    intersect      : (R,R) -> R
      -- intersection of two intervals
    union          : (R,R) -> R
      -- union of two intervals
    join           : (L C,I) -> R
    parametricRange: % -> R
    select         : (L P,P -> F,(F,F) -> F) -> F
    rangeRefine    : (C,R) -> C
    adaptivePlot   : (C,R,R,R,I) -> C
    basicPlot      : (F -> P,R) -> C
    basicRefine    : (C,R) -> C
    pt             : (F,F) -> P
    Fnan?           : F -> Boolean
    Pnan?           : P -> Boolean

--% representation

    Rep := Record( parametric: B, _
                   display: L R, _
                   bounds: L R, _
                   axisLabels: L S, _
                   functions: L C )

--% global constants

    ADAPTIVE: B := true
    MINPOINTS: I := 49
    MAXPOINTS: I := 1000
    NUMFUNEVALS: I := 0
    SCREENRES: I := 500
    ANGLEBOUND: F := cos inv (4::F)
    DEBUG: B := false

    Fnan?(x) == x ~= x

    Pnan?(x) == any?(Fnan?,x)

--% graphics output

    listBranches plot ==
      outList : L L P := nil()
      for curve in plot.functions repeat
        -- curve is C
        newl:L P:=nil()
        for p in curve.points repeat
          if not Pnan? p then newl:=cons(p,newl)
          else if not empty? newl then 
                outList := concat(newl:=reverse! newl,outList)
                newl:=nil()
        if not empty? newl then outList := concat(newl:=reverse! newl,outList)
      outList

    checkRange r == (lo r > hi r => error "ranges cannot be negative"; r)

    intersect(s,t) == checkRange (max(lo s,lo t) .. min(hi s,hi t))

    union(s,t) == min(lo s,lo t) .. max(hi s,hi t)

    join(l,i) ==
      rr := first l
      u : R :=
        i = 0 => first(rr.ranges)
        i = 1 => second(rr.ranges)
        third(rr.ranges)
      for r in rest l repeat
        i = 0 => u := union(u,first(r.ranges))
        i = 1 => u := union(u,second(r.ranges))
        u := union(u,third(r.ranges))
      u
    parametricRange r == first(r.bounds)

    minPoints() == MINPOINTS

    setMinPoints n ==
      if n < 3 then error "three points minimum required"
      if MAXPOINTS < n then MAXPOINTS := n
      MINPOINTS := n

    maxPoints() == MAXPOINTS

    setMaxPoints n ==
      if n < 3 then error "three points minimum required"
      if MINPOINTS > n then MINPOINTS := n
      MAXPOINTS := n

    screenResolution() == SCREENRES

    setScreenResolution n ==
      if n < 2 then error "buy a new terminal"
      SCREENRES := n

    adaptive?() == ADAPTIVE

    setAdaptive b == ADAPTIVE := b

    parametric? p == p.parametric

    numFunEvals() == NUMFUNEVALS

    debug b == DEBUG := b

    xRange plot == second plot.bounds

    yRange plot == third plot.bounds

    tRange plot == first plot.bounds

    select(l,f,g) ==
      m := f first l
      if Fnan? m then m := 0
      for p in rest l repeat
        n := m
        m := g(m, f p)
        if Fnan? m then m := n
      m

    rangeRefine(curve,nRange) ==
      checkRange nRange; l := lo nRange; h := hi nRange
      t := curve.knots; p := curve.points; f := curve.source
      while not null t and first t < l repeat
        (t := rest t; p := rest p)
      c: L F := nil(); q: L P := nil()
      while not null t and (first t) <= h repeat
        c := concat(first t,c); q := concat(first p,q)
        t := rest t; p := rest p
      if null c then return basicPlot(f,nRange)
      if first c < h then
        c := concat(h,c)
        q := concat(f h,q)
        NUMFUNEVALS := NUMFUNEVALS + 1
      t := c := reverse_! c; p := q := reverse_! q
      s := (h-l)/(minPoints()::F-1)
      if (first t) ^= l then
        t := c := concat(l,c)
        p := q := concat(f l,p)
        NUMFUNEVALS := NUMFUNEVALS + 1
      while not null rest t repeat
        n := wholePart((second(t) - first(t))/s)
        d := (second(t) - first(t))/((n+1)::F)
        for i in 1..n repeat
          t.rest := concat(first(t) + d,rest t)
          p.rest := concat(f second t,rest p)
          NUMFUNEVALS := NUMFUNEVALS + 1
          t := rest t; p := rest p
        t := rest t
        p := rest p
      xRange := select(q,xCoord,min) .. select(q,xCoord,max)
      yRange := select(q,yCoord,min) .. select(q,yCoord,max)
      [ f, [nRange,xRange,yRange], c, q]

    adaptivePlot(curve,tRange,xRange,yRange,pixelfraction) ==
      xDiff := hi xRange - lo xRange
      yDiff := hi yRange - lo yRange
      xDiff = 0 or yDiff = 0 => curve
      l := lo tRange; h := hi tRange
      (tDiff := h-l) = 0 => curve
      t := curve.knots
      #t < 3 => curve
      p := curve.points; f := curve.source
      minLength:F := 4::F/500::F
      maxLength:F := 1::F/6::F
      tLimit := tDiff/(pixelfraction*500)::F
      while not null t and first t < l repeat (t := rest t; p := rest p)
      #t < 3 => curve
      headert := t; headerp := p
      st := t; sp := p
      todot : L L F := nil()
      todop : L L P := nil()
      while not null rest rest st repeat
        todot := concat_!(todot, st)
        todop := concat_!(todop, sp)
        st := rest st; sp := rest sp
      st := headert; sp := headerp
      todo1 := todot; todo2 := todop
      n : I := 0
      while not null todo1 repeat
        st := first(todo1)
        t0 := first(st); t1 := second(st); t2 := third(st)
        if t2 > h then leave
        t2 - t0 < tLimit =>
            todo1 := rest todo1
            todo2 := rest todo2
            if not null todo1 then (t := first(todo1); p := first(todo2))
        sp := first(todo2)
        x0 := xCoord first(sp); y0 := yCoord first(sp)
        x1 := xCoord second(sp); y1 := yCoord second(sp)
        x2 := xCoord third(sp); y2 := yCoord third(sp)
        a1 := (x1-x0)/xDiff; b1 := (y1-y0)/yDiff
        a2 := (x2-x1)/xDiff; b2 := (y2-y1)/yDiff
        s1 := sqrt(a1**2+b1**2); s2 := sqrt(a2**2+b2**2)
        dp := a1*a2+b1*b2

        s1 < maxLength and s2 < maxLength and _
          (s1 = 0::F or s2 = 0::F or
             s1 < minLength and s2 < minLength or _
             dp/s1/s2 > ANGLEBOUND) =>
                todo1 := rest todo1
                todo2 := rest todo2
                if not null todo1 then (t := first(todo1); p := first(todo2))
        if n > MAXPOINTS then leave else n := n + 1
        st := rest t
        if not null rest rest st then
          tm := (t0+t1)/2::F
          tj := tm
          t.rest := concat(tj,rest t)
          p.rest := concat(f tj, rest p)
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          t := rest t; p := rest p
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          t := rest t; p := rest p
          todo1 := rest todo1; todo2 := rest todo2

          tm := (t1+t2)/2::F
          tj := tm
          t.rest := concat(tj, rest t)
          p.rest := concat(f tj, rest p)
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          t := rest t; p := rest p
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          todo1 := rest todo1
          todo2 := rest todo2
          if not null todo1 then (t := first(todo1); p := first(todo2))
        else
          tm := (t0+t1)/2::F
          tj := tm
          t.rest := concat(tj,rest t)
          p.rest := concat(f tj, rest p)
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          t := rest t; p := rest p
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          t := rest t; p := rest p

          tm := (t1+t2)/2::F
          tj := tm
          t.rest := concat(tj, rest t)
          p.rest := concat(f tj, rest p)
          todo1 := concat_!(todo1, t)
          todo2 := concat_!(todo2, p)
          todo1 := rest todo1
          todo2 := rest todo2
          if not null todo1 then (t := first(todo1); p := first(todo2))
      n > 0 =>
        NUMFUNEVALS := NUMFUNEVALS + n
        t := curve.knots; p := curve.points
        xRange := select(p,xCoord,min) .. select(p,xCoord,max)
        yRange := select(p,yCoord,min) .. select(p,yCoord,max)
        [ curve.source, [tRange,xRange,yRange], t, p ]
      curve

    basicPlot(f,tRange) ==
      checkRange tRange
      l := lo tRange
      h := hi tRange
      t : L F := list l
      p : L P := list f l
      s := (h-l)/(minPoints()-1)::F
      for i in 2..minPoints()-1 repeat
        l := l+s 
        t := concat(l,t) 
        p := concat(f l,p)
      t := reverse_! concat(h,t)
      p := reverse_! concat(f h,p)
      xRange : R := select(p,xCoord,min) .. select(p,xCoord,max)
      yRange : R := select(p,yCoord,min) .. select(p,yCoord,max)
      [ f, [tRange,xRange,yRange], t, p ]

    zoom(p,xRange) ==
      [p.parametric, [xRange,third(p.display)], p.bounds, _
       p.axisLabels, p.functions]

    zoom(p,xRange,yRange) ==
      [p.parametric, [xRange,yRange], p.bounds, _
       p.axisLabels, p.functions]

    basicRefine(curve,nRange) ==
      tRange:R := first curve.ranges
      -- curve := copy$C curve  -- Yet another compiler bug
      curve: C := [curve.source,curve.ranges,curve.knots,curve.points]
      t := curve.knots := copy curve.knots
      p := curve.points := copy curve.points
      l := lo nRange; h := hi nRange
      f := curve.source
      while not null rest t and first t < h repeat
        second(t) < l => (t := rest t; p := rest p)
        -- insert new point between t.0 and t.1
        tm : F := (first(t) + second(t))/2::F
        pm := f tm
        NUMFUNEVALS := NUMFUNEVALS + 1
        t.rest := concat(tm,rest t); t := rest rest t
        p.rest := concat(pm,rest p); p := rest rest p
      t := curve.knots; p := curve.points
      xRange := select(p,xCoord,min) .. select(p,xCoord,max)
      yRange := select(p,yCoord,min) .. select(p,yCoord,max)
      [ curve.source, [tRange,xRange,yRange], t, p ]

    refine p == refine(p,parametricRange p)

    refine(p,nRange) ==
      NUMFUNEVALS := 0
      tRange := parametricRange p
      nRange := intersect(tRange,nRange)
      curves: L C := [basicRefine(c,nRange) for c in p.functions]
      xRange := join(curves,1); yRange := join(curves,2)
      if adaptive? then
        tlimit := if parametric? p then 8 else 1
        curves := [adaptivePlot(c,nRange,xRange,yRange, _
                   tlimit) for c in curves]
        xRange := join(curves,1); yRange := join(curves,2)
      [p.parametric, p.display, [tRange,xRange,yRange], _
       p.axisLabels, curves ]

    plot(p:%,tRange:R) ==
      -- re plot p on a new range making use of the points already
      -- computed if possible
      NUMFUNEVALS := 0
      curves: L C := [rangeRefine(c,tRange) for c in p.functions]
      xRange := join(curves,1); yRange := join(curves,2)
      if adaptive? then
        tlimit := if parametric? p then 8 else 1
        curves := [adaptivePlot(c,tRange,xRange,yRange,tlimit) for c in curves]
        xRange := join(curves,1); yRange := join(curves,2)
      [ p.parametric, [xRange,yRange], [tRange,xRange,yRange],
        p.axisLabels, curves ]

    pt(xx,yy) == point(l : L F := [xx,yy])

    myTrap: (F-> F, F) -> F
    myTrap(ff:F-> F, f:F):F ==
      s := trapNumericErrors(ff(f))$Lisp :: Union(F, "failed")
      s case "failed" => _$NaNvalue$Lisp
      r:F:=s::F
      r > max()$F or r < min()$F => _$NaNvalue$Lisp
      r

    plot(f:F -> F,xRange:R) ==
      p := basicPlot((u1:F):P +-> pt(u1,myTrap(f,u1)),xRange)
      r := p.ranges
      NUMFUNEVALS := minPoints()
      if adaptive? then
        p := adaptivePlot(p,first r,second r,third r,1)
        r := p.ranges
      [ false, rest r, r, nil(), [ p ] ]

    plot(f:F -> F,xRange:R,yRange:R) ==
      p := plot(f,xRange)
      p.display := [xRange,checkRange yRange]
      p

    plot(f:F -> F,g:F -> F,tRange:R) ==
      p := basicPlot((z1:F):P +-> pt(myTrap(f,z1),myTrap(g,z1)),tRange)
      r := p.ranges
      NUMFUNEVALS := minPoints()
      if adaptive? then
        p := adaptivePlot(p,first r,second r,third r,8)
        r := p.ranges
      [ true, rest r, r, nil(), [ p ] ]

    plot(f:F -> F,g:F -> F,tRange:R,xRange:R,yRange:R) ==
      p := plot(f,g,tRange)
      p.display := [checkRange xRange,checkRange yRange]
      p

    pointPlot(f:F -> P,tRange:R) ==
      p := basicPlot(f,tRange)
      r := p.ranges
      NUMFUNEVALS := minPoints()
      if adaptive? then
        p := adaptivePlot(p,first r,second r,third r,8)
        r := p.ranges
      [ true, rest r, r, nil(), [ p ] ]

    pointPlot(f:F -> P,tRange:R,xRange:R,yRange:R) ==
      p := pointPlot(f,tRange)
      p.display := [checkRange xRange,checkRange yRange]
      p

    plot(l:L(F -> F),xRange:R) ==
      if null l then error "empty list of functions"
      t: L C := 
        [ basicPlot((z1:F):P +-> pt(z1,myTrap(f,z1)),xRange) for f in l ]
      yRange := join(t,2)
      NUMFUNEVALS := # l * minPoints()
      if adaptive? then
        t := [adaptivePlot(p,xRange,xRange,yRange,1) _
                for f in l for p in t]
        yRange := join(t,2)
      [false, [xRange,yRange], [xRange,xRange,yRange], nil(), t ]

    plot(l:L(F -> F),xRange:R,yRange:R) ==
      p := plot(l,xRange)
      p.display := [xRange,checkRange yRange]
      p

    plotPolar(f,thetaRange) ==
      plot((u1:F):F +-> f(u1) * cos(u1),
           (v1:F):F +-> f(v1) * sin(v1),thetaRange)

    plotPolar f == plotPolar(f,segment(0,2*pi()))

--% terminal output

    coerce r ==
      spaces: OUT := coerce "   "
      xSymbol := "x = " :: OUT
      ySymbol := "y = " :: OUT
      tSymbol := "t = " :: OUT
      plotSymbol := "PLOT" :: OUT
      tRange := (parametricRange r) :: OUT
      f : L OUT := nil()
      for curve in r.functions repeat
        xRange := second(curve.ranges) :: OUT
        yRange := third(curve.ranges) :: OUT
        l : L OUT := [xSymbol,xRange,spaces,ySymbol,yRange]
        if parametric? r then
          l := concat_!([tSymbol,tRange,spaces],l)
        h : OUT := hconcat l
        l := [p::OUT for p in curve.points]
        f := concat(vconcat concat(h,l),f)
      prefix("PLOT" :: OUT, reverse_! f)