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

)abbrev domain PLOT3D Plot3D
++ Author: Clifton J. Williamson based on code by Michael Monagan
++ Date Created: Jan 1989
++ Date Last Updated: 22 November 1990 (Jon Steinbach)
++ Description:
++ Plot3D supports parametric plots 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 are
++ real number systems. The facilities at this point are limited
++ to 3-dimensional parametric plots.

Plot3D() : SIG == CODE where

  B   ==> Boolean
  F   ==> DoubleFloat
  I   ==> Integer
  L   ==> List
  N   ==> NonNegativeInteger
  OUT ==> OutputForm
  P   ==> Point F
  S   ==> String
  R   ==> Segment F
  O   ==> OutputPackage
  C   ==> Record(source: F -> P,ranges: L R, knots: L F, points: L P)
 
  SIG ==> PlottableSpaceCurveCategory with
 
    pointPlot : (F -> P,R) -> %
      ++ pointPlot(f,g,h,a..b) plots {/emx = f(t), y = g(t), z = h(t)} as
      ++ t ranges over {/em[a,b]}.

    pointPlot : (F -> P,R,R,R,R) -> %
        ++ pointPlot(f,x,y,z,w) is not documented 

    plot : (F -> F,F -> F,F -> F,F -> F,R) -> %
      ++ plot(f,g,h,a..b) plots {/emx = f(t), y = g(t), z = h(t)} as
      ++ t ranges over {/em[a,b]}.

    plot : (F -> F,F -> F,F -> F,F -> F,R,R,R,R) -> %
        ++ plot(f1,f2,f3,f4,x,y,z,w) is not documented
 
    plot : (%,R) -> %                       -- change the range
         ++ plot(x,r) is not documented

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

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

    refine : % -> %
        ++ refine(x) is not documented
 
    tRange : % -> R
      ++ tRange(p) returns the range of the parameter in a parametric plot p.

    tValues : % -> L L F
      ++ tValues(p) returns a list of lists of the values of the parameter for
      ++ which a point is computed, one list for each curve in the plot p.
 
    minPoints3D : () -> I
      ++ minPoints3D() returns the minimum number of points in a plot.

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

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

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

    screenResolution3D : () -> I
      ++ screenResolution3D() returns the screen resolution for a 3d graph.

    setScreenResolution3D : I -> I
      ++ setScreenResolution3D(i) sets the screen resolution for a 3d graph
      ++ to i.

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

    setAdaptive3D : B -> B
      ++ setAdaptive3D(true) turns adaptive plotting on;
      ++ setAdaptive3D(false) turns adaptive plotting off.

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

    debug3D : B -> B
      ++ debug3D(true) turns debug mode on;
      ++ debug3D(false) turns debug mode off.
 
  CODE ==> add

    import PointPackage(F)
 
    --% local functions
 
    fourth         : L R -> R
    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,R,I,I) -> C
    basicPlot      : (F -> P,R) -> C
    basicRefine    : (C,R) -> C
    point          : (F,F,F,F) -> P
 
--% representation
 
    Rep := Record( display: L R, _
                   bounds: L R, _
                   screenres: I, _
                   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
 
    point(xx,yy,zz,col) == point(l : L F := [xx,yy,zz,col])
 
    fourth list == first rest rest rest list
 
    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:R,t:R) == 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)
        i = 2 => third(rr.ranges)
        fourth(rr.ranges)
      for r in rest l repeat
        i = 0 => union(u,first(r.ranges))
        i = 1 => union(u,second(r.ranges))
        i = 2 => union(u,third(r.ranges))
        union(u,fourth(r.ranges))
      u

    parametricRange r == first(r.bounds)
 
    minPoints3D() == MINPOINTS

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

    maxPoints3D() == MAXPOINTS

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

    screenResolution3D() == SCREENRES

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

    adaptive3D?() == ADAPTIVE

    setAdaptive3D b == ADAPTIVE := b
 
    numFunEvals3D() == NUMFUNEVALS

    debug3D b == DEBUG := b
 
    xRange plot == second plot.bounds

    yRange plot == third plot.bounds

    zRange plot == fourth plot.bounds

    tRange plot == first plot.bounds
 
    tValues plot ==
      outList : L L F := nil()
      for curve in plot.functions repeat
        outList := concat(curve.knots,outList)
      outList
 
    select(l,f,g) ==
      m := f first l
      if (EQL(m, _$NaNvalue$Lisp)$Lisp) then m := 0
      for p in rest l repeat
        fp : F := f p
        if (EQL(fp, _$NaNvalue$Lisp)$Lisp) then fp := 0
        m := g(m,fp)
      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); t1 := second t
          p.rest := concat(f t1,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)
      zRange := select(q,zCoord,min) .. select(q,zCoord,max)
      [f,[nRange,xRange,yRange,zRange],c,q]
 

    adaptivePlot(curve,tRg,xRg,yRg,zRg,pixelfraction,resolution) ==
      xDiff := hi xRg - lo xRg
      yDiff := hi yRg - lo yRg
      zDiff := hi zRg - lo zRg
      if xDiff = 0::F then xDiff := 1::F
      if yDiff = 0::F then yDiff := 1::F
      if zDiff = 0::F then zDiff := 1::F
      l := lo tRg; h := hi tRg
      (tDiff := h-l) = 0 => curve
      t := curve.knots
      #t < 3 => curve
      p := curve.points; f := curve.source
      minLength:F := 4::F/resolution::F
      maxLength := 1/4::F
      tLimit := tDiff/(pixelfraction*resolution)::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); 
        z0 := zCoord first(sp) 
        x1 := xCoord second(sp); 
        y1 := yCoord second(sp); 
        z1 := zCoord second(sp) 
        x2 := xCoord third(sp); y2 := yCoord third(sp); z2 := zCoord third(sp)
        a1 := (x1-x0)/xDiff; b1 := (y1-y0)/yDiff; c1 := (z1-z0)/zDiff
        a2 := (x2-x1)/xDiff; b2 := (y2-y1)/yDiff; c2 := (z2-z1)/zDiff
        s1 := sqrt(a1**2+b1**2+c1**2); s2 := sqrt(a2**2+b2**2+c2**2)
        dp := a1*a2+b1*b2+c1*c2
        s1 < maxLength and s2 < maxLength and _
           (s1 = 0 or s2 = 0 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
        --if DEBUG then
           --r : L F := [minLength,maxLength,s1,s2,dp/s1/s2,ANGLEBOUND]
           --output(r::E)$O
        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))
      if n > 0 then
        NUMFUNEVALS := NUMFUNEVALS + n
        t := curve.knots; p := curve.points
        xRg := select(p,xCoord,min) .. select(p,xCoord,max)
        yRg := select(p,yCoord,min) .. select(p,yCoord,max)
        zRg := select(p,zCoord,min) .. select(p,zCoord,max)
        [curve.source,[tRg,xRg,yRg,zRg],t,p]
      else 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)
      zRange : R := select(p,zCoord,min) .. select(p,zCoord,max)
      [f,[tRange,xRange,yRange,zRange],t,p]
 
    zoom(p,xRange,yRange,zRange) ==
     [[xRange,yRange,zRange],p.bounds,
      p.screenres,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
        -- if DEBUG then output$O (tm::E)
        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)
      zRange := select(p,zCoord,min) .. select(p,zCoord,max)
      [curve.source,[tRange,xRange,yRange,zRange],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)
      zRange := join(curves,3)
      scrres := p.screenres
      if adaptive3D? then
        tlimit := 8
        curves := [adaptivePlot(c,nRange,xRange,yRange,zRange, _
                   tlimit,scrres := 2*scrres) for c in curves]
        xRange := join(curves,1); yRange := join(curves,2)
        zRange := join(curves,3)
      [p.display,[tRange,xRange,yRange,zRange], _
       scrres,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)
      zRange := join(curves,3)
      if adaptive3D? then
        tlimit := 8
        curves := [adaptivePlot(c,tRange,xRange,yRange,zRange,tlimit, _
                      p.screenres) for c in curves]
        xRange := join(curves,1); yRange := join(curves,2)
        zRange := join(curves,3)
      [[xRange,yRange,zRange],[tRange,xRange,yRange,zRange],
        p.screenres,p.axisLabels,curves]
 
    pointPlot(f:F -> P,tRange:R) ==
      p := basicPlot(f,tRange)
      r := p.ranges
      NUMFUNEVALS := MINPOINTS
      if adaptive3D? then
       p := adaptivePlot(p,first r,second r,third r,fourth r,8,SCREENRES)
      [ rest r, r, SCREENRES, nil(), [ p ] ]
 
    pointPlot(f:F -> P,tRange:R,xRange:R,yRange:R,zRange:R) ==
      p := pointPlot(f,tRange)
      p.display:= [checkRange xRange,checkRange yRange,checkRange zRange]
      p
 
    myTrap: (F-> F, F) -> F
    myTrap(ff:F-> F, f:F):F ==
      s := trapNumericErrors(ff(f))$Lisp :: Union(F, "failed")
      if (s) case "failed" then
        r:F := _$NaNvalue$Lisp
      else
        r:F := s
      r

    plot(f1:F -> F,f2:F -> F,f3:F -> F,col:F -> F,tRange:R) ==
      p := basicPlot(
       (z:F):P+->point(myTrap(f1,z),myTrap(f2,z),myTrap(f3,z),col(z)),tRange)
      r := p.ranges
      NUMFUNEVALS := MINPOINTS
      if adaptive3D? then
       p := adaptivePlot(p,first r,second r,third r,fourth r,8,SCREENRES)
      [ rest r, r, SCREENRES, nil(), [ p ] ]
 
    plot(f1:F -> F,f2:F -> F,f3:F -> F,col:F -> F,_
              tRange:R,xRange:R,yRange:R,zRange:R) ==
      p := plot(f1,f2,f3,col,tRange)
      p.display:= [checkRange xRange,checkRange yRange,checkRange zRange]
      p
 
--% terminal output
 
    coerce r ==
      spaces := "   " :: OUT
      xSymbol := "x = " :: OUT; ySymbol := "y = " :: OUT
      zSymbol := "z = " :: OUT; tSymbol := "t = " :: OUT
      tRange := (parametricRange r) :: OUT
      f : L OUT := nil()
      for curve in r.functions repeat
        xRange := coerce curve.ranges.1
        yRange := coerce curve.ranges.2
        zRange := coerce curve.ranges.3
        l : L OUT := [xSymbol,xRange,spaces,ySymbol,yRange,_
                       spaces,zSymbol,zRange]
        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)
 
----% graphics output
 
    listBranches plot ==
      outList : L L P := nil()
      for curve in plot.functions repeat
        outList := concat(curve.points,outList)
      outList