asymptote 1.83

git-svn-id: svn://tug.org/texlive/trunk@14696 c570f23f-e606-0410-a88d-b1316a301751

32 files changed, 1883 insertions, 1182 deletions

diff --git a/Master/texmf/asymptote/GUI/xasyVersion.py b/Master/texmf/asymptote/GUI/xasyVersion.py
index a1ab6b203c5..5d5a51724f2 100755
--- a/Master/texmf/asymptote/GUI/xasyVersion.py
+++ b/Master/texmf/asymptote/GUI/xasyVersion.py
@@ -1,2 +1,2 @@
 #!/usr/bin/env python
-xasyVersion = "1.82"
+xasyVersion = "1.83"
diff --git a/Master/texmf/asymptote/animation.asy b/Master/texmf/asymptote/animation.asy
index 6fcd2ad70e2..27f7171a639 100644
--- a/Master/texmf/asymptote/animation.asy
+++ b/Master/texmf/asymptote/animation.asy
@@ -8,16 +8,17 @@
 // animation delay is in milliseconds
 real animationdelay=50;
 
-typedef frame fit(string prefix="",picture);
+typedef frame enclosure(frame);
 
-frame NoBox(string prefix="", picture pic) {
-  return pic.fit(prefix);
+frame NoBox(frame f) {
+  return f;
 }
 
-fit BBox(string prefix="", real xmargin=0, real ymargin=xmargin,
-         pen p=currentpen, filltype filltype=NoFill) {
-  return new frame(string prefix, picture pic) {
-    return bbox(prefix,pic,xmargin,ymargin,p,filltype);
+enclosure BBox(real xmargin=0, real ymargin=xmargin,
+               pen p=currentpen, filltype filltype=NoFill) {
+  return new frame(frame f) {
+    box(f,xmargin,ymargin,p,filltype,above=false);
+    return f;
   };
 }
 
@@ -59,11 +60,11 @@ struct animation {
     shipped=false;
   }
   
-  void add(picture pic=currentpicture, fit fit=NoBox) {
+  void add(picture pic=currentpicture, enclosure enclosure=NoBox) {
     if(global) {
       ++index;
       pictures.push(pic.copy());
-    } else this.shipout(fit(prefix,pic));
+    } else this.shipout(enclosure(pic.fit()));
   }
   
   void purge(bool keep=settings.keep) {
@@ -86,30 +87,33 @@ struct animation {
     return rc;
   }
 
+  void glmovie(string prefix=prefix, projection P=currentprojection) {
+    if(!view() || settings.render == 0) return;
+    fit(prefix,pictures,view=true,P);
+  }
+
   // Export all frames with the same scaling.
-  void export(string prefix=prefix, fit fit=NoBox,
-              bool multipage=false, bool view=false) {
+  void export(string prefix=prefix, enclosure enclosure=NoBox,
+              bool multipage=false, bool view=false,
+              projection P=currentprojection) {
     if(pictures.length == 0) return;
     if(!global) multipage=false;
-    rescale(pictures);
-    frame multi;
     bool inlinetex=settings.inlinetex;
     if(multipage)
       settings.inlinetex=false;
-    for(int i=0; i < pictures.length; ++i) {
+    frame multi;
+    frame[] fits=fit(prefix,pictures,view=false,P);
+    for(int i=0; i < fits.length; ++i) {
+      string s=name(prefix,i);
       if(multipage) {
-        add(multi,fit(prefix,pictures[i]));
+        add(multi,enclosure(fits[i]));
         newpage(multi);
+        files.push(s+"."+nativeformat());
       } else {
-        if(pictures[i].empty3() || settings.render <= 0) {
-          real render=settings.render;
-          settings.render=0;
-          this.shipout(name(prefix,i),fit(prefix,pictures[i]));
-          settings.render=render;
-        } else { // Render 3D frames
-          files.push(name(prefix,i)+"."+nativeformat());
-          fit(prefix,pictures[i]);
-        }
+        if(pictures[i].empty3() || settings.render <= 0)
+          this.shipout(s,enclosure(fits[i]));
+        else // 3D frames
+          files.push(s+"."+nativeformat());
       }
     }
     if(multipage) {
@@ -129,8 +133,8 @@ struct animation {
     return s;
   }
 
-  string pdf(fit fit=NoBox, real delay=animationdelay, string options="",
-             bool keep=settings.keep, bool multipage=true) {
+  string pdf(enclosure enclosure=NoBox, real delay=animationdelay,
+             string options="", bool keep=settings.keep, bool multipage=true) {
     if(settings.inlinetex) multipage=true;
     if(!global) multipage=false;
     if(settings.tex != "pdflatex")
@@ -141,16 +145,17 @@ struct animation {
     string pdfname=filename+".pdf";
 
     if(global)
-      export(filename,fit,multipage=multipage);
+      export(filename,enclosure,multipage=multipage);
     
     shipped=false;
 
-    if(!keep && !settings.inlinetex) {
+    if(!keep) {
       exitfcn currentexitfunction=atexit();
       void exitfunction() {
         if(currentexitfunction != null) currentexitfunction();
-        this.purge();
-        if(multipage)
+        if(multipage || !settings.inlinetex)
+          this.purge();
+        if(multipage && !settings.inlinetex)
           delete(pdfname);
       }
       atexit(exitfunction);
@@ -162,15 +167,15 @@ struct animation {
     return load(index,delay,options,multipage);
   }
 
-  int movie(fit fit=NoBox, int loops=0, real delay=animationdelay,
+  int movie(enclosure enclosure=NoBox, int loops=0, real delay=animationdelay,
             string format=settings.outformat == "" ? "gif" : settings.outformat,
             string options="", bool keep=settings.keep) {
     if(global) {
       if(format == "pdf") {
-        export(fit,multipage=true,view=true);
+        export(enclosure,multipage=true,view=true);
         return 0;
       }
-      export(fit);
+      export(enclosure);
     }
     return merge(loops,delay,format,options,keep);
   }
diff --git a/Master/texmf/asymptote/bezulate.asy b/Master/texmf/asymptote/bezulate.asy
index 2b852623faf..d935d358989 100644
--- a/Master/texmf/asymptote/bezulate.asy
+++ b/Master/texmf/asymptote/bezulate.asy
@@ -238,7 +238,7 @@ path[] bezulate(path[] p)
             // avoid infinite recursion
             ++refinements;
             if(refinements > mantissaBits) {
-              write("warning: too many subdivisions");
+              warning("subdivisions","too many subdivisions",position=true);
             } else {
               p=subdivide(p);
               i=-1;
diff --git a/Master/texmf/asymptote/geometry.asy b/Master/texmf/asymptote/geometry.asy
index 12e35339dcc..0db3172d365 100644
--- a/Master/texmf/asymptote/geometry.asy
+++ b/Master/texmf/asymptote/geometry.asy
@@ -228,24 +228,6 @@ pair[] intersectionpoints(pair A, pair B, real[] equation)
 
 // *=======================================================*
 // *......................COORDINATES......................*
-// Copyright (c) 2007, Philippe Ivaldi.
-// Version: $Id: coordinates.asy,v 0.0 2007/02/03 16:06:23 Philippe Ivaldi Exp$
-// Last modified: Wed Aug 15 15:53:01 CEST 2007
-
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 3 of the License, or
-// any later version.
-
-// This program is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-// General Public License for more details.
-
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-// 02110-1301, USA.
 
 real EPS=sqrt(realEpsilon);
 
@@ -492,8 +474,9 @@ bool samecoordsys(bool warn=true ... point[] M)
     t=M[i].coordsys;
   }
   if(warn && !ret)
-    write("Warning, the coordinate system of two objects are not the same.
-The operation will be done relatively to the default coordinate system.");
+    warning("coodinatesystem",
+            "the coordinate system of two objects are not the same.
+The operation will be done relative to the default coordinate system.");
   return ret;
 }
 
@@ -2418,8 +2401,10 @@ bool samecoordsys(bool warn=true ... bqe[] bqes)
     t=bqes[i].coordsys;
   }
   if(warn && !ret)
-    write("Warning, the coordinate system of two  bivariate quadratic equations are not the same.
-The operation will be done relatively to the default coordinate system.");
+    warning("coodinatesystem",
+            "the coordinate system of two bivariate quadratic equations are not
+the same. The operation will be done relatively to the default coordinate
+system.");
   return ret;
 }
 
@@ -2710,7 +2695,9 @@ int circlenodesnumberfactor=100;/*<asyxml></code><documentation>Factor for the n
 /*<asyxml><function type="int" signature="circlenodesnumber(real)"><code></asyxml>*/
 int circlenodesnumber(real r)
 {/*<asyxml></code><documentation>Return the number of nodes for drawing a circle of radius 'r'.</documentation></function></asyxml>*/
-  if (circlenodesnumberfactor < 100) write("Warning: variable 'circlenodesnumberfactor' maybe too small.");
+  if (circlenodesnumberfactor < 100)
+    warning("circlenodesnumberfactor",
+            "variable 'circlenodesnumberfactor' may be too small.");
   int oi=ceil(circlenodesnumberfactor*abs(r)^0.1);
   oi=45*floor(oi/45);
   return oi == 0 ? 4 : conicnodesfactor*oi;
@@ -2729,7 +2716,9 @@ int ellipsenodesnumberfactor=250;/*<asyxml></code><documentation>Factor for the
 /*<asyxml><function type="int" signature="ellipsenodesnumber(real,real)"><code></asyxml>*/
 int ellipsenodesnumber(real a, real b)
 {/*<asyxml></code><documentation>Return the number of nodes to draw a ellipse of axis 'a' and 'b'.</documentation></function></asyxml>*/
-  if (ellipsenodesnumberfactor < 250) write("Warning: variable 'ellipsenodesnumberfactor' maybe too small.");
+  if (ellipsenodesnumberfactor < 250)
+    write("ellipsenodesnumberfactor",
+          "variable 'ellipsenodesnumberfactor' maybe too small.");
   int tmp=circlenodesnumberfactor;
   circlenodesnumberfactor=ellipsenodesnumberfactor;
   int oi=circlenodesnumber(max(abs(a),abs(b))/min(abs(a),abs(b)));
@@ -5481,7 +5470,7 @@ circle excircle(point A, point B, point C)
 }
 
 private int[] numarray={1,2,3};
-numarray.cyclic(true);
+numarray.cyclic=true;
 
 /*<asyxml><struct signature="triangle"><code></asyxml>*/
 struct triangle {/*<asyxml></code><documentation></documentation></asyxml>*/
@@ -6440,13 +6429,19 @@ point operator *(inversion i, point P)
   return inverse(i.k,i.C,P);
 }
 
+void lineinversion()
+{
+    warning("lineinversion","the inversion of the line is not a circle.
+The returned circle has an infinite radius, circle.l has been set.");
+}
+
+
 /*<asyxml><function type="circle" signature="inverse(real,point,line)"><code></asyxml>*/
 circle inverse(real k, point A, line l)
 {/*<asyxml></code><documentation>Return the inverse circle of 'l' with
    respect to point 'A' and inversion radius 'k'.</documentation></function></asyxml>*/
   if(A @ l) {
-    write("Warning: the inversion of the line is not a circle.");
-    write("The returned circle has an infinite radius, cirlce.l have been set.");
+    lineinversion();
     circle C=circle(A, infinity);
     C.l=l;
     return C;
@@ -6467,8 +6462,7 @@ circle inverse(real k, point A, circle c)
    respect to point A and inversion radius 'k'.</documentation></function></asyxml>*/
   if(degenerate(c)) return inverse(k,A,c.l);
   if(A @ c) {
-    write("Warning: the inversion of the circle is not a circle.");
-    write("The returned circle has an infinite radius, cirlce.l have been set.");
+    lineinversion();
     point M=rotate(180,c.C)*A, Mp=rotate(90,c.C)*A;
     circle oc=circle(A,infinity);
     oc.l=line(inverse(k,A,M),inverse(k,A,Mp));
@@ -7123,7 +7117,7 @@ path square(pair z1, pair z2)
 // relative to the path z--z+dir.
 void perpendicular(picture pic=currentpicture, pair z, pair align,
                    pair dir=E, real size=0, pen p=currentpen,
-                   margin margin=NoMargin, filltype filltype=NoFill) 
+                   margin margin=NoMargin, filltype filltype=NoFill)
 {
   perpendicularmark(pic,(point) z,align,dir,size,p,margin,filltype);
 }
@@ -7133,7 +7127,7 @@ void perpendicular(picture pic=currentpicture, pair z, pair align,
 // relative to the path z--z+dir(g,0)
 void perpendicular(picture pic=currentpicture, pair z, pair align, path g,
                    real size=0, pen p=currentpen, margin margin=NoMargin,
-                   filltype filltype=NoFill) 
+                   filltype filltype=NoFill)
 {
   perpendicularmark(pic,(point) z,align,dir(g,0),size,p,margin,filltype);
 }
diff --git a/Master/texmf/asymptote/graph.asy b/Master/texmf/asymptote/graph.asy
index 2983137b612..55974c63f7d 100644
--- a/Master/texmf/asymptote/graph.asy
+++ b/Master/texmf/asymptote/graph.asy
@@ -1,6 +1,6 @@
 private import math;
 
-import splinetype;
+import graph_splinetype;
 import graph_settings;
 
 scaleT Linear;
@@ -1893,7 +1893,7 @@ interpolate Hermite(splinetype splinetype)
   };
 }
 
-interpolate Hermite=Hermite(defaultspline);
+interpolate Hermite=Hermite(Spline);
 
 guide graph(picture pic=currentpicture, real f(real), real a, real b,
             int n=ngraph, real T(real)=identity, interpolate join=operator --)
diff --git a/Master/texmf/asymptote/graph3.asy b/Master/texmf/asymptote/graph3.asy
index 8fe6c9ed5c6..c2630e10f11 100644
--- a/Master/texmf/asymptote/graph3.asy
+++ b/Master/texmf/asymptote/graph3.asy
@@ -1509,6 +1509,27 @@ path3[] segment(triple[] v, bool[] cond, interpolate3 join=operator --)
                   segment.length);
 }
 
+bool uperiodic(triple[][] a) {
+  int n=a.length;
+  if(n == 0) return false;
+  int m=a[0].length;
+  triple[] a0=a[0];
+  triple[] a1=a[n-1];
+  real epsilon=sqrtEpsilon*norm(a);
+  for(int j=0; j < m; ++j)
+    if(abs(a0[j]-a1[j]) > epsilon) return false;
+  return true;
+}
+bool vperiodic(triple[][] a) {
+  int n=a.length;
+  if(n == 0) return false;
+  int m=a[0].length-1;
+  real epsilon=sqrtEpsilon*norm(a);
+  for(int i=0; i < n; ++i)
+    if(abs(a[i][0]-a[i][m]) > epsilon) return false;
+  return true;
+}
+
 // return the surface described by a matrix f
 surface surface(triple[][] f, bool[][] cond={})
 {
@@ -1533,6 +1554,7 @@ surface surface(triple[][] f, bool[][] cond={})
   }
 
   surface s=surface(count);
+  s.index=new int[nx][ny];
   int k=-1;
   for(int i=0; i < nx; ++i) {
     bool[] condi,condp;
@@ -1542,122 +1564,18 @@ surface surface(triple[][] f, bool[][] cond={})
     }
     triple[] fi=f[i];
     triple[] fp=f[i+1];
+    int[] indexi=s.index[i];
     for(int j=0; j < ny; ++j) {
       if(all || (condi[j] && condi[j+1] && condp[j] && condp[j+1]))
         s.s[++k]=patch(new triple[] {fi[j],fp[j],fp[j+1],fi[j+1]});
-    }
-  }
-  return s;
-}
-
-private surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
-                         real[] x, real[] y, bool[][] cond={})
-{ // z[i][j] is the value at (x[i],y[j])
-  // p and q are the first derivatives with respect to x and y, respectively
-  // r is the second derivative ddu/dxdy
-  int n=x.length-1;
-  int m=y.length-1;
-
-  bool all=cond.length == 0;
-
-  int count;
-  if(all)
-    count=n*m;
-  else {
-    count=0;
-    for(int i=0; i < n; ++i) {
-      bool[] condi=cond[i];
-      for(int j=0; j < m; ++j)
-        if(condi[j]) ++count;
-    }
-  }
-
-  surface g=surface(count);
-  int k=-1;
-  for(int i=0; i < n; ++i) {
-    bool[] condi=all ? null : cond[i];
-    real xi=x[i];
-    real[] zi=z[i];
-    real[] zp=z[i+1];
-    real[] ri=r[i];
-    real[] rp=r[i+1];
-    real[] pi=p[i];
-    real[] pp=p[i+1];
-    real[] qi=q[i];
-    real[] qp=q[i+1];
-    real xp=x[i+1];
-    real hx=(xp-xi)/3;
-    for(int j=0; j < m; ++j) {
-      real yj=y[j];
-      real yp=y[j+1];
-      if(all || condi[j]) {
-        triple[][] P=array(4,array(4,O));
-        real hy=(yp-yj)/3;
-        real hxy=hx*hy;
-        // first x and y  directions
-        for(int k=0; k < 4; ++k) {
-          P[k][0] += xi*X;
-          P[0][k] += yj*Y;
-          P[k][1] += (xp+2*xi)/3*X;
-          P[1][k] += (yp+2*yj)/3*Y;
-          P[k][2] += (2*xp+xi)/3*X;
-          P[2][k] += (2*yp+yj)/3*Y;
-          P[k][3] += xp*X;
-          P[3][k] += yp*Y;
-        }
-        // now z, first the value 
-        P[0][0] += zi[j]*Z;
-        P[0][3] += zp[j]*Z;
-        P[3][0] += zi[j+1]*Z;
-        P[3][3] += zp[j+1]*Z;
-        // 2nd, first derivative
-        P[0][1] += (P[0][0].z+hx*pi[j])*Z;
-        P[3][1] += (P[3][0].z+hx*pi[j+1])*Z;
-        P[0][2] += (P[0][3].z-hx*pp[j])*Z;
-        P[3][2] += (P[3][3].z-hx*pp[j+1])*Z;
-        P[1][0] += (P[0][0].z+hy*qi[j])*Z;
-        P[1][3] += (P[0][3].z+hy*qp[j])*Z;
-        P[2][0] += (P[3][0].z-hy*qi[j+1])*Z;
-        P[2][3] += (P[3][3].z-hy*qp[j+1])*Z;
-        // 3nd, second derivative
-        P[1][1] += (P[1][0].z+P[0][1].z-P[0][0].z+hxy*ri[j])*Z;
-        P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*ri[j+1])*Z;
-        P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*rp[j])*Z;
-        P[2][2] += (P[3][2].z+P[2][3].z-P[3][3].z+hxy*rp[j+1])*Z;
-        g.s[++k]=patch(P);
+        indexi[j]=k;
       }
-    }
   }
-  return g;
-}
 
-// return the surface described by a real matrix f, interpolated with
-// xsplinetype and ysplinetype.
-surface surface(real[][] f, real[] x, real[] y,
-                splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype,
-                bool[][] cond={})
-{
-  real epsilon=sqrtEpsilon*max(abs(y));
-  if(xsplinetype == null)
-    xsplinetype=(abs(x[0]-x[x.length-1]) <= epsilon) ? periodic : notaknot;
-  if(ysplinetype == null)
-    ysplinetype=(abs(y[0]-y[y.length-1]) <= epsilon) ? periodic : notaknot;
-  int n=x.length; int m=y.length;
-  real[][] ft=transpose(f);
-  real[][] tp=new real[m][];
-  for(int j=0; j < m; ++j)
-    tp[j]=xsplinetype(x,ft[j]);
-  real[][] q=new real[n][];
-  for(int i=0; i < n; ++i)
-    q[i]=ysplinetype(y,f[i]);
-  real[][] qt=transpose(q);
-  real[] d1=xsplinetype(x,qt[0]);
-  real[] d2=xsplinetype(x,qt[m-1]);
-  real[][] r=new real[n][];
-  real[][] p=transpose(tp);
-  for(int i=0; i < n; ++i)
-    r[i]=clamped(d1[i],d2[i])(y,p[i]);
-  return bispline(f,p,q,r,x,y,cond);
+  if(uperiodic(f)) s.ucyclic(true);
+  if(vperiodic(f)) s.vcyclic(true);
+  
+  return s;
 }
 
 // return the surface described by a real matrix f, interpolated with
@@ -1738,74 +1656,8 @@ surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu,
                 splinetype[] usplinetype, splinetype[] vsplinetype=Spline,
                 bool cond(pair z)=null)
 {
-  real[] upt=uniform(a.x,b.x,nu);
-  real[] vpt=uniform(a.y,b.y,nv);
-  real[] ipt=sequence(nu+1);
-  real[] jpt=sequence(nv+1);
-  real[][] fx=new real[nu+1][nv+1];
-  real[][] fy=new real[nu+1][nv+1];
-  real[][] fz=new real[nu+1][nv+1];
-
-  bool[][] active;
-  bool all=cond == null;
-  if(!all) active=new bool[nu+1][nv+1];
-
-  real norm;
-  for(int i=0; i <= nu; ++i) {
-    real upti=upt[i];
-    real[] fxi=fx[i];
-    real[] fyi=fy[i];
-    real[] fzi=fz[i];
-    bool[] activei=all ? null : active[i];
-    for(int j=0; j <= nv; ++j) {
-      pair z=(upti,vpt[j]);
-      triple f=(all || (activei[j]=cond(z))) ? f(z) : O;
-      norm=max(norm,abs(f));
-      fxi[j]=f.x;
-      fyi[j]=f.y;
-      fzi[j]=f.z;
-    }
-  }
-
-  real epsilon=sqrtEpsilon*norm;
-
-  if(usplinetype.length == 0) {
-    bool uperiodic(real[][] a) {
-      for(int j=0; j < nv; ++j)
-        if(abs(a[0][j]-a[nu][j]) > epsilon) return false;
-      return true;
-    }
-    usplinetype=new splinetype[] {uperiodic(fx) ? periodic : notaknot,
-                                  uperiodic(fy) ? periodic : notaknot,
-                                  uperiodic(fz) ? periodic : notaknot};
-  } else if(usplinetype.length != 3) abort("usplinetype must have length 3");
-
-  if(vsplinetype.length == 0) {
-    bool vperiodic(real[][] a) {
-      for(int i=0; i < nu; ++i)
-        if(abs(a[i][0]-a[i][nv]) > epsilon) return false;
-      return true;
-    }
-    vsplinetype=new splinetype[] {vperiodic(fx) ? periodic : notaknot,
-                                  vperiodic(fy) ? periodic : notaknot,
-                                  vperiodic(fz) ? periodic : notaknot};
-  } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3");
-  
-  surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
-  surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
-  surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
-
-  surface s=surface(sx.s.length);
-  for(int k=0; k < sx.s.length; ++k) {
-    triple[][] Q=new triple[4][4];
-    for(int i=0; i < 4 ; ++i) {
-      for(int j=0; j < 4 ; ++j) {
-        Q[i][j]=(sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z);
-      }
-      s.s[k]=patch(Q);
-    }
-  }
-  return s;
+  return surface(f,uniform(a.x,b.x,nu),uniform(a.y,b.y,nv),
+                 usplinetype,vsplinetype,cond);
 }
 
 // return the surface described by a real function f over box(a,b),
diff --git a/Master/texmf/asymptote/graph_splinetype.asy b/Master/texmf/asymptote/graph_splinetype.asy
new file mode 100644
index 00000000000..fbcb1b00c75
--- /dev/null
+++ b/Master/texmf/asymptote/graph_splinetype.asy
@@ -0,0 +1,251 @@
+typedef real[] splinetype(real[], real[]);
+
+restricted real[] Spline(real[] x, real[] y);
+restricted splinetype[] Spline;
+
+string morepoints="interpolation requires at least 2 points";
+string differentlengths="arrays have different lengths";
+void checklengths(int x, int y, string text=differentlengths)
+{
+  if(x != y)
+    abort(text+": "+string(x)+" != "+string(y));
+}
+
+// Linear interpolation
+real[] linear(real[] x, real[] y)
+{
+  int n=x.length;
+  checklengths(n,y.length);
+  real[] d=new real[n];
+  for(int i=0; i < n-1; ++i)
+    d[i]=(y[i+1]-y[i])/(x[i+1]-x[i]);
+  d[n-1]=d[n-2];
+  return d;
+}
+
+// Standard cubic spline interpolation with not-a-knot condition:
+// s'''(x_2^-)=s'''(x_2^+) et s'''(x_(n_2)^-)=s'''(x_(n-2)^+)
+// if n=2, linear interpolation is returned
+// if n=3, an interpolation polynomial of degree <= 2 is returned:
+// p(x_1)=y_1, p(x_2)=y_2, p(x_3)=y_3
+real[] notaknot(real[] x, real[] y)
+{
+  int n=x.length;
+  checklengths(n,y.length);
+  real[] d;
+  if(n > 3) {
+    real[] a=new real[n];
+    real[] b=new real[n];
+    real[] c=new real[n];
+    real[] g=new real[n];
+    b[0]=x[2]-x[1];
+    c[0]=x[2]-x[0];
+    a[0]=0;
+    g[0]=((x[1]-x[0])^2*(y[2]-y[1])/b[0]+b[0]*(2*b[0]+3*(x[1]-x[0]))*
+          (y[1]-y[0])/(x[1]-x[0]))/c[0];
+    for(int i=1; i < n-1; ++i) {
+      a[i]=x[i+1]-x[i];
+      c[i]=x[i]-x[i-1];
+      b[i]=2*(a[i]+c[i]);
+      g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]);
+    }
+    c[n-1]=0;
+    b[n-1]=x[n-2]-x[n-3];
+    a[n-1]=x[n-1]-x[n-3];
+    g[n-1]=((x[n-1]-x[n-2])^2*(y[n-2]-y[n-3])/b[n-1]+
+            b[n-1]*(2*b[n-1]+3(x[n-1]-x[n-2]))*
+            (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))/a[n-1];
+    d=tridiagonal(a,b,c,g);
+  } else if(n == 2) {
+    real val=(y[1]-y[0])/(x[1]-x[0]);
+    d=new real[] {val,val};
+  } else if(n == 3) {
+    real a=(y[1]-y[0])/(x[1]-x[0]);
+    real b=(y[2]-y[1])/(x[2]-x[1]);
+    real c=(b-a)/(x[2]-x[0]);
+    d=new real[] {a+c*(x[0]-x[1]),a+c*(x[1]-x[0]),a+c*(2*x[2]-x[0]-x[1])};
+  } else abort(morepoints);
+  return d;
+}
+
+// Standard cubic spline interpolation with periodic condition
+// s'(a)=s'(b), s''(a)=s''(b), assuming that f(a)=f(b)
+// if n=2, linear interpolation is returned
+real[] periodic(real[] x, real[] y)
+{
+  int n=x.length;
+  checklengths(n,y.length);
+  if(abs(y[n-1]-y[0]) > sqrtEpsilon*norm(y))
+    abort("function values are not periodic");
+  real[] d;
+  if(n > 2) {
+    real[] a=new real[n-1];
+    real[] b=new real[n-1];
+    real[] c=new real[n-1];
+    real[] g=new real[n-1];
+    c[0]=x[n-1]-x[n-2];
+    a[0]=x[1]-x[0];
+    b[0]=2*(a[0]+c[0]);
+    g[0]=3*c[0]*(y[1]-y[0])/a[0]+3*a[0]*(y[n-1]-y[n-2])/c[0];
+    for(int i=1; i < n-1; ++i) {
+      a[i]=x[i+1]-x[i];
+      c[i]=x[i]-x[i-1];
+      b[i]=2*(a[i]+c[i]);
+      g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]);
+    }
+    d=tridiagonal(a,b,c,g);
+    d.push(d[0]);
+  } else if(n == 2) {
+    d=new real[] {0,0};
+  } else abort(morepoints);
+  return d;
+}
+
+// Standard cubic spline interpolation with the natural condition
+// s''(a)=s''(b)=0.
+// if n=2, linear interpolation is returned
+// Don't use the natural type unless the underlying function 
+// has zero second end points derivatives.
+real[] natural(real[] x, real[] y)
+{
+  int n=x.length;
+  checklengths(n,y.length);
+  real[] d;
+  if(n > 2) {
+    real[] a=new real[n];
+    real[] b=new real[n];
+    real[] c=new real[n];
+    real[] g=new real[n];
+    b[0]=2*(x[1]-x[0]);
+    c[0]=x[1]-x[0];
+    a[0]=0;
+    g[0]=3*(y[1]-y[0]);
+    for(int i=1; i < n-1; ++i) {
+      a[i]=x[i+1]-x[i];
+      c[i]=x[i]-x[i-1];
+      b[i]=2*(a[i]+c[i]);
+      g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]);
+    }
+    c[n-1]=0;
+    a[n-1]=x[n-1]-x[n-2];
+    b[n-1]=2*a[n-1];
+    g[n-1]=3*(y[n-1]-y[n-2]);
+    d=tridiagonal(a,b,c,g);
+  } else if(n == 2) {
+    real val=(y[1]-y[0])/(x[1]-x[0]);
+    d=new real[] {val,val};
+  } else abort(morepoints);
+  return d;
+}
+
+// Standard cubic spline interpolation with clamped conditions f'(a), f'(b)
+splinetype clamped(real slopea, real slopeb)
+{
+  return new real[] (real[] x, real[] y) {
+    int n=x.length;
+    checklengths(n,y.length);
+    real[] d;
+    if(n > 2) {
+      real[] a=new real[n];
+      real[] b=new real[n];
+      real[] c=new real[n];
+      real[] g=new real[n];
+      b[0]=x[1]-x[0];
+      g[0]=b[0]*slopea;
+      c[0]=0;
+      a[0]=0;
+      for(int i=1; i < n-1; ++i) {
+        a[i]=x[i+1]-x[i];
+        c[i]=x[i]-x[i-1];
+        b[i]=2*(a[i]+c[i]);
+        g[i]=3*(c[i]*(y[i+1]-y[i])/a[i]+a[i]*(y[i]-y[i-1])/c[i]);
+      }
+      c[n-1]=0;
+      a[n-1]=0;
+      b[n-1]=x[n-1]-x[n-2];
+      g[n-1]=b[n-1]*slopeb;
+      d=tridiagonal(a,b,c,g);
+    } else if(n == 2) {
+      d=new real[] {slopea,slopeb};
+    } else abort(morepoints);
+    return d;
+  };
+}
+
+// Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)
+// Modified MATLAB code
+// [1] Fritsch, F. N. and R. E. Carlson, 
+//      "Monotone Piecewise Cubic Interpolation," 
+//      SIAM J. Numerical Analysis, Vol. 17, 1980, pp.238-246.
+// [2] Kahaner, David, Cleve Moler, Stephen Nash, 
+//      Numerical Methods and Software, Prentice Hall, 1988.
+real[] monotonic(real[] x, real[] y) 
+{
+  int n=x.length; 
+  checklengths(n,y.length);
+  real[] d=new real[n]; 
+  if(n > 2) {
+    real[] h=new real[n-1];
+    real[] del=new real[n-1];
+    for(int i=0; i < n-1; ++i) {
+      h[i]=x[i+1]-x[i]; 
+      del[i]=(y[i+1]-y[i])/h[i]; 
+    } 
+    int j=0; 
+    int k[]=new int[];
+    for(int i=0; i < n-2; ++i)
+      if((sgn(del[i])*sgn(del[i+1])) > 0) {k[j]=i; j=j+1;}
+
+    real[] hs=new real[j];
+    for(int i=0; i < j; ++i) hs[i]=h[k[i]]+h[k[i]+1];
+    real w1[]=new real[j];
+    real w2[]=new real[j];
+    real dmax[]=new real[j];
+    real dmin[]=new real[j];
+    for(int i=0; i < j; ++i) {
+      w1[i]=(h[k[i]]+hs[i])/(3*hs[i]);
+      w2[i]=(h[k[i]+1]+hs[i])/(3*hs[i]);
+      dmax[i]=max(abs(del[k[i]]),abs(del[k[i]+1]));
+      dmin[i]=min(abs(del[k[i]]),abs(del[k[i]+1]));
+    } 
+    for(int i=0; i < n; ++i) d[i]=0;
+    for(int i=0; i < j; ++i)
+      d[k[i]+1]=dmin[i]/(w1[i]*(del[k[i]]/dmax[i])+w2[i]*(del[k[i]+1]/dmax[i])); 
+    d[0]=((2*h[0]+h[1])*del[0]-h[0]*del[1])/(h[0]+h[1]);
+    if(sgn(d[0]) != sgn(del[0])) {d[0]=0;}
+    else if((sgn(del[0]) != sgn(del[1])) && (abs(d[0]) > abs(3*del[0])))
+      d[0]=3*del[0];
+
+    d[n-1]=((2*h[n-2]+h[n-3])*del[n-2]-h[n-2]*del[n-2])/(h[n-2]+h[n-3]);
+    if(sgn(d[n-1]) != sgn(del[n-2])) {d[n-1]=0;}
+    else if((sgn(del[n-2]) != sgn(del[n-3])) &&
+            (abs(d[n-1]) > abs(3*del[n-2])))
+      d[n-1]=3*del[n-2];
+  } else if(n == 2) {
+    d[0]=d[1]=(y[1]-y[0])/(x[1]-x[0]);
+  } else abort(morepoints);
+  return d;
+} 
+
+// Return standard cubic spline interpolation as a guide
+guide hermite(real[] x, real[] y, splinetype splinetype=null)
+{
+  int n=x.length;
+  if(n == 0) return nullpath;
+
+  guide g=(x[0],y[0]);
+  if(n == 1) return g;
+  if(n == 2) return g--(x[1],y[1]);
+
+  if(splinetype == null)
+    splinetype=(x[0] == x[x.length-1] && y[0] == y[y.length-1]) ?
+      periodic : notaknot;
+
+  real[] dy=splinetype(x,y);
+  for(int i=1; i < n; ++i) {
+    pair z=(x[i],y[i]);
+    real dx=x[i]-x[i-1];
+    g=g..controls((x[i-1],y[i-1])+dx*(1,dy[i-1])/3) and (z-dx*(1,dy[i])/3)..z;
+  }
+  return g;
+}
diff --git a/Master/texmf/asymptote/interpolate.asy b/Master/texmf/asymptote/interpolate.asy
index a6b820bc7d5..a600c94c125 100644
--- a/Master/texmf/asymptote/interpolate.asy
+++ b/Master/texmf/asymptote/interpolate.asy
@@ -29,7 +29,7 @@
 // Here s is a real function that is constant on (-infinity,a] and [b,infinity).
 
 private import math;
-import splinetype;
+import graph_splinetype;
 
 typedef real fhorner(real);
 
diff --git a/Master/texmf/asymptote/labelpath.asy b/Master/texmf/asymptote/labelpath.asy
index fe597a1608d..39c1908efd3 100644
--- a/Master/texmf/asymptote/labelpath.asy
+++ b/Master/texmf/asymptote/labelpath.asy
@@ -12,7 +12,7 @@ void labelpath(frame f, Label L, path g, string justify=Centered,
     _labelpath(f,L.s,L.size,g,justify,(L.T.x,L.T.y+0.5linewidth(p)),p);
     return;
   }
-  write("warning: labelpath requires -tex latex");
+  warning("labelpathlatex","labelpath requires -tex latex");
 }
 
 void labelpath(picture pic=currentpicture, Label L, path g,
diff --git a/Master/texmf/asymptote/labelpath3.asy b/Master/texmf/asymptote/labelpath3.asy
index bf2e6d4c084..8483a76a047 100644
--- a/Master/texmf/asymptote/labelpath3.asy
+++ b/Master/texmf/asymptote/labelpath3.asy
@@ -41,11 +41,11 @@ triple[] nextframe(path3 p, real reltimestart, triple[] start, real
   return bf[subdiv];
 }
   
-surface labelpath(string txt, path3 p, real angle=90, triple optional=O)
+surface labelpath(string s, path3 p, real angle=90, triple optional=O)
 {
   real Cos=Cos(angle);
   real Sin=Sin(angle);
-  path[] text=texpath(Label(txt,(0,0),Align,basealign));
+  path[] text=texpath(Label(s,(0,0),Align,basealign));
   text=scale(1/(max(text).x-min(text).x))*text;
   path[][] decompose=containmentTree(text);
 	
diff --git a/Master/texmf/asymptote/obj.asy b/Master/texmf/asymptote/obj.asy
index 02f839b594d..158d417c136 100644
--- a/Master/texmf/asymptote/obj.asy
+++ b/Master/texmf/asymptote/obj.asy
@@ -21,7 +21,7 @@ struct obj {
   pen[] meshpen;
 
   path3[][] read(string datafile, bool verbose=false) {
-    file in=word(line(input(datafile)));
+    file in=input(datafile).word().line();
     triple[] vert;
     path3[][] g;
     g[0]=new path3[] ;
@@ -92,8 +92,8 @@ struct obj {
                      material surfacepen, pen meshpen=nullpen) {
     material[] surfacepen={surfacepen};
     pen[] meshpen={meshpen};
-    surfacepen.cyclic(true);
-    meshpen.cyclic(true);
+    surfacepen.cyclic=true;
+    meshpen.cyclic=true;
     operator init(read(datafile,verbose),surfacepen,meshpen);
   }
 }
diff --git a/Master/texmf/asymptote/plain.asy b/Master/texmf/asymptote/plain.asy
index 8424c767ffa..f9e18b883ea 100644
--- a/Master/texmf/asymptote/plain.asy
+++ b/Master/texmf/asymptote/plain.asy
@@ -20,8 +20,8 @@ include plain_constants;
 
 access version;             
 if(version.VERSION != VERSION) {
-  write(stdout,"warning: using possibly incompatible version "+
-        version.VERSION+" of plain.asy"+'\n');
+  warning("version","using possibly incompatible version "+
+          version.VERSION+" of plain.asy"+'\n');
 }
    
 include plain_pens;
@@ -68,7 +68,7 @@ saveFunction[] saveFunctions={};
 // When save is called, this will be redefined to do the corresponding restore.
 void restore()
 {
-  write("warning: restore called with no matching save");
+  warning("nomatchingsave","restore called with no matching save");
 }
 
 void addSaveFunction(saveFunction s)
@@ -113,7 +113,8 @@ restoreThunk save()
 
 void restoredefaults()
 {
-  write("warning: restoredefaults called with no matching savedefaults");
+  warning("nomatchingsavedefaults",
+          "restoredefaults called with no matching savedefaults");
 }
 
 restoreThunk buildRestoreDefaults()
diff --git a/Master/texmf/asymptote/plain_Label.asy b/Master/texmf/asymptote/plain_Label.asy
index 3c11a2cedbd..b75d0de3967 100644
--- a/Master/texmf/asymptote/plain_Label.asy
+++ b/Master/texmf/asymptote/plain_Label.asy
@@ -335,7 +335,13 @@ struct Label {
     pic.add(new void (frame f, transform t) {
         out(f,t,point(g,position),
             alignrelative ? -Align*dir(t*g,position)*I : Align);
-      },true);
+      },!alignrelative);
+
+    frame f;
+    pair align=alignrelative ? -Align*dir(g,position)*I : Align;
+    label(f,(0,0),align);
+    pair position=point(g,position);
+    pic.addBox(position,position,min(f),max(f));
   }
   
   void write(file file=stdout, suffix suffix=endl) {
diff --git a/Master/texmf/asymptote/plain_arrows.asy b/Master/texmf/asymptote/plain_arrows.asy
index 0ef8c2b2d5c..18618c2b5b8 100644
--- a/Master/texmf/asymptote/plain_arrows.asy
+++ b/Master/texmf/asymptote/plain_arrows.asy
@@ -543,27 +543,36 @@ void arrow(picture pic=currentpicture, Label L="", pair b, pair dir,
   draw(b,pic,L,a--(0,0),align,p,arrow,margin);
 }
 
-// Force an array of 2D pictures to be as least as large as picture all.
-void rescale2(picture[] pictures, picture all)
-{
-  if(!all.empty2()) {
-    transform t=inverse(all.calculateTransform()*pictures[0].T);
-    pair m=t*min(all);
-    pair M=t*max(all);
-    for(int i=0; i < pictures.length; ++i) {
-      draw(pictures[i],m,nullpen);
-      draw(pictures[i],M,nullpen);
-    }
-  }
-}
-
-// Force an array of pictures to have a uniform scaling.
-void rescale(picture[] pictures)
-{
-  if(pictures.length == 0) return;
+// Fit an array of pictures simultaneously using the sizing of picture all.
+frame[] fit2(picture[] pictures, picture all)
+{
+ frame[] out;
+ if(!all.empty2()) {
+   transform t=all.calculateTransform();
+   pair m=all.min(t);
+   pair M=all.max(t);
+   for(picture pic : pictures) {
+     frame f=pic.fit(t);
+     draw(f,m,nullpen);
+     draw(f,M,nullpen);
+     out.push(f);
+   }
+ }
+  return out;
+}
+
+// Fit an array of pictures simultaneously using the size of the first picture.
+frame[] fit(string prefix="", picture[] pictures, string format="",
+            bool view=true, string options="", string script="",
+            projection P=currentprojection)
+{
+  if(pictures.length == 0)
+    return new frame[];
+ 
   picture all;
   size(all,pictures[0]);
   for(picture pic : pictures)
     add(all,pic);
-  rescale2(pictures,all);
+
+  return fit2(pictures,all);
 }
diff --git a/Master/texmf/asymptote/plain_boxes.asy b/Master/texmf/asymptote/plain_boxes.asy
index 7fb82c37c1f..ff82a0b0b32 100644
--- a/Master/texmf/asymptote/plain_boxes.asy
+++ b/Master/texmf/asymptote/plain_boxes.asy
@@ -126,11 +126,11 @@ pair point(object F, pair dir, transform t=identity())
   return z+point(F.g,T[0]);
 }
 
-frame bbox(string prefix="", picture pic=currentpicture,
+frame bbox(picture pic=currentpicture,
            real xmargin=0, real ymargin=xmargin,
            pen p=currentpen, filltype filltype=NoFill)
 {
-  frame f=pic.fit(prefix,max(pic.xsize-2*xmargin,0),max(pic.ysize-2*ymargin,0));
+  frame f=pic.fit(max(pic.xsize-2*xmargin,0),max(pic.ysize-2*ymargin,0));
   box(f,xmargin,ymargin,p,filltype,above=false);
   return f;
 }
diff --git a/Master/texmf/asymptote/plain_constants.asy b/Master/texmf/asymptote/plain_constants.asy
index c1d8c086f71..32138756891 100644
--- a/Master/texmf/asymptote/plain_constants.asy
+++ b/Master/texmf/asymptote/plain_constants.asy
@@ -67,6 +67,13 @@ void tab(file file) {write(file,'\t');}
 void comma(file file) {write(file,',');}
 typedef void suffix(file);
 
+// Used by interactive write to warn that the outputted type is the resolution
+// of an overloaded name.
+void overloadedMessage(file file) {
+  write(file,'    <overloaded>');
+  endl(file);
+}
+
 void write(suffix suffix=endl) {suffix(stdout);}
 void write(file file, suffix suffix=none) {suffix(file);}
 
diff --git a/Master/texmf/asymptote/plain_paths.asy b/Master/texmf/asymptote/plain_paths.asy
index d457fa5bc28..bf3f2790c1d 100644
--- a/Master/texmf/asymptote/plain_paths.asy
+++ b/Master/texmf/asymptote/plain_paths.asy
@@ -178,7 +178,7 @@ slice cut(path p, path knife, int n)
   slice s;
   real[][] T=intersections(p,knife);
   if(T.length == 0) {s.before=p; s.after=nullpath; return s;}
-  T.cyclic(true);
+  T.cyclic=true;
   real t=T[n][0];
   s.before=subpath(p,0,t);
   s.after=subpath(p,t,length(p));
diff --git a/Master/texmf/asymptote/plain_pens.asy b/Master/texmf/asymptote/plain_pens.asy
index fd2dfadde68..a8a0d0c7461 100644
--- a/Master/texmf/asymptote/plain_pens.asy
+++ b/Master/texmf/asymptote/plain_pens.asy
@@ -139,11 +139,11 @@ pen[] colorPen={red,blue,green,magenta,cyan,orange,purple,brown,
                 deepblue,deepgreen,chartreuse,fuchsia,lightred,
                 lightblue,black,pink,yellow,gray};
 
-colorPen.cyclic(true);
+colorPen.cyclic=true;
 
 pen[] monoPen={solid,dashed,dotted,longdashed,dashdotted,
                longdashdotted};
-monoPen.cyclic(true);
+monoPen.cyclic=true;
 
 pen Pen(int n) 
 {
diff --git a/Master/texmf/asymptote/plain_picture.asy b/Master/texmf/asymptote/plain_picture.asy
index b3836d7ef51..432a4a6ae6d 100644
--- a/Master/texmf/asymptote/plain_picture.asy
+++ b/Master/texmf/asymptote/plain_picture.asy
@@ -401,10 +401,18 @@ struct transformation {
   transform3 compute() {
     return projection*modelview;
   }
+  transformation copy() {
+    transformation T=new transformation;
+    T.modelview=copy(modelview);
+    T.projection=copy(projection);
+    T.infinity=infinity;
+    T.oblique=oblique;
+    return T;
+  }
 }
 
 struct projection {
-  transform3 t; // projection*modelview (cached)
+  transform3 t;         // projection*modelview (cached)
   bool infinity;
   bool oblique;
   bool absolute=false;
@@ -420,21 +428,16 @@ struct projection {
   bool autoadjust=true; // Adjust camera to lie outside bounding volume?
   bool center=false;    // Center target within bounding volume?
   int ninterpolate;     // Used for projecting nurbs to 2D Bezier curves.
+  transformation T;
 
   void calculate() {
-    transformation T=projector(camera,up,target);
+    T=projector(camera,up,target);
     t=T.compute();
     infinity=T.infinity;
     oblique=T.oblique;
     ninterpolate=infinity ? 1 : 16;
   }
 
-  transformation transformation() {
-    return projector(camera,up,target);
-  }
-
-  transform3 modelview() {return transformation().modelview;}
-
   triple vector() {
     return camera-target;
   }
@@ -473,6 +476,7 @@ struct projection {
     P.center=center;
     P.projector=projector;
     P.ninterpolate=ninterpolate;
+    P.T=T.copy();
     return P;
   }
 
@@ -500,6 +504,71 @@ struct projection {
 
 projection currentprojection;
 
+struct light {
+  real[][] diffuse;
+  real[][] ambient;
+  real[][] specular;
+  pen background=nullpen; // Background color of the 3D canvas.
+  real specularfactor;
+  bool viewport; // Are the lights specified (and fixed) in the viewport frame?
+  triple[] position; // Only directional lights are currently implemented.
+
+  transform3 T=identity(4); // Transform to apply to normal vectors.
+
+  bool on() {return position.length > 0;}
+  
+  void operator init(pen[] diffuse,
+                     pen[] ambient=array(diffuse.length,black),
+                     pen[] specular=diffuse, pen background=nullpen,
+                     real specularfactor=1,
+                     bool viewport=false, triple[] position) {
+    int n=diffuse.length;
+    assert(ambient.length == n && specular.length == n && position.length == n);
+    
+    this.diffuse=new real[n][];
+    this.ambient=new real[n][];
+    this.specular=new real[n][];
+    this.background=background;
+    this.position=new triple[n];
+    for(int i=0; i < position.length; ++i) {
+      this.diffuse[i]=rgba(diffuse[i]);
+      this.ambient[i]=rgba(ambient[i]);
+      this.specular[i]=rgba(specular[i]);
+      this.position[i]=unit(position[i]);
+    }
+    this.specularfactor=specularfactor;
+    this.viewport=viewport;
+  }
+
+  void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
+                     pen background=nullpen, real specularfactor=1,
+                     bool viewport=false...triple[] position) {
+    int n=position.length;
+    operator init(array(n,diffuse),array(n,ambient),array(n,specular),
+                  background,specularfactor,viewport,position);
+  }
+
+  void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
+                     pen background=nullpen, bool viewport=false,
+                     real x, real y, real z) {
+    operator init(diffuse,ambient,specular,background,viewport,(x,y,z));
+  }
+
+  void operator init(explicit light light) {
+    diffuse=copy(light.diffuse);
+    ambient=copy(light.ambient);
+    specular=copy(light.specular);
+    background=light.background;
+    specularfactor=light.specularfactor;
+    viewport=light.viewport;
+    position=copy(light.position);
+  }
+
+  real[] background() {return rgba(background == nullpen ? white : background);}
+}
+
+light currentlight;
+
 triple min3(pen p)
 {
   return linewidth(p)*(-0.5,-0.5,-0.5);
@@ -990,7 +1059,7 @@ struct picture {
     if(status == simplex.problem.OPTIMAL) {
       return scaling.build(p.a(),p.b()).a;
     } else if(status == simplex.problem.UNBOUNDED) {
-      if(warn) write("warning: "+dir+" scaling in picture unbounded");
+      if(warn) warning("unbounded",dir+" scaling in picture unbounded");
       return 0;
     } else {
       if(!warn) return 1;
@@ -1001,8 +1070,8 @@ struct picture {
           abort("unbounded picture");
       }
       if(userzero) return 1;
-      write("warning: cannot fit picture to "+dir+"size "+(string) size
-            +"...enlarging...");
+      warning("cannotfit","cannot fit picture to "+dir+"size "+(string) size
+              +"...enlarging...");
       return calculateScaling(dir,coords,sqrt(2)*size,warn);
     }
   }
@@ -1184,11 +1253,11 @@ struct picture {
     pair d=size(f);
     static real epsilon=100*realEpsilon;
     if(d.x > xsize*(1+epsilon)) 
-      write("warning: frame exceeds xlimit: "+(string) d.x+" > "+
-            (string) xsize);
+      warning("xlimit","frame exceeds xlimit: "+(string) d.x+" > "+
+              (string) xsize);
     if(d.y > ysize*(1+epsilon))
-      write("warning: frame exceeds ylimit: "+(string) d.y+" > "+
-            (string) ysize);
+      warning("ylimit","frame exceeds ylimit: "+(string) d.y+" > "+
+              (string) ysize);
     return f;
   }
   
@@ -1290,14 +1359,15 @@ struct picture {
   }
 
   static frame fitter(string,picture,string,real,real,bool,bool,string,string,
-                      projection);
+                      light,projection);
   frame fit(string prefix="", string format="",
             real xsize=this.xsize, real ysize=this.ysize,
             bool keepAspect=this.keepAspect, bool view=false,
-            string options="", string script="",
+            string options="", string script="", light light=currentlight,
             projection P=currentprojection) {
     return fitter == null ? fit2(xsize,ysize,keepAspect) :
-      fitter(prefix,this,format,xsize,ysize,keepAspect,view,options,script,P);
+      fitter(prefix,this,format,xsize,ysize,keepAspect,view,options,script,
+             light,P);
   }
   
   // In case only an approximate picture size estimate is available, return the
diff --git a/Master/texmf/asymptote/plain_shipout.asy b/Master/texmf/asymptote/plain_shipout.asy
index 22aeef7121f..2e86822cad0 100644
--- a/Master/texmf/asymptote/plain_shipout.asy
+++ b/Master/texmf/asymptote/plain_shipout.asy
@@ -27,7 +27,7 @@ typedef frame orientation(frame);
 orientation orientation=Portrait;
 
 // Forward references to functions defined in module three.
-object embed3(string, frame, string, string, string, projection);
+object embed3(string, frame, string, string, string, light, projection);
 string Embed(string name, string options="", real width=0, real height=0);
 string Link(string label, string text, string options="");
 
@@ -59,10 +59,10 @@ include plain_xasy;
 void shipout(string prefix=defaultfilename, frame f,
              string format="", bool wait=false, bool view=true,
 	     string options="", string script="",
-	     projection P=currentprojection)
+	     light light=currentlight, projection P=currentprojection)
 {
   if(is3D(f)) {
-    f=enclose(prefix,embed3(prefix,f,format,options,script,P));
+    f=enclose(prefix,embed3(prefix,f,format,options,script,light,P));
     if(settings.render != 0 && !prc(format)) {
       shipped=true;
       return;
@@ -88,14 +88,14 @@ void shipout(string prefix=defaultfilename, picture pic=currentpicture,
 	     orientation orientation=orientation,
 	     string format="", bool wait=false, bool view=true,
 	     string options="", string script="",
-	     projection P=currentprojection)
+	     light light=currentlight, projection P=currentprojection)
 {
   if(!uptodate()) {
     bool inlinetex=settings.inlinetex;
     bool prc=prc(format);
     if(prc && !pic.empty3())
       settings.inlinetex=settings.inlineimage;
-    frame f=pic.fit(prefix,format,view=view,options,script,P);    
+    frame f=pic.fit(prefix,format,view=view,options,script,light,P);
     if(!pic.empty2() || settings.render == 0 || prc)
       shipout(prefix,orientation(f),format,wait,view);
     settings.inlinetex=inlinetex;
diff --git a/Master/texmf/asymptote/slide.asy b/Master/texmf/asymptote/slide.asy
index b154c1d4746..b23a352fc5d 100644
--- a/Master/texmf/asymptote/slide.asy
+++ b/Master/texmf/asymptote/slide.asy
@@ -217,11 +217,10 @@ bool checkposition()
   return true;
 }
 
-void step()
-{
+void erasestep(int erasenode) {
   if(!stepping || !allowstepping) return;
   if(!checkposition()) return;
-  lastnode.push(currentpicture.nodes.length-1);
+  lastnode.push(erasenode);
   nextpage(steppagenumberpen);
   for(int i=0; i < firstnode.length; ++i) {
     for(int j=firstnode[i]; j <= lastnode[i]; ++j) {
@@ -233,6 +232,12 @@ void step()
   tex(bulletcolor(newbulletcolor));
 }
 
+void step()
+{
+  // Step without erasing anything.
+  erasestep(currentpicture.nodes.length-1);
+}
+
 void incrementposition(pair z)
 {
   currentposition += z;
@@ -279,11 +284,13 @@ void remark(bool center=false, string s, pair align=0, pen p=itempen,
   pair M=tinv*min(f);
   
   if(abs(offset.x+M.x) > 1)
-    write("warning: slide too wide on page "+(string) page+':\n'+(string) s);
+    warning("slidetoowide","slide too wide on page "+(string) page+':\n'+
+            (string) s);
 
   if(abs(offset.y+M.y) > 1) {
     void toohigh() {
-      write("warning: slide too high on page "+(string) page+':\n'+(string) s);
+      warning("slidetoohigh","slide too high on page "+(string) page+':\n'+
+              (string) s);
     }
     if(M.y-m.y < 2) {
       newslide(); offset=(offset.x,currentposition.y);
@@ -326,7 +333,7 @@ void skip(real n=1)
 }
 
 void display(frame[] f, real margin=0, pair align=S, pen p=itempen,
-             pen figuremattpen=figuremattpen)
+             pen figuremattpen=figuremattpen, bool final=true)
 {
   if(f.length == 0) return;
   real[] width=new real[f.length];
@@ -336,7 +343,7 @@ void display(frame[] f, real margin=0, pair align=S, pen p=itempen,
     sum += width[i];
   }
   if(sum > pagewidth)
-    write("warning: slide too wide on page "+(string) page);
+    warning("toowide","slide too wide on page "+(string) page);
   else margin=(pagewidth-sum)/(f.length+1);
   real pos;
   frame F;
@@ -347,19 +354,22 @@ void display(frame[] f, real margin=0, pair align=S, pen p=itempen,
     pos += w;
   }
   add(F,(0,currentposition.y),align);
-  real a=0.5(unit(align).y-1);
-  incrementposition((0,(tinv*(a*(max(F)-min(F))-itemskip*I*lineskip(p)*pt)).y));
+  if (final) {
+    real a=0.5(unit(align).y-1);
+    incrementposition(
+        (0, (tinv*(a*(max(F)-min(F))-itemskip*I*lineskip(p)*pt)).y));
+  }
 }
 
 void display(frame f, real margin=0, pair align=S, pen p=itempen,
-             pen figuremattpen=figuremattpen)
+             pen figuremattpen=figuremattpen, bool final=true)
 {
-  display(new frame[] {f},margin,align,p,figuremattpen);
+  display(new frame[] {f},margin,align,p,figuremattpen, final);
 }
 
 void display(string[] s, real margin=0, string[] captions=new string[],
              string caption="", pair align=S, pen p=itempen,
-             pen figuremattpen=figuremattpen)
+             pen figuremattpen=figuremattpen, bool final=true)
 {
   frame[] f=new frame[s.length];
   frame F;
@@ -374,32 +384,59 @@ void display(string[] s, real margin=0, string[] captions=new string[],
     if(captions[i] != "")
       label(f[i],captions[i],point(f[i],S).x+I*y,S);
   }
-  display(f,margin,align,p,figuremattpen);
+  display(f,margin,align,p,figuremattpen, final);
   if(caption != "") center(caption,p);
 }
 
 void display(string s, string caption="", pair align=S, pen p=itempen,
-             pen figuremattpen=figuremattpen)
+             pen figuremattpen=figuremattpen, bool final=true)
 {
-  display(new string[] {s},caption,align,p,figuremattpen);
+  display(new string[] {s},caption,align,p,figuremattpen, final);
 }
 
 void figure(string[] s, string options="", real margin=0, 
             string[] captions=new string[], string caption="",
-            pair align=S, pen p=itempen, pen figuremattpen=figuremattpen)
+            pair align=S, pen p=itempen, pen figuremattpen=figuremattpen,
+            bool final=true)
 {
   string[] S;
   for(int i=0; i < s.length; ++i) {
     S[i]=graphic(s[i],options);
   }
 
-  display(S,margin,captions,caption,align,itempen,figuremattpen);
+  display(S,margin,captions,caption,align,itempen,figuremattpen, final);
 }
 
 void figure(string s, string options="", string caption="", pair align=S,
-            pen p=itempen, pen figuremattpen=figuremattpen)
+            pen p=itempen, pen figuremattpen=figuremattpen, bool final=true)
+{
+  figure(new string[] {s},options,caption,align,p,figuremattpen, final);
+}
+
+void multifigure(string[] slist, string options="", string caption="",
+                 pair align=S, pen p=itempen, pen figuremattpen=figuremattpen)
+{
+  if (stepping) {
+    int lastnode = currentpicture.nodes.length-1;
+    for (int i=0; i<slist.length-1; ++i) {
+      figure(slist[i], options, caption, align, p, figuremattpen, final=false);
+      erasestep(lastnode);
+    }
+  }
+  figure(slist[slist.length-1], options, caption, align, p, figuremattpen,
+         final=true);
+
+  if(!firststep) step();
+  firststep=false;
+}
+
+void indexedfigure(string prefix, int n, string options="", string caption="",
+                 pair align=S, pen p=itempen, pen figuremattpen=figuremattpen)
 {
-  figure(new string[] {s},options,caption,align,p,figuremattpen);
+  string[] s;
+  for (int i=0; i<n; ++i)
+    s.push(prefix+"+"+string(i));
+  multifigure(s, options, caption, align, p, figuremattpen);
 }
 
 string[] codefile;
diff --git a/Master/texmf/asymptote/solids.asy b/Master/texmf/asymptote/solids.asy
index 0b2ffa2a508..db98af3088c 100644
--- a/Master/texmf/asymptote/solids.asy
+++ b/Master/texmf/asymptote/solids.asy
@@ -212,7 +212,8 @@ struct revolution {
   // must be recomputed if camera is adjusted
   path3[] silhouette(int m=64, projection P=currentprojection) {
     if(is3D())
-      write("warning: silhouette routine is intended only for 2d projections");
+      warning("2Dsilhouette",
+              "silhouette routine is intended only for 2d projections");
     path3 G,H;
     int N=size(g);
     int M=(m == 0) ? N : m;
diff --git a/Master/texmf/asymptote/syzygy.asy b/Master/texmf/asymptote/syzygy.asy
index dce87805c9d..93889b252e5 100644
--- a/Master/texmf/asymptote/syzygy.asy
+++ b/Master/texmf/asymptote/syzygy.asy
@@ -669,16 +669,6 @@ Braid apply(Relation r, Braid b, int step, int place) {
 }
 
 // Tableau {{{1
-frame[] fit(picture[] pics) {
-  frame[] f;
-  for (int i=0; i<pics.length; ++i) {
-    frame ff=pics[i].fit();
-    //label(ff, (string)i, (10,10));
-    //f.push(pics[i].fit());
-    f.push(ff);
-  }
-  return f;
-}
 
 // Draw a number of frames in a nice circular arrangement.
 picture tableau(frame[] cards, bool number=false) {
diff --git a/Master/texmf/asymptote/three.asy b/Master/texmf/asymptote/three.asy
index 39075193073..1f9392ee094 100644
--- a/Master/texmf/asymptote/three.asy
+++ b/Master/texmf/asymptote/three.asy
@@ -13,13 +13,11 @@ if(prc0()) {
 
 real defaultshininess=0.25;
 real defaultgranularity=0;
-real linegranularity=0.01;
-real tubegranularity=0.003;
+real linegranularity=0.005;
+int linesectors=8;        // Number of angular sectors.
 real dotgranularity=0.0001;
-real viewportfactor=1.002;   // Factor used to expand orthographic viewport.
-real angleprecision=1e-3;    // Precision for centering perspective projections.
-real anglefactor=max(1.01,1+angleprecision);
-// Factor used to expand perspective viewport.
+real angleprecision=1e-5; // Precision for centering perspective projections.
+real rendermargin=0.02;
 
 string defaultembed3Doptions;
 string defaultembed3Dscript;
@@ -137,7 +135,9 @@ triple project(triple u, triple v)
 triple perp(triple v)
 {
   triple u=cross(v,Y);
-  return (abs(u) > sqrtEpsilon) ? unit(u) : unit(cross(v,Z));
+  if(abs(u) > sqrtEpsilon) return unit(u);
+  u=cross(v,Z);
+  return (abs(u) > sqrtEpsilon) ? unit(u) : X;
 }
 
 // Return the transformation corresponding to moving the camera from the target
@@ -757,7 +757,7 @@ real[] theta(triple[] v, real[] alpha, real[] beta,
 {
   real[] a,b,c,f,l,psi;
   int n=alpha.length;
-  bool cyclic=v.cyclicflag;
+  bool cyclic=v.cyclic;
   for(int i=0; i < n; ++i)
     l[i]=1/length(v[i+1]-v[i]);
   int i0,in;
@@ -766,8 +766,8 @@ real[] theta(triple[] v, real[] alpha, real[] beta,
   for(int i=0; i < in; ++i)
     psi[i]=angle(v[i+1]-v[i],v[i+2]-v[i+1],reference);
   if(cyclic) {
-    l.cyclic(true);
-    psi.cyclic(true);
+    l.cyclic=true;
+    psi.cyclic=true;
   } else {
     psi[n-1]=0;
     if(dir0 == O) {
@@ -869,9 +869,9 @@ void aim(flatguide3 g, int N)
     beta[k]=g.Tension[K].in;
   }
   if(cyclic) {
-    v.cyclic(true);
-    alpha.cyclic(true);
-    beta.cyclic(true);
+    v.cyclic=true;
+    alpha.cyclic=true;
+    beta.cyclic=true;
   } else v[n]=g.nodes[(start+n) % N];
   int final=(end-1) % N;
 
@@ -883,8 +883,8 @@ void aim(flatguide3 g, int N)
   real[] theta=theta(v,alpha,beta,d0,d1,g.out[start].gamma,g.in[final].gamma,
                      reference);
 
-  v.cyclic(true);
-  theta.cyclic(true);
+  v.cyclic=true;
+  theta.cyclic=true;
     
   for(int k=1; k < (cyclic ? n+1 : n); ++k) {
     triple w=dir(theta[k],v[k]-v[k-1],v[k+1]-v[k],reference);
@@ -1029,6 +1029,17 @@ path3[] path3(explicit path[] g, triple plane(pair)=XYplane)
   return sequence(new path3(int i) {return path3(g[i],plane);},g.length);
 }
 
+path3 invert(path p, triple normal, triple point,
+             projection P=currentprojection)
+{
+  return path3(p,new triple(pair z) {return invert(z,normal,point,P);});
+}
+
+path3 invert(path p, projection P=currentprojection)
+{
+  return path3(p,new triple(pair z) {return invert(z,P);});
+}
+
 // Construct a path from a path3 by applying P to each control point.
 path path(path3 p, pair P(triple)=xypart)
 {
@@ -1801,6 +1812,16 @@ triple[] operator * (transform3 t, triple[] v)
   return sequence(new triple(int i) {return t*v[i];},v.length);
 }
 
+triple[][] operator * (transform3 t, triple[][] v) 
+{
+  triple[][] V=new triple[v.length][];
+  for(int i=0; i < v.length; ++i) {
+    triple[] vi=v[i];
+    V[i]=sequence(new triple(int j) {return t*vi[j];},vi.length);
+  }
+  return V;
+}
+
 triple min(explicit path3[] p)
 {
   checkEmpty(p.length);
@@ -1978,6 +1999,38 @@ void addPath(picture pic, path3 g, pen p)
 include three_surface;
 include three_margins;
 
+pair min(path3 p, projection P) 
+{
+  path3 q=P.T.modelview*p;
+  if(P.infinity)
+    return xypart(min(q));
+  return maxratio(q)/P.T.projection[3][2];
+}
+
+pair max(path3 p, projection P) 
+{
+  path3 q=P.T.modelview*p;
+  if(P.infinity)
+    return xypart(max(q));
+  return minratio(q)/P.T.projection[3][2];
+}
+
+pair min(frame f, projection P) 
+{
+  frame g=P.T.modelview*f;
+  if(P.infinity)
+    return xypart(min3(g));
+  return maxratio(g)/P.T.projection[3][2];
+}
+
+pair max(frame f, projection P) 
+{
+  frame g=P.T.modelview*f;
+  if(P.infinity)
+    return xypart(max3(g));
+  return minratio(g)/P.T.projection[3][2];
+}
+
 void draw(picture pic=currentpicture, Label L="", path3 g,
           align align=NoAlign, material p=currentpen, margin3 margin=NoMargin3,
           light light=nolight)
@@ -1988,8 +2041,8 @@ void draw(picture pic=currentpicture, Label L="", path3 g,
       if(is3D()) {
         draw(f,G,p,light,null);
         if(pic != null && size(G) > 0) {
-          pic.addPoint(min(G,P.t));
-          pic.addPoint(max(G,P.t));
+          pic.addPoint(min(G,P));
+          pic.addPoint(max(G,P));
         }
       }
       if(pic != null)
@@ -2004,7 +2057,7 @@ void draw(picture pic=currentpicture, Label L="", path3 g,
   addPath(pic,g,q);
 }
 
-include three_arrows;
+include three_tube;
 
 draw=new void(frame f, path3 g, material p=currentpen,
               light light=nolight, projection P=currentprojection) {
@@ -2015,7 +2068,7 @@ draw=new void(frame f, path3 g, material p=currentpen,
       if(settings.thick) {
         real width=linewidth(q);
         if(width > 0) {
-          tube T=tube(g,width,p.granularity);
+          tube T=tube(g,width,linesectors);
           int L=length(g);
           if(L >= 0) {
             if(!cyclic(g)) {
@@ -2050,7 +2103,7 @@ draw=new void(frame f, path3 g, material p=currentpen,
     else {
       real[] dash=(real[]) split(type," ");
       if(sum(dash) > 0) {
-        dash.cyclic(true);
+        dash.cyclic=true;
         real offset=offset(q);
         real L=arclength(g);
         int i=0;
@@ -2081,6 +2134,8 @@ void draw(picture pic=currentpicture, explicit path3[] g,
   for(int i=0; i < g.length; ++i) draw(pic,g[i],p,margin,light);
 }
 
+include three_arrows;
+
 void draw(picture pic=currentpicture, Label L="", path3 g, 
           align align=NoAlign, material p=currentpen, arrowbar3 arrow,
           arrowbar3 bar=None, margin3 margin=NoMargin3, light light=nolight,
@@ -2321,7 +2376,7 @@ pair viewportmargin(pair lambda)
 
 string embed3D(string label="", string text=label, string prefix,
                frame f, string format="",
-               real width=0, real height=0, real angle=30,
+               real width=0, real height=0,
                string options="", string script="",
                light light=currentlight, projection P=currentprojection)
 {
@@ -2341,7 +2396,7 @@ string embed3D(string label="", string text=label, string prefix,
     pair margin=viewportmargin((lambda.x,lambda.y));
     viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/(cm*P.zoom);
   } else
-    if(!P.absolute) angle=2*aTan(Tan(0.5*angle));
+    if(!P.absolute) P.angle=2*aTan(Tan(0.5*P.angle));
 
   string name=prefix+".js";
   writeJavaScript(name,lightscript+projection(P.infinity,viewplanesize),script);
@@ -2375,7 +2430,7 @@ string embed3D(string label="", string text=label, string prefix,
     options3 += ",poster";
   options3 += ",text={"+text+"},label="+label+
     ",toolbar="+(settings.toolbar ? "true" : "false")+
-    ",3Daac="+Format(P.absolute ? P.angle : angle)+
+    ",3Daac="+Format(P.angle)+
     ",3Dc2c="+Format(u)+
     ",3Dcoo="+Format(target/cm)+
     ",3Droll="+Format(roll)+
@@ -2387,302 +2442,350 @@ string embed3D(string label="", string text=label, string prefix,
   return Embed(stripdirectory(prefix),options3,width,height);
 }
 
-object embed(string label="", string text=label, 
-             string prefix=defaultfilename, 
-             frame f, string format="",
-             real width=0, real height=0, real angle=30,
-             string options="", string script="", 
-             light light=currentlight, projection P=currentprojection)
-{
-  object F;
-
-  if(is3D(format))
-    F.L=embed3D(label,text,prefix,f,format,width,height,angle,options,script,
-                light,P);
-  else
-    F.f=f;
-  return F;
-}
-
-object embed(string label="", string text=label,
-             string prefix=defaultfilename,
-             picture pic, string format="",
-             real xsize=pic.xsize, real ysize=pic.ysize,
-             bool keepAspect=pic.keepAspect, bool view=true, string options="",
-             string script="", real angle=0,
-             light light=currentlight, projection P=currentprojection)
+struct scene
 {
-  object F;
-  real xsize3=pic.xsize3, ysize3=pic.ysize3, zsize3=pic.zsize3;
-  bool warn=true;
-  transform3 modelview;
-        
-  if(xsize3 == 0 && ysize3 == 0 && zsize3 == 0) {
-    xsize3=ysize3=zsize3=max(xsize,ysize);
-    warn=false;
+  frame f;
+  transform3 t;
+  projection P;
+  bool adjusted;
+  real width,height;
+  transform3 T=identity4;
+  picture pic2;
+  
+  void operator init(frame f, projection P=currentprojection) {
+    this.f=f;
+    this.t=identity4;
+    this.P=P;
   }
+  
+  void operator init(picture pic, real xsize=pic.xsize, real ysize=pic.ysize,
+                     bool keepAspect=pic.keepAspect, bool is3D=true,
+                     projection P=currentprojection) {
+    real xsize3=pic.xsize3, ysize3=pic.ysize3, zsize3=pic.zsize3;
+    bool warn=true;
+        
+    if(xsize3 == 0 && ysize3 == 0 && zsize3 == 0) {
+      xsize3=ysize3=zsize3=max(xsize,ysize);
+      warn=false;
+    }
 
-  projection P=P.copy();
+    if(P.absolute)
+      this.P=P.copy();
+    else if(P.showtarget)
+      draw(pic,P.target,nullpen);
 
-  if(!P.absolute && P.showtarget)
-    draw(pic,P.target,nullpen);
+    t=pic.scaling(xsize3,ysize3,zsize3,keepAspect,warn);
+    adjusted=false;
+    triple m=pic.min(t);
+    triple M=pic.max(t);
 
-  transform3 t=pic.scaling(xsize3,ysize3,zsize3,keepAspect,warn);
-  bool adjusted=false;
-  transform3 tinv=inverse(t);
-  triple m=pic.min(t);
-  triple M=pic.max(t);
+    if(!P.absolute) {
+      this.P=t*P;
+      if(this.P.center) {
+        bool recalculate=false;
+        triple target=0.5*(m+M);
+        this.P.target=target;
+        recalculate=true;
+        if(recalculate) this.P.calculate();
+      }
+      if(this.P.autoadjust || this.P.infinity) 
+        adjusted=adjusted | this.P.adjust(m,M);
+    }
 
-  if(!P.absolute) {
-    P=t*P;
-    if(P.center) {
-      bool recalculate=false;
-      triple target=0.5*(m+M);
-      P.target=target;
-      recalculate=true;
-      if(recalculate) P.calculate();
+    f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,this.P);
+
+    if(!pic.bounds3.exact) {
+      transform3 s=pic.scale3(f,xsize3,ysize3,zsize3,keepAspect);
+      t=s*t;
+      this.P=s*this.P;
+      f=pic.fit3(t,pic2,this.P);
     }
-    if(P.autoadjust || P.infinity) 
-      adjusted=adjusted | P.adjust(m,M);
-  }
 
-  picture pic2;
-  
-  frame f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,P);
+    bool scale=xsize != 0 || ysize != 0;
+
+    if(is3D || scale) {
+      pic2.bounds.exact=true;
+      transform s=pic2.scaling(xsize,ysize,keepAspect);
+
+      pair m2=pic2.min(s);
+      pair M2=pic2.max(s);
+      pair lambda=M2-m2;
+      pair viewportmargin=viewportmargin(lambda);
+      width=ceil(lambda.x+2*viewportmargin.x);
+      height=ceil(lambda.y+2*viewportmargin.y);
+
+      if(!this.P.absolute) {
+        if(scale && this.P.autoadjust) {
+          pair v=(s.xx,s.yy);
+          transform3 T=this.P.t;
+          pair x=project(X,T);
+          pair y=project(Y,T);
+          pair z=project(Z,T);
+          real f(pair a, pair b) {
+            return b == 0 ? (0.5*(a.x+a.y)) :
+              (b.x^2*a.x+b.y^2*a.y)/(b.x^2+b.y^2);
+          }
+          pic2.erase();
+          transform3 s=keepAspect ? scale3(min(f(v,x),f(v,y),f(v,z))) :
+            xscale3(f(v,x))*yscale3(f(v,y))*zscale3(f(v,z));
+          s=shift(this.P.target)*s*shift(-this.P.target);
+          t=s*t;
+          this.P=s*this.P;
+          f=pic.fit3(t,is3D ? null : pic2,this.P);
+        }
 
-  if(!pic.bounds3.exact) {
-    transform3 s=pic.scale3(f,xsize3,ysize3,zsize3,keepAspect);
-    t=s*t;
-    tinv=inverse(t);
-    P=s*P;
-    f=pic.fit3(t,pic2,P);
+        if(this.P.autoadjust || this.P.infinity)
+          adjusted=adjusted | this.P.adjust(min3(f),max3(f));
+      }
+    }
   }
 
-  bool is3D=is3D(format);
-  bool scale=xsize != 0 || ysize != 0;
-
-  if(is3D || scale) {
-    pic2.bounds.exact=true;
-    transform s=pic2.scaling(xsize,ysize,keepAspect);
-    pair m2=pic2.min(s);
-    pair M2=pic2.max(s);
-    pair lambda=M2-m2;
-    pair viewportmargin=viewportmargin(lambda);
-    real width=ceil(lambda.x+2*viewportmargin.x);
-    real height=ceil(lambda.y+2*viewportmargin.y);
-
-    projection Q;
-    if(!P.absolute) {
-      if(scale && P.autoadjust) {
-        pair v=(s.xx,s.yy);
-        transform3 T=P.t;
-        pair x=project(X,T);
-        pair y=project(Y,T);
-        pair z=project(Z,T);
-        real f(pair a, pair b) {
-          return b == 0 ? (0.5*(a.x+a.y)) : (b.x^2*a.x+b.y^2*a.y)/(b.x^2+b.y^2);
+  // Choose the angle to be just large enough to view the entire image.
+  real angle(projection P) {
+    T=identity4;
+    int maxiterations=100;
+    real h=-0.5*P.target.z;
+    pair r,R;
+    real diff=realMax;
+    pair s;
+    int i;
+    do {
+      r=minratio(f);
+      R=maxratio(f);
+      pair lasts=s;
+      if(P.autoadjust) {
+        s=r+R;
+        if(s != 0) {
+          transform3 t=shift(h*s.x,h*s.y,0);
+          f=t*f;
+          T=t*T;
+          adjusted=true;
         }
-        pic2.erase();
-        transform3 s=keepAspect ? scale3(min(f(v,x),f(v,y),f(v,z))) :
-          xscale3(f(v,x))*yscale3(f(v,y))*zscale3(f(v,z));
-        s=shift(P.target)*s*shift(-P.target);
-        t=s*t;
-        P=s*P;
-        f=pic.fit3(t,is3D ? null : pic2,P);
       }
+      diff=abs(s-lasts);
+      ++i;
+    } while (diff > angleprecision && i < maxiterations);
+    real aspect=width > 0 ? height/width : 1;
+    real rx=-r.x*aspect;
+    real Rx=R.x*aspect;
+    real ry=-r.y;
+    real Ry=R.y;
+    if(!P.autoadjust) {
+      if(rx > Rx) Rx=rx;
+      else rx=Rx;
+      if(ry > Ry) Ry=ry;
+      else ry=Ry;
+    }
+    return (1+angleprecision)*max(aTan(rx)+aTan(Rx),aTan(ry)+aTan(Ry));
+  }
+}
 
-      if(P.autoadjust || P.infinity)
-        adjusted=adjusted | P.adjust(min3(f),max3(f));
+object embed(string label="", string text=label, string prefix=defaultfilename,
+             scene S, string format="", bool view=true, string options="",
+             string script="", light light=currentlight)
+{
+  object F;
+  transform3 modelview;
+  projection P=S.P;
+  transform3 tinv=inverse(S.t);
 
-      triple target=P.target;
-      modelview=P.modelview();
-      f=modelview*f;
-      P=modelview*P;
-      Q=P.copy();
-      light=modelview*light;
-
-      if(P.infinity) {
-        triple m=min3(f);
-        triple M=max3(f);
-        triple lambda=M-m;
-        viewportmargin=viewportmargin((lambda.x,lambda.y));
-        width=lambda.x+2*viewportmargin.x;
-        height=lambda.y+2*viewportmargin.y;
-        f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*f;  // Eye will be at (0,0,0).
-      } else {
-        transform3 T=identity4;
-        // Choose the angle to be just large enough to view the entire image:
-        if(angle == 0) angle=P.angle;
-        int maxiterations=100;
-        if(is3D && angle == 0) {
-          real h=-0.5*P.target.z;
-          pair r,R;
-          real diff=realMax;
-          pair s;
-          int i;
-          do {
-            r=minratio(f);
-            R=maxratio(f);
-            pair lasts=s;
-            if(P.autoadjust) {
-              s=r+R;
-              if(s != 0) {
-                transform3 t=shift(h*s.x,h*s.y,0);
-                f=t*f;
-                T=t*T;
-                adjusted=true;
-              }
-            }
-            diff=abs(s-lasts);
-            ++i;
-          } while (diff > angleprecision && i < maxiterations);
-          real aspect=width > 0 ? height/width : 1;
-          real rx=-r.x*aspect;
-          real Rx=R.x*aspect;
-          real ry=-r.y;
-          real Ry=R.y;
-          if(!P.autoadjust) {
-            if(rx > Rx) Rx=rx;
-            else rx=Rx;
-            if(ry > Ry) Ry=ry;
-            else ry=Ry;
-          }
-          
-          angle=anglefactor*max(aTan(rx)+aTan(Rx),aTan(ry)+aTan(Ry));
-          if(viewportmargin.y != 0)
-            angle=2*aTan(Tan(0.5*angle)-viewportmargin.y/P.target.z);
-          
-          modelview=T*modelview;
-        }
-        if(settings.verbose > 0) {
-          transform3 inv=inverse(modelview);
-          if(adjusted) 
-            write("adjusting camera to ",tinv*inv*P.camera);
-          target=inv*P.target;
-        }
-        P=T*P;
+  projection Q;
+  modelview=P.T.modelview;
+  if(P.absolute) {
+    Q=modelview*P;
+  } else {
+    triple target=P.target;
+    S.f=modelview*S.f;
+    P=modelview*P;
+    Q=P.copy();
+    light=modelview*light;
+
+    pair viewportmargin=viewportmargin;
+    if(P.infinity) {
+      triple m=min3(S.f);
+      triple M=max3(S.f);
+      triple lambda=M-m;
+      viewportmargin=viewportmargin((lambda.x,lambda.y));
+      S.width=lambda.x+2*viewportmargin.x;
+      S.height=lambda.y+2*viewportmargin.y;
+      S.f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*S.f; // Eye will be at (0,0,0)
+    } else {
+      if(P.angle == 0) {
+        Q.angle=P.angle=S.angle(P);
+        modelview=S.T*modelview;
+        if(viewportmargin.y != 0)
+          P.angle=2*aTan(Tan(0.5*P.angle)-viewportmargin.y/P.target.z);
       }
       if(settings.verbose > 0) {
-        if(P.center || (!P.infinity && P.autoadjust))
-          write("adjusting target to ",tinv*target);
+        transform3 inv=inverse(modelview);
+        if(S.adjusted) 
+          write("adjusting camera to ",tinv*inv*P.camera);
+        target=inv*P.target;
       }
+      P=S.T*P;
     }
+    if(settings.verbose > 0) {
+      if(P.center || (!P.infinity && P.autoadjust))
+        write("adjusting target to ",tinv*target);
+    }
+  }
     
-    if(prefix == "") prefix=outprefix();
-    bool prc=prc(format);
-    bool preview=settings.render > 0;
-    if(prc) {
-      // The movie15.sty package cannot handle spaces or dots in filenames.
-      prefix=replace(prefix,new string[][]{{" ","_"},{".","_"}});
-      if(settings.embed || nativeformat() == "pdf")
-        prefix += "+"+(string) file3.length;
-    } else
-      preview=false;
-    if(preview || (!prc && settings.render != 0)) {
-      frame f=f;
-      triple m,M;
-      real zcenter;
-      if(P.absolute) {
-        modelview=P.modelview();
-        f=modelview*f;
-        Q=modelview*P;
-        m=min3(f);
-        M=max3(f);
-        real r=0.5*abs(M-m);
-        zcenter=0.5*(M.z+m.z);
-        M=(M.x,M.y,zcenter+r);
-        m=(m.x,m.y,zcenter-r);
-        angle=P.angle;
-      } else {
-        m=min3(f);
-        M=max3(f);
-        zcenter=P.target.z;
-        real d=P.distance(m,M);
-        M=(M.x,M.y,zcenter+d);
-        m=(m.x,m.y,zcenter-d);
-      }
-
-      if(P.infinity) {
-        triple margin=(viewportfactor-1.0)*(abs(M.x-m.x),abs(M.y-m.y),0)
-          +(viewportmargin.x,viewportmargin.y,0);
-        M += margin; 
-        m -= margin;
-      } else if(M.z >= 0) abort("camera too close");
-
-      shipout3(prefix,f,preview ? nativeformat() : format,
-               width,height,P.infinity ? 0 : 2aTan(Tan(0.5*angle)*P.zoom),
-               P.zoom,m,M,P.viewportshift,
-               tinv*inverse(modelview)*shift(0,0,zcenter),light.background(),
-               P.absolute ? (modelview*light).position : light.position,
-               light.diffuse,light.ambient,light.specular,
-               light.viewport,view && !preview);
-      if(!preview) return F;
+  if(prefix == "") prefix=outprefix();
+  bool prc=prc(format);
+  bool preview=settings.render > 0;
+  if(prc) {
+    // The movie15.sty package cannot handle spaces or dots in filenames.
+    prefix=replace(prefix,new string[][]{{" ","_"},{".","_"}});
+    if(settings.embed || nativeformat() == "pdf")
+      prefix += "+"+(string) file3.length;
+  } else
+    preview=false;
+  if(preview || (!prc && settings.render != 0)) {
+    frame f=S.f;
+    triple m,M;
+    real zcenter;
+    real r;
+    if(P.absolute) {
+      f=modelview*f;
+      m=min3(f);
+      M=max3(f);
+      r=0.5*abs(M-m);
+      zcenter=0.5*(M.z+m.z);
+    } else {
+      m=min3(f);
+      M=max3(f);
+      zcenter=P.target.z;
+      r=P.distance(m,M);
     }
+    M=(M.x,M.y,zcenter+r);
+    m=(m.x,m.y,zcenter-r);
+
+    if(P.infinity) {
+      triple margin=(viewportmargin.x,viewportmargin.y,0);
+      M += margin; 
+      m -= margin;
+    } else if(M.z >= 0) abort("camera too close");
+
+    shipout3(prefix,f,preview ? nativeformat() : format,
+             S.width-rendermargin,S.height-rendermargin,
+             P.infinity ? 0 : 2aTan(Tan(0.5*P.angle)*P.zoom),
+             P.zoom,m,M,P.viewportshift,
+             tinv*inverse(modelview)*shift(0,0,zcenter),light.background(),
+             P.absolute ? (modelview*light).position : light.position,
+             light.diffuse,light.ambient,light.specular,
+             light.viewport,view && !preview);
+    if(!preview) return F;
+  }
+
+  string image;
+  if(preview && settings.embed) {
+    image=prefix;
+    if(settings.inlinetex) image += "_0";
+    image += "."+nativeformat();
+    if(!settings.inlinetex) file3.push(image);
+    image=graphic(image,"hiresbb");
+  }
+  if(prc) {
+    if(!P.infinity && P.viewportshift != 0)
+      warning("offaxis",
+              "PRC does not support off-axis projections; use pan instead of
+shift");
+    F.L=embed3D(label,text=image,prefix,S.f,format,
+                S.width-2,S.height-2,options,script,light,Q);
+  }
+  return F;
+}
 
-    string image;
-    if(preview && settings.embed) {
-      image=prefix;
-      if(settings.inlinetex) image += "_0";
-      image += "."+nativeformat();
-      if(!settings.inlinetex) file3.push(image);
-      image=graphic(image,"hiresbb");
-    }
-    if(prc) {
-      if(!P.infinity && P.viewportshift != 0)
-        write("warning: PRC does not support off-axis projections; use pan instead of shift");
-      F.L=embed3D(label,text=image,prefix,f,format,
-                  width,height,angle,options,script,light,Q);
-    }
-    
+object embed(string label="", string text=label,
+             string prefix=defaultfilename,
+             picture pic, string format="",
+             real xsize=pic.xsize, real ysize=pic.ysize,
+             bool keepAspect=pic.keepAspect, bool view=true, string options="",
+             string script="", light light=currentlight,
+             projection P=currentprojection)
+{
+  bool is3D=is3D(format);
+  scene S=scene(pic,xsize,ysize,keepAspect,is3D,P);
+  if(is3D)
+    return embed(label,text,prefix,S,format,view,options,script,light);
+  else {
+    object F;
+    transform T=S.pic2.scaling(xsize,ysize,keepAspect);
+    F.f=pic.fit(scale(S.t[0][0])*T);
+    add(F.f,S.pic2.fit(T));
+    return F;
   }
+}
 
-  if(!is3D) {
-    transform T=pic2.scaling(xsize,ysize,keepAspect);
-    F.f=pic.fit(scale(t[0][0])*T);
-    add(F.f,pic2.fit(T));
+object embed(string label="", string text=label,
+             string prefix=defaultfilename,
+             frame f, string format="", real width=0, real height=0,
+             bool view=true, string options="", string script="",
+             light light=currentlight, projection P=currentprojection)
+{
+  if(is3D(format))
+    return embed(label,text,prefix,scene(f,P),format,view,options,script,light);
+  else {
+    object F;
+    F.f=f;
+    return F;
   }
-      
-  return F;
 }
 
 embed3=new object(string prefix, frame f, string format, string options,
-                  string script, projection P) {
-  return embed(prefix=prefix,f,format,options,script,P);
+                  string script, light light, projection P) {
+  return embed(prefix=prefix,f,format,options,script,light,P);
 };
 
+frame embedder(object embedder(string prefix, string format),
+               string prefix, string format, bool view, light light)
+{
+  frame f;
+  bool prc=prc(format);
+  if(!prc && settings.render != 0 && !view) {
+    static int previewcount=0;
+    bool keep=prefix != "";
+    prefix=outprefix(prefix)+"+"+(string) previewcount;
+    ++previewcount;
+    format=nativeformat();
+    if(!keep) file3.push(prefix+"."+format);
+  }
+  object F=embedder(prefix,format);
+  if(prc)
+    label(f,F.L);
+  else {
+    if(settings.render == 0) {
+      add(f,F.f);
+      if(light.background != nullpen)
+        box(f,light.background,Fill,above=false);
+    } else if(!view)
+      label(f,graphic(prefix,"hiresbb"));
+  }
+  return f;
+}
+
 currentpicture.fitter=new frame(string prefix, picture pic, string format,
-                                real xsize, real ysize,
-                                bool keepAspect, bool view,
-                                string options, string script, projection P) {
+                                real xsize, real ysize, bool keepAspect,
+                                bool view, string options, string script,
+                                light light, projection P) {
   frame f;
   bool empty3=pic.empty3();
+  if(!empty3) f=embedder(new object(string prefix, string format) {
+      return embed(prefix=prefix,pic,format,xsize,ysize,keepAspect,view,
+                   options,script,light,P);
+    },prefix,format,view,light);
   if(is3D(format) || empty3) add(f,pic.fit2(xsize,ysize,keepAspect));
-  if(!empty3) {
-    bool prc=prc(format);
-    if(!prc && settings.render != 0 && !view) {
-      static int previewcount=0;
-      bool keep=prefix != "";
-      prefix=outprefix(prefix)+"+"+(string) previewcount;
-      ++previewcount;
-      format=nativeformat();
-      if(!keep) file3.push(prefix+"."+format);
-    }
-    object F=embed(prefix=prefix,pic,format,xsize,ysize,keepAspect,view,
-                   options,script,currentlight,P);
-    if(prc)
-      label(f,F.L);
-    else {
-      if(settings.render == 0) {
-        add(f,F.f);
-        if(currentlight.background != nullpen)
-          box(f,currentlight.background,Fill,above=false);
-      } else if(!view)
-        label(f,graphic(prefix,"hiresbb"));
-    }
-  }
   return f;
 };
 
+frame embedder(string prefix, frame f, string format, bool view,
+               string options, string script, light light, projection P)
+{
+  return embedder(new object(string prefix, string format) {
+      return embed(prefix=prefix,f,format,view,options,script,light,P);
+    },prefix,format,view,light);
+}
+
 void addViews(picture dest, picture src, bool group=true,
               filltype filltype=NoFill)
 {
@@ -2733,29 +2836,53 @@ void addAllViews(picture src, bool group=true, filltype filltype=NoFill)
   addAllViews(currentpicture,src,group,filltype);
 }
 
-// Force an array of 3D pictures to be as least as large as picture all.
-void rescale3(picture[] pictures, picture all, projection P=currentprojection)
+// Fit an array of 3D pictures simultaneously using the sizing of picture all.
+frame[] fit3(string prefix="", picture[] pictures, picture all,
+             string format="", bool view=true, string options="",
+             string script="", light light=currentlight,
+             projection P=currentprojection)
 {
-  if(!all.empty3()) {
-    transform3 t=inverse(all.calculateTransform3(P)*pictures[0].T3);
-    triple m=t*min3(all);
-    triple M=t*max3(all);
+  frame[] out;
+  scene S=scene(all,P);
+  triple m=all.min(S.t);
+  triple M=all.max(S.t);
+  out=new frame[pictures.length];
+  int i=0;
+  bool loop=settings.loop;
+  for(picture pic : pictures) {
+    picture pic2;
+    frame f=pic.fit3(S.t,pic2,S.P);
+    if(loop && !settings.loop) break;
+    add(f,pic2.fit2());
+    draw(f,m,nullpen);
+    draw(f,M,nullpen);
+    out[i]=loop ? f : embedder(prefix,f,format,view,options,script,light,S.P);
+    ++i;
+  }
+
+  while(settings.loop)
     for(int i=0; i < pictures.length; ++i) {
-      draw(pictures[i],m,nullpen);
-      draw(pictures[i],M,nullpen);
+      if(!settings.loop) break;
+      embedder(prefix,out[i],format,view,options,script,light,S.P);
     }
-  }
+  
+  return out;
 }
 
-// Force an array of pictures to have a uniform scaling using currenprojection.
-rescale=new void(picture[] pictures) {
-  if(pictures.length == 0) return;
+// Fit an array of pictures simultaneously using the size of the first picture.
+fit=new frame[](string prefix="", picture[] pictures, string format="",
+                bool view=true, string options="", string script="",
+                projection P=currentprojection) {
+  if(pictures.length == 0)
+    return new frame[];
+
   picture all;
   size(all,pictures[0]);
   for(picture pic : pictures)
     add(all,pic);
-  rescale2(pictures,all);
-  rescale3(pictures,all);
+
+  return all.empty3() ? fit2(pictures,all) :
+  fit3(prefix,pictures,all,format,view,options,script,P);
 };
 
 exitfcn currentexitfunction=atexit();
diff --git a/Master/texmf/asymptote/three_arrows.asy b/Master/texmf/asymptote/three_arrows.asy
index 841c4f731f4..161a716bca6 100644
--- a/Master/texmf/asymptote/three_arrows.asy
+++ b/Master/texmf/asymptote/three_arrows.asy
@@ -48,91 +48,20 @@ void bend(surface s, path3 g, real L)
   }
 }
 
-void render(path3 s, real granularity=tubegranularity, void f(path3, real))
-{
-  static int maxdepth=ceil(-log(realEpsilon)/log(2))+1;
-  void Split(triple z0, triple c0, triple c1, triple z1, real t0=0, real t1=1,
-             real depth=mantissaBits) {
-    if(depth > 0) {
-      real S=straightness(z0,c0,c1,z1);
-      real R=infinity;
-      int nintervals=3;
-      real fuzz=sqrtEpsilon*max(abs(z0),abs(z1));
-      if(S > 0) {
-        --depth;
-        for(int i=0; i <= nintervals; ++i) {
-          R=min(R,radius(z0,c0,c1,z1,i/nintervals));
-          if(S > max(granularity*R,fuzz)) {
-            triple m0=0.5*(z0+c0);
-            triple m1=0.5*(c0+c1);
-            triple m2=0.5*(c1+z1);
-            triple m3=0.5*(m0+m1);
-            triple m4=0.5*(m1+m2);
-            triple m5=0.5*(m3+m4);
-            real tm=0.5*(t0+t1);
-            Split(z0,m0,m3,m5,t0,tm,depth);
-            Split(m5,m4,m2,z1,tm,t1,depth);
-            return;
-          }
-        }
-      }
-    }
-    f(z0..controls c0 and c1..z1,t0);
-  }
-  Split(point(s,0),postcontrol(s,0),precontrol(s,1),point(s,1));
-}
-
-struct tube 
-{
-  surface s;
-  path3 center;
-
-  void operator init(path3 g, real width, real granularity=tubegranularity) {
-    real r=0.5*width;
-
-    transform3 t=scale3(r);
-
-    int n=length(g);
-    for(int i=0; i < n; ++i) {
-      path3 gi=subpath(g,i,i+1);
-      real S=straightness(g,i);
-      if(S < sqrtEpsilon*r) {
-        triple v=point(g,i);
-        triple u=point(g,i+1)-v;
-        s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder);
-        center=center&gi;
-      } else {
-        render(gi,granularity,new void(path3 q, real) {
-            real L=arclength(q);
-            surface segment=scale(r,r,L)*unitcylinder;
-            bend(segment,q,L);
-            s.s.append(segment.s);
-            triple a=L*Z/3;
-            center=center&bend(O,q,L)..controls bend(a,q,L) and bend(2a,q,L)..
-              bend(3a,q,L);
-          });
-      }
-
-      if((cyclic(g) || i > 0) && abs(dir(g,i,1)-dir(g,i,-1)) > sqrtEpsilon)
-        s.append(shift(point(g,i))*t*align(dir(g,i,-1))*unithemisphere);
-    }
-  }
-}
-
 // Refine a noncyclic path3 g so that it approaches its endpoint in
 // geometrically spaced steps.
-path3 approach(path3 g, int n, real radix=3)
+path3 approach(path3 p, int n, real radix=3)
 {
   guide3 G;
-  real L=length(g);
+  real L=length(p);
   real tlast=0;
   real r=1/radix;
   for(int i=1; i < n; ++i) {
     real t=L*(1-r^i);
-    G=G&subpath(g,tlast,t);
+    G=G&subpath(p,tlast,t);
     tlast=t;
   }
-  return G&subpath(g,tlast,L);
+  return G&subpath(p,tlast,L);
 }
 
 struct arrowhead3
@@ -176,7 +105,7 @@ struct arrowhead3
     }
     if(draw)
       for(path3 g : H)
-        s.append(tube(g,width).s);
+        s.append(tube(g,width,linesectors).s);
     return shift(v)*s;
   }
 
@@ -231,8 +160,7 @@ DefaultHead3.head=new surface(path3 g, position position=EndPoint,
           for(int i=0; i < 4; ++i) {
             for(int j=0; j < 4; ++j) {
               real k=1-p.P[i][j].z/remainL;
-              p.P[i][j]=(k*p.P[i][j].x,k*p.P[i][j].y,p.P[i][j].z);
-              p.P[i][j]=bend(p.P[i][j],q,l);
+              p.P[i][j]=bend((k*p.P[i][j].x,k*p.P[i][j].y,p.P[i][j].z),q,l);
             }
           }
         }
diff --git a/Master/texmf/asymptote/three_light.asy b/Master/texmf/asymptote/three_light.asy
index 8e5a6281795..36ba9a34c19 100644
--- a/Master/texmf/asymptote/three_light.asy
+++ b/Master/texmf/asymptote/three_light.asy
@@ -75,88 +75,27 @@ material emissive(material m, real granularity=m.granularity)
                   granularity);
 }
 
-struct light {
-  real[][] diffuse;
-  real[][] ambient;
-  real[][] specular;
-  pen background; // Background color of the 3D canvas.
-  real specularfactor;
-  bool viewport; // Are the lights specified (and fixed) in the viewport frame?
-  triple[] position; // Only directional lights are currently implemented.
-
-  transform3 T=identity(4); // Transform to apply to normal vectors.
-
-  bool on() {return position.length > 0;}
-  
-  void operator init(pen[] diffuse,
-                     pen[] ambient=array(diffuse.length,black),
-                     pen[] specular=diffuse, pen background=nullpen,
-                     real specularfactor=1,
-                     bool viewport=true, triple[] position) {
-    int n=diffuse.length;
-    assert(ambient.length == n && specular.length == n && position.length == n);
-    
-    this.diffuse=new real[n][];
-    this.ambient=new real[n][];
-    this.specular=new real[n][];
-    this.background=background;
-    this.position=new triple[n];
-    for(int i=0; i < position.length; ++i) {
-      this.diffuse[i]=rgba(diffuse[i]);
-      this.ambient[i]=rgba(ambient[i]);
-      this.specular[i]=rgba(specular[i]);
-      this.position[i]=unit(position[i]);
-    }
-    this.specularfactor=specularfactor;
-    this.viewport=viewport;
-  }
-
-  void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
-                     pen background=nullpen, real specularfactor=1,
-                     bool viewport=true...triple[] position) {
-    int n=position.length;
-    operator init(array(n,diffuse),array(n,ambient),array(n,specular),
-                  background,specularfactor,viewport,position);
-  }
-
-  void operator init(pen diffuse=white, pen ambient=black, pen specular=diffuse,
-                     pen background=nullpen, bool viewport=true,
-                     real x, real y, real z) {
-    operator init(diffuse,ambient,specular,background,viewport,(x,y,z));
+pen color(triple normal, material m, light light, transform3 T=light.T) {
+  triple[] position=light.position;
+  if(invisible((pen) m)) return invisible;
+  if(position.length == 0) return m.diffuse();
+  normal=unit(T*normal);
+  if(settings.twosided) normal *= sgn(normal.z);
+  real s=m.shininess*128;
+  real[] Diffuse=rgba(m.diffuse());
+  real[] Ambient=rgba(m.ambient());
+  real[] Specular=rgba(m.specular());
+  real[] p=rgba(m.emissive());
+  for(int i=0; i < position.length; ++i) {
+    triple L=light.viewport ? position[i] : T*position[i];
+    real Ldotn=max(dot(normal,L),0);
+    p += light.ambient[i]*Ambient+Ldotn*light.diffuse[i]*Diffuse;
+    // Apply specularfactor to partially compensate non-pixel-based rendering.
+    if(Ldotn > 0) // Phong-Blinn model of specular reflection
+      p += dot(normal,unit(L+Z))^s*light.specularfactor*
+        light.specular[i]*Specular;
   }
-
-  void operator init(explicit light light) {
-    diffuse=copy(light.diffuse);
-    ambient=copy(light.ambient);
-    specular=copy(light.specular);
-    background=light.background;
-    specularfactor=light.specularfactor;
-    viewport=light.viewport;
-    position=copy(light.position);
-  }
-
-  pen color(triple normal, material m, transform3 T=T) {
-    if(invisible((pen) m)) return invisible;
-    if(position.length == 0) return m.diffuse();
-    normal=unit(T*normal);
-    if(settings.twosided) normal *= sgn(normal.z);
-    real s=m.shininess*128;
-    real[] Diffuse=rgba(m.diffuse());
-    real[] Ambient=rgba(m.ambient());
-    real[] Specular=rgba(m.specular());
-    real[] p=rgba(m.emissive());
-    for(int i=0; i < position.length; ++i) {
-      triple L=viewport ? position[i] : T*position[i];
-      real Ldotn=max(dot(normal,L),0);
-      p += ambient[i]*Ambient+Ldotn*diffuse[i]*Diffuse;
-      // Apply specularfactor to partially compensate non-pixel-based rendering.
-      if(Ldotn > 0) // Phong-Blinn model of specular reflection
-        p += dot(normal,unit(L+Z))^s*specularfactor*specular[i]*Specular;
-    }
-    return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse()));
-  }
-
-  real[] background() {return rgba(background == nullpen ? white : background);}
+  return rgb(p[0],p[1],p[2])+opacity(opacity(m.diffuse()));
 }
 
 light operator * (transform3 t, light light)
@@ -168,14 +107,16 @@ light operator * (transform3 t, light light)
 
 light operator cast(triple v) {return light(v);}
 
-light currentlight=light(ambient=gray(0.1),specularfactor=3,
-                         (0.25,-0.25,1));
+light Viewport=light(ambient=gray(0.1),specularfactor=3,viewport=true,
+                     (0.25,-0.25,1));
 
 light White=light(new pen[] {rgb(0.38,0.38,0.45),rgb(0.6,0.6,0.67),
-                             rgb(0.5,0.5,0.57)},specularfactor=3,viewport=false,
+                             rgb(0.5,0.5,0.57)},specularfactor=3,
   new triple[] {(-2,-1.5,-0.5),(2,1.1,-2.5),(-0.5,0,2)});
 
-light Headlamp=light(gray(0.9),ambient=gray(0.1),specularfactor=3,
-                     (0.5,0.5,1/sqrt(2)),specular=gray(0.7));
+light Headlamp=light(gray(0.8),ambient=gray(0.1),specular=gray(0.7),
+                     specularfactor=3,viewport=true,dir(42,48));
+
+currentlight=Headlamp;
 
 light nolight;
diff --git a/Master/texmf/asymptote/three_surface.asy b/Master/texmf/asymptote/three_surface.asy
index f2ad86dfe7e..e545d87cf7b 100644
--- a/Master/texmf/asymptote/three_surface.asy
+++ b/Master/texmf/asymptote/three_surface.asy
@@ -1,4 +1,5 @@
 import bezulate;
+private import interpolate;
 
 int nslice=12;
 real camerafactor=1.2;
@@ -15,24 +16,24 @@ struct patch {
 
   path3 external() {
     return
-      P[0][0]..controls P[0][1] and P[0][2]..
-      P[0][3]..controls P[1][3] and P[2][3]..
-      P[3][3]..controls P[3][2] and P[3][1]..
-      P[3][0]..controls P[2][0] and P[1][0]..cycle;
+      P[0][0]..controls P[1][0] and P[2][0]..
+      P[3][0]..controls P[3][1] and P[3][2]..
+      P[3][3]..controls P[2][3] and P[1][3]..
+      P[0][3]..controls P[0][2] and P[0][1]..cycle;
   }
 
   triple[] internal() {
-    return new triple[] {P[1][1],P[1][2],P[2][2],P[2][1]};
+    return new triple[] {P[1][1],P[2][1],P[2][2],P[1][2]};
   }
 
   triple cornermean() {
-    return 0.25*(P[0][0]+P[0][3]+P[3][3]+P[3][0]);
+    return 0.25*(P[0][0]+P[0][3]+P[3][0]+P[3][3]);
   }
 
-  triple[] corners() {return new triple[] {P[0][0],P[0][3],P[3][3],P[3][0]};}
+  triple[] corners() {return new triple[] {P[0][0],P[3][0],P[3][3],P[0][3]};}
 
   real[] map(real f(triple)) {
-    return new real[] {f(P[0][0]),f(P[0][3]),f(P[3][3]),f(P[3][0])};
+    return new real[] {f(P[0][0]),f(P[3][0]),f(P[3][3]),f(P[0][3])};
   }
 
   triple Bu(int j, real u) {return bezier(P[0][j],P[1][j],P[2][j],P[3][j],u);}
@@ -43,8 +44,10 @@ struct patch {
   triple BuPPP(int j) {return bezierPPP(P[0][j],P[1][j],P[2][j],P[3][j]);}
 
   path3 uequals(real u) {
-    return straight ? Bu(0,u)--Bu(3,u) :
-      Bu(0,u)..controls Bu(1,u) and Bu(2,u)..Bu(3,u);
+    triple z0=Bu(0,u);
+    triple z1=Bu(3,u);
+    return path3(new triple[] {z0,Bu(2,u)},new triple[] {z0,z1},
+                 new triple[] {Bu(1,u),z1},new bool[] {straight,false},false);
   }
 
   triple Bv(int i, real v) {return bezier(P[i][0],P[i][1],P[i][2],P[i][3],v);}
@@ -55,8 +58,10 @@ struct patch {
   triple BvPPP(int i) {return bezierPPP(P[i][0],P[i][1],P[i][2],P[i][3]);}
 
   path3 vequals(real v) {
-    return straight ? Bv(0,v)--Bv(3,v) :
-      Bv(0,v)..controls Bv(1,v) and Bv(2,v)..Bv(3,v);
+    triple z0=Bv(0,v);
+    triple z1=Bv(3,v);
+    return path3(new triple[] {z0,Bv(2,v)},new triple[] {z0,z1},
+                 new triple[] {Bv(1,v),z1},new bool[] {straight,false},false);
   }
 
   triple point(real u, real v) {        
@@ -65,76 +70,76 @@ struct patch {
 
   // compute normal vectors for degenerate cases
   private triple normal0(real u, real v, real epsilon) {
-    triple n=0.5*(cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),
-                        bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+
-                  cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),   
-                        bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v)));
+    triple n=0.5*(cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v),
+                        bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u))+
+                  cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v),   
+                        bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u)));
     return abs(n) > epsilon ? n :
-      0.25*cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),   
-                 bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+
-      1/6*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),
-                 bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))+
-           cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),   
-                 bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+
-      1/12*(cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),
-                  bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v))+
-            cross(bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u),   
-                  bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v)))+
-      1/36*cross(bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u),   
-                 bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v));
+      0.25*cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v),   
+                 bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u))+
+      1/6*(cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v),   
+                 bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u))+
+           cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v),
+                 bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u)))+
+      1/12*(cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v),
+                  bezier(BvPP(0,v),BvPP(1,v),BvPP(2,v),BvPP(3,v),u))+
+            cross(bezier(BuPP(0,u),BuPP(1,u),BuPP(2,u),BuPP(3,u),v),   
+                  bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u)))+
+      1/36*cross(bezier(BuPPP(0),BuPPP(1),BuPPP(2),BuPPP(3),v),   
+                 bezier(BvPPP(0),BvPPP(1),BvPPP(2),BvPPP(3),u));
   }
 
   static real fuzz=1000*realEpsilon;
 
   triple normal(real u, real v) {
-    triple n=cross(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u),
-                   bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v));
+    triple n=cross(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v),
+                   bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u));
     real epsilon=fuzz*change2(P);
     return (abs(n) > epsilon) ? n : normal0(u,v,epsilon);
   }
   
   triple normal00() {
-    triple n=9*cross(P[0][1]-P[0][0],P[1][0]-P[0][0]);
+    triple n=9*cross(P[1][0]-P[0][0],P[0][1]-P[0][0]);
     real epsilon=fuzz*change2(P);
     return abs(n) > epsilon ? n : normal0(0,0,epsilon);
   }
 
-  triple normal01() {
-    triple n=9*cross(P[0][3]-P[0][2],P[1][3]-P[0][3]);
+  triple normal10() {
+    triple n=9*cross(P[3][0]-P[2][0],P[3][1]-P[3][0]);
     real epsilon=fuzz*change2(P);
-    return abs(n) > epsilon ? n : normal0(0,1,epsilon);
+    return abs(n) > epsilon ? n : normal0(1,0,epsilon);
   }
 
   triple normal11() {
-    triple n=9*cross(P[3][3]-P[3][2],P[3][3]-P[2][3]);
+    triple n=9*cross(P[3][3]-P[2][3],P[3][3]-P[3][2]);
     real epsilon=fuzz*change2(P);
     return abs(n) > epsilon ? n : normal0(1,1,epsilon);
   }
 
-  triple normal10() {
-    triple n=9*cross(P[3][1]-P[3][0],P[3][0]-P[2][0]);
+  triple normal01() {
+    triple n=9*cross(P[1][3]-P[0][3],P[0][3]-P[0][2]);
     real epsilon=fuzz*change2(P);
-    return abs(n) > epsilon ? n : normal0(1,0,epsilon);
+    return abs(n) > epsilon ? n : normal0(0,1,epsilon);
   }
 
   pen[] colors(material m, light light=currentlight) {
     bool nocolors=colors.length == 0;
     if(normals.length > 0)
-      return new pen[] {light.color(normals[0],nocolors ? m : colors[0]),
-          light.color(normals[1],nocolors ? m : colors[1]),
-          light.color(normals[2],nocolors ? m : colors[2]),
-          light.color(normals[3],nocolors ? m : colors[3])};
+      return new pen[] {color(normals[0],nocolors ? m : colors[0],light),
+          color(normals[1],nocolors ? m : colors[1],light),
+          color(normals[2],nocolors ? m : colors[2],light),
+          color(normals[3],nocolors ? m : colors[3],light)};
     if(planar) {
       triple normal=normal(0.5,0.5);
-      return new pen[] {light.color(normal,nocolors ? m : colors[0]),
-          light.color(normal,nocolors ? m : colors[1]),
-          light.color(normal,nocolors ? m : colors[2]),
-          light.color(normal,nocolors ? m : colors[3])};
+      return new pen[] {color(normal,nocolors ? m : colors[0],light),
+          color(normal,nocolors ? m : colors[1],light),
+          color(normal,nocolors ? m : colors[2],light),
+          color(normal,nocolors ? m : colors[3],light)};
     }
-    return new pen[] {light.color(normal00(),nocolors ? m : colors[0]),
-        light.color(normal01(),nocolors ? m : colors[1]),
-        light.color(normal11(),nocolors ? m : colors[2]),
-        light.color(normal10(),nocolors ? m : colors[3])};
+    return new pen[] {color(normal00(),nocolors ? m : colors[0],light),
+        color(normal10(),nocolors ? m : colors[1],light),
+        color(normal11(),nocolors ? m : colors[2],light),
+        color(normal01(),nocolors ? m : colors[3],light)};
   }
   
   triple min3,max3;
@@ -162,11 +167,19 @@ struct patch {
   }
 
   pair min(projection P, pair bound=project(this.P[0][0],P.t)) {
-    return minbezier(this.P,P.t,bound);
+    triple[][] Q=P.T.modelview*this.P;
+    if(P.infinity)
+      return xypart(minbezier(Q,(bound.x,bound.y,0)));
+    real d=P.T.projection[3][2];
+    return maxratio(Q,d*bound)/d; // d is negative
   }
 
   pair max(projection P, pair bound=project(this.P[0][0],P.t)) {
-    return maxbezier(this.P,P.t,bound);
+    triple[][] Q=P.T.modelview*this.P;
+    if(P.infinity)
+      return xypart(maxbezier(Q,(bound.x,bound.y,0)));
+    real d=P.T.projection[3][2];
+    return minratio(Q,d*bound)/d; // d is negative
   }
 
   void operator init(triple[][] P, triple[] normals=new triple[],
@@ -242,11 +255,11 @@ struct patch {
     } else straight=false;
 
     P=new triple[][] {
-      {point(external,0),postcontrol(external,0),precontrol(external,1),
-       point(external,1)},
-      {precontrol(external,0),internal[0],internal[1],postcontrol(external,1)},
-      {postcontrol(external,3),internal[3],internal[2],precontrol(external,2)},
-      {point(external,3),precontrol(external,3),postcontrol(external,2),
+      {point(external,0),precontrol(external,0),postcontrol(external,3),
+       point(external,3)},
+      {postcontrol(external,0),internal[0],internal[3],precontrol(external,3)},
+      {precontrol(external,1),internal[1],internal[2],postcontrol(external,2)},
+      {point(external,1),postcontrol(external,1),precontrol(external,2),
        point(external,2)}
     };
   }
@@ -280,10 +293,10 @@ struct patch {
       delta[j]=(external[(j+1)% 4]-external[j])/3;
 
     P=new triple[][] {
-      {external[0],external[0]+delta[0],external[1]-delta[0],external[1]},
-      {external[0]-delta[3],internal[0],internal[1],external[1]+delta[1]},
-      {external[3]+delta[3],internal[3],internal[2],external[2]-delta[1]},
-      {external[3],external[3]-delta[2],external[2]+delta[2],external[2]}
+      {external[0],external[0]-delta[3],external[3]+delta[3],external[3]},
+      {external[0]+delta[0],internal[0],internal[3],external[3]-delta[2]},
+      {external[1]-delta[0],internal[1],internal[2],external[2]+delta[2]},
+      {external[1],external[1]+delta[1],external[2]-delta[1],external[2]}
     };
   }
 }
@@ -322,8 +335,243 @@ patch reverse(patch s)
   return S;
 }
 
+// Return the Coons patch control points corresponding to path p.
+pair[][] coons(path p)
+{
+  int L=length(p);
+  if(L == 1)
+    p=p--cycle--cycle--cycle;
+  else if(L == 2)
+    p=p--cycle--cycle;
+  else if(L == 3)
+    p=p--cycle;
+ 
+  pair[] internal=new pair[4];
+  for(int j=0; j < 4; ++j) {
+    internal[j]=nineth*(-4*point(p,j)
+                        +6*(precontrol(p,j)+postcontrol(p,j))
+                        -2*(point(p,j-1)+point(p,j+1))
+                        +3*(precontrol(p,j-1)+postcontrol(p,j+1))
+                        -point(p,j+2));
+  }
+    
+  return new pair[][] {
+    {point(p,0),precontrol(p,0),postcontrol(p,3),point(p,3)},
+      {postcontrol(p,0),internal[0],internal[3],precontrol(p,3)},
+        {precontrol(p,1),internal[1],internal[2],postcontrol(p,2)},
+          {point(p,1),postcontrol(p,1),precontrol(p,2),point(p,2)}
+  };
+}
+
+// Decompose a possibly nonconvex cyclic path into an array of paths that
+// yield nondegenerate Coons patches.
+path[] regularize(path p, bool checkboundary=true)
+{
+  path[] s;
+
+  if(!cyclic(p))
+    abort("cyclic path expected");
+
+  int L=length(p);
+
+  if(L > 4) {
+    for(path g : bezulate(p))
+      s.append(regularize(g,checkboundary));
+    return s;
+  }
+        
+  bool straight=piecewisestraight(p);
+  if(L <= 3 && straight) {
+    return new path[] {p};
+  }
+    
+  // Split p along the angle bisector at t.
+  bool split(path p, real t) {
+    pair dir=dir(p,t);
+    if(dir != 0) {
+      path g=subpath(p,t,t+length(p));
+      int L=length(g);
+      pair z=point(g,0);
+      real[] T=intersections(g,z,z+I*dir);
+      for(int i=0; i < T.length; ++i) {
+        real cut=T[i];
+        if(cut > sqrtEpsilon && cut < L-sqrtEpsilon) {
+          pair w=point(g,cut);
+          if(!inside(p,0.5*(z+w),zerowinding)) continue;
+          pair delta=sqrtEpsilon*(w-z);
+          if(intersections(g,z-delta--w+delta).length != 2) continue;
+          s.append(regularize(subpath(g,0,cut)--cycle,checkboundary));
+          s.append(regularize(subpath(g,cut,L)--cycle,checkboundary));
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+
+  // Ensure that all interior angles are less than 180 degrees.
+  real fuzz=1e-4;
+  int sign=sgn(windingnumber(p,inside(p,zerowinding)));
+  for(int i=0; i < L; ++i) {
+    if(sign*(conj(dir(p,i,-1))*dir(p,i,1)).y < -fuzz) {
+      if(split(p,i)) return s;
+    }
+  }
+
+  if(straight)
+    return new path[] {p};
+    
+  pair[][] P=coons(p);
+
+  // Check for degeneracy.
+  pair[][] U=new pair[3][4];
+  pair[][] V=new pair[4][3];
+
+  for(int i=0; i < 3; ++i) {
+    for(int j=0; j < 4; ++j)
+      U[i][j]=P[i+1][j]-P[i][j];
+  }
+      
+  for(int i=0; i < 4; ++i) {
+    for(int j=0; j < 3; ++j)
+      V[i][j]=P[i][j+1]-P[i][j];
+  }
+      
+  int[] choose2={1,2,1};
+  int[] choose3={1,3,3,1};
+
+  real T[][]=new real[6][6];
+  for(int p=0; p < 6; ++p) {
+    int kstart=max(p-2,0);
+    int kstop=min(p,3);
+    real[] Tp=T[p];
+    for(int q=0; q < 6; ++q) {
+      real Tpq;
+      int jstop=min(q,3);
+      int jstart=max(q-2,0);
+      for(int k=kstart; k <= kstop; ++k) {
+        int choose3k=choose3[k];
+        for(int j=jstart; j <= jstop; ++j) {
+          int i=p-k;
+          int l=q-j;
+          Tpq += (conj(U[i][j])*V[k][l]).y*
+            choose2[i]*choose3k*choose3[j]*choose2[l];
+        }
+      }
+      Tp[q]=Tpq;
+    }
+  }
+
+  bool3 aligned=default;
+  bool degenerate=false;
+
+  for(int p=0; p < 6; ++p) {
+    for(int q=0; q < 6; ++q) {
+      if(aligned == default) {
+        if(T[p][q] > sqrtEpsilon) aligned=true;
+        if(T[p][q] < -sqrtEpsilon) aligned=false;
+      } else {
+        if((T[p][q] > sqrtEpsilon && aligned == false) ||
+           (T[p][q] < -sqrtEpsilon && aligned == true)) degenerate=true;
+      }
+    }
+  }
+
+  if(!degenerate) {
+    if(aligned == (sign >= 0))
+      return new path[] {p};
+    return s;
+  }
+
+  if(checkboundary) {
+    // Polynomial coefficients of (B_i'' B_j + B_i' B_j')/3.
+    static real[][][] fpv0={
+      {{5, -20, 30, -20, 5},
+       {-3, 24, -54, 48, -15},
+       {0, -6, 27, -36, 15},
+       {0, 0, -3, 8, -5}},
+      {{-7, 36, -66, 52, -15},
+       {3, -36, 108, -120, 45},
+       {0, 6, -45, 84, -45},
+       {0, 0, 3, -16, 15}},
+      {{2, -18, 45, -44, 15},
+       {0, 12, -63, 96, -45},
+       {0, 0, 18, -60, 45},
+       {0, 0, 0, 8, -15}},
+      {{0, 2, -9, 12, -5},
+       {0, 0, 9, -24, 15},
+       {0, 0, 0, 12, -15},
+       {0, 0, 0, 0, 5}}
+    };
+
+    // Compute one-ninth of the derivative of the Jacobian along the boundary.
+    real[][] c=array(4,array(5,0.0));
+    for(int i=0; i < 4; ++i) {
+      real[][] fpv0i=fpv0[i];
+      for(int j=0; j < 4; ++j) {
+        real[] w=fpv0i[j];
+        c[0] += w*(conj(P[i][0])*(P[j][1]-P[j][0])).y; // v=0
+        c[1] += w*(conj(P[3][j]-P[2][j])*P[3][i]).y;   // u=1
+        c[2] += w*(conj(P[i][3])*(P[j][3]-P[j][2])).y; // v=1
+        c[3] += w*(conj(P[0][j]-P[1][j])*P[0][i]).y;   // u=0
+      }
+    }
+    
+    pair BuP(int j, real u) {
+      return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u);
+    }
+    pair BvP(int i, real v) {
+      return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v);
+    }
+    real normal(real u, real v) {
+      return (conj(bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v))*
+              bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u)).y;
+    }
+
+    // Use Rolle's theorem to check for degeneracy on the boundary.
+    real M=0;
+    real cut;
+    for(int i=0; i < 4; ++i) {
+      if(!straight(p,i)) {
+        real[] ci=c[i];
+        pair[] R=quarticroots(ci[4],ci[3],ci[2],ci[1],ci[0]);
+        for(pair r : R) {
+          if(fabs(r.y) < sqrtEpsilon) {
+            real t=r.x;
+            if(0 <= t && t <= 1) {
+              real[] U={t,1,t,0};
+              real[] V={0,t,1,t};
+              real[] T={t,t,1-t,1-t};
+              real N=sign*normal(U[i],V[i]);
+              if(N < M) {
+                M=N; cut=i+T[i];
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // Split at the worst boundary degeneracy.
+    if(M < 0 && split(p,cut)) return s;
+  }
+    
+  // Split arbitrarily to resolve any remaining (internal) degeneracy.
+  checkboundary=false;
+  for(int i=0; i < L; ++i)
+    if(!straight(p,i) && split(p,i+0.5)) return s;
+
+  while(true)
+    for(int i=0; i < L; ++i)
+      if(!straight(p,i) && split(p,i+unitrand())) return s;
+
+  return s;
+}
+
 struct surface {
   patch[] s;
+  int index[][];
+  bool vcyclic;
   
   bool empty() {
     return s.length == 0;
@@ -341,6 +589,8 @@ struct surface {
     this.s=new patch[s.s.length];
     for(int i=0; i < s.s.length; ++i)
       this.s[i]=patch(s.s[i]);
+    this.index=copy(s.index);
+    this.vcyclic=s.vcyclic;
   }
 
   void operator init(triple[][][] P, triple[][] normals=new triple[][],
@@ -352,10 +602,8 @@ struct surface {
   }
 
   void colors(pen[][] palette) {
-    for(int i=0; i < s.length; ++i) {
-      pen[] palettei=palette[i];
-      s[i].colors=new pen[] {palettei[0],palettei[1],palettei[2],palettei[3]};
-    }
+    for(int i=0; i < s.length; ++i)
+      s[i].colors=copy(palette[i]);
   }
 
   triple[][] corners() {
@@ -376,232 +624,66 @@ struct surface {
     return sequence(new triple(int i) {return s[i].cornermean();},s.length);
   }
 
-  // A constructor for a possibly nonconvex cyclic path in a given plane.
-  void operator init (path p, triple plane(pair)=XYplane,
-                      bool checkboundary=true) {
-    if(!cyclic(p))
-      abort("cyclic path expected");
-
-    int L=length(p);
+  triple point(real u, real v) {
+    int U=floor(u);
+    int V=floor(v);
+    int index=index.length == 0 ? U+V : index[U][V];
+    return s[index].point(u-U,v-V);
+  }    
 
-    if(L > 4) {
-      for(path g : bezulate(p))
-        s.append(surface(g,plane,checkboundary).s);
-      return;
-    }
-        
-    pair[][] P(path p) {
-      if(L == 1)
-        p=p--cycle--cycle--cycle;
-      else if(L == 2)
-        p=p--cycle--cycle;
-      else if(L == 3)
-        p=p--cycle;
- 
-      pair[] internal=new pair[4];
-      for(int j=0; j < 4; ++j) {
-        internal[j]=nineth*(-4*point(p,j)
-                            +6*(precontrol(p,j)+postcontrol(p,j))
-                            -2*(point(p,j-1)+point(p,j+1))
-                            +3*(precontrol(p,j-1)+postcontrol(p,j+1))
-                            -point(p,j+2));
-      }
-    
-      return new pair[][] {
-        {point(p,0),postcontrol(p,0),precontrol(p,1),point(p,1)},
-          {precontrol(p,0),internal[0],internal[1],postcontrol(p,1)},
-            {postcontrol(p,3),internal[3],internal[2],precontrol(p,2)},
-              {point(p,3),precontrol(p,3),postcontrol(p,2),point(p,2)}
-      };
-    }
+  triple normal(real u, real v) {
+    int U=floor(u);
+    int V=floor(v);
+    int index=index.length == 0 ? U+V : index[U][V];
+    return s[index].normal(u-U,v-V);
+  }
+  
+  void ucyclic(bool f) 
+  {
+    index.cyclic=f;
+  }
+  
+  void vcyclic(bool f) 
+  {
+    for(int[] i : index)
+      i.cyclic=f;
+    vcyclic=f;
+  }
+  
+  bool ucyclic() 
+  {
+    return index.cyclic;
+  }
+  
+  bool vcyclic() 
+  {
+    return vcyclic;
+  }
 
+  path3 uequals(real u) {
+    if(index.length == 0) return nullpath3;
+    int U=floor(u);
+    int[] index=index[U];
+    path3 g;
+    for(int i : index)
+      g=g&s[i].uequals(u-U);
+    return vcyclic() ? g&cycle : g;
+  }
+  
+  path3 vequals(real v) {
+    if(index.length == 0) return nullpath3;
+    int V=floor(v);
+    path3 g;
+    for(int[] i : index)
+      g=g&s[i[V]].vequals(v-V);
+    return ucyclic() ? g&cycle : g;
+  }
+  
+  // A constructor for a possibly nonconvex cyclic path in a given plane.
+  void operator init(path p, triple plane(pair)=XYplane) {
     bool straight=piecewisestraight(p);
-    if(L <= 3 && straight) {
-      s=new patch[] {patch(P(p),plane,straight)};
-      return;
-    }
-    
-    // Split p along the angle bisector at t.
-    bool split(path p, real t) {
-      pair dir=dir(p,t);
-      if(dir != 0) {
-        path g=subpath(p,t,t+length(p));
-        int L=length(g);
-        pair z=point(g,0);
-        real[] T=intersections(g,z,z+I*dir);
-        for(int i=0; i < T.length; ++i) {
-          real cut=T[i];
-          if(cut > sqrtEpsilon && cut < L-sqrtEpsilon) {
-            pair w=point(g,cut);
-            if(!inside(p,0.5*(z+w),zerowinding)) continue;
-            pair delta=sqrtEpsilon*(w-z);
-            if(intersections(g,z-delta--w+delta).length != 2) continue;
-            s=surface(subpath(g,0,cut)--cycle,plane,checkboundary).s;
-            s.append(surface(subpath(g,cut,L)--cycle,plane,checkboundary).s);
-            return true;
-          }
-        }
-      }
-      return false;
-    }
-
-    // Ensure that all interior angles are less than 180 degrees.
-    real fuzz=1e-4;
-    int sign=sgn(windingnumber(p,inside(p,zerowinding)));
-    for(int i=0; i < L; ++i) {
-      if(sign*(conj(dir(p,i,-1))*dir(p,i,1)).y < -fuzz) {
-        if(split(p,i)) return;
-      }
-    }
-
-    pair[][] P=P(p);
-
-    if(straight) {
-      s=new patch[] {patch(P,plane,straight)};
-      return;
-    }
-    
-    // Check for degeneracy.
-    pair[][] U=new pair[3][4];
-    pair[][] V=new pair[4][3];
-
-    for(int i=0; i < 3; ++i) {
-      for(int j=0; j < 4; ++j)
-        U[i][j]=P[i+1][j]-P[i][j];
-    }
-      
-    for(int i=0; i < 4; ++i) {
-      for(int j=0; j < 3; ++j)
-        V[i][j]=P[i][j+1]-P[i][j];
-    }
-      
-    int[] choose2={1,2,1};
-    int[] choose3={1,3,3,1};
-
-    real T[][]=new real[6][6];
-    for(int p=0; p < 6; ++p) {
-      int kstart=max(p-2,0);
-      int kstop=min(p,3);
-      real[] Tp=T[p];
-      for(int q=0; q < 6; ++q) {
-        real Tpq;
-        int jstop=min(q,3);
-        int jstart=max(q-2,0);
-        for(int k=kstart; k <= kstop; ++k) {
-          int choose3k=choose3[k];
-          for(int j=jstart; j <= jstop; ++j) {
-            int i=p-k;
-            int l=q-j;
-            Tpq += (conj(U[i][j])*V[k][l]).y*
-              choose2[i]*choose3k*choose3[j]*choose2[l];
-          }
-        }
-        Tp[q]=Tpq;
-      }
-    }
-
-    bool3 aligned=default;
-    bool degenerate=false;
-
-    for(int p=0; p < 6; ++p) {
-      for(int q=0; q < 6; ++q) {
-        if(aligned == default) {
-          if(T[p][q] < -sqrtEpsilon) aligned=true;
-          if(T[p][q] > sqrtEpsilon) aligned=false;
-        } else {
-          if((T[p][q] < -sqrtEpsilon && aligned == false) ||
-             (T[p][q] > sqrtEpsilon && aligned == true)) degenerate=true;
-        }
-      }
-    }
-
-    if(!degenerate) {
-      if(aligned == (sign >= 0))
-        s=new patch[] {patch(P,plane)};
-      return;
-    }
-
-    if(checkboundary) {
-      // Polynomial coefficients of (B_i'' B_j + B_i' B_j')/3.
-      static real[][][] fpv0={
-        {{5, -20, 30, -20, 5},
-         {-3, 24, -54, 48, -15},
-         {0, -6, 27, -36, 15},
-         {0, 0, -3, 8, -5}},
-        {{-7, 36, -66, 52, -15},
-         {3, -36, 108, -120, 45},
-         {0, 6, -45, 84, -45},
-         {0, 0, 3, -16, 15}},
-        {{2, -18, 45, -44, 15},
-         {0, 12, -63, 96, -45},
-         {0, 0, 18, -60, 45},
-         {0, 0, 0, 8, -15}},
-        {{0, 2, -9, 12, -5},
-         {0, 0, 9, -24, 15},
-         {0, 0, 0, 12, -15},
-         {0, 0, 0, 0, 5}}
-      };
-
-      // Compute one-ninth of the derivative of the Jacobian along the boundary.
-      real[][] c=array(4,array(5,0.0));
-      for(int i=0; i < 4; ++i) {
-        real[][] fpv0i=fpv0[i];
-        for(int j=0; j < 4; ++j) {
-          real[] w=fpv0i[j];
-          c[0] += w*(conj(P[0][j]-P[1][j])*P[0][i]).y;   // u=0
-          c[1] += w*(conj(P[i][3])*(P[j][3]-P[j][2])).y; // v=1
-          c[2] += w*(conj(P[3][j]-P[2][j])*P[3][i]).y;   // u=1
-          c[3] += w*(conj(P[i][0])*(P[j][1]-P[j][0])).y; // v=0
-        }
-      }
-    
-      pair BuP(int j, real u) {
-        return bezierP(P[0][j],P[1][j],P[2][j],P[3][j],u);
-      }
-      pair BvP(int i, real v) {
-        return bezierP(P[i][0],P[i][1],P[i][2],P[i][3],v);
-      }
-      real normal(real u, real v) {
-        return (conj(bezier(BvP(0,v),BvP(1,v),BvP(2,v),BvP(3,v),u))*
-                bezier(BuP(0,u),BuP(1,u),BuP(2,u),BuP(3,u),v)).y;
-      }
-
-      // Use Rolle's theorem to check for degeneracy on the boundary.
-      real M=0;
-      real cut;
-      for(int i=0; i < 4; ++i) {
-        if(!straight(p,i)) {
-          real[] ci=c[i];
-          pair[] R=quarticroots(ci[4],ci[3],ci[2],ci[1],ci[0]);
-          for(pair r : R) {
-            if(fabs(r.y) < sqrtEpsilon) {
-              real t=r.x;
-              if(0 <= t && t <= 1) {
-                real[] U={0,t,1,t};
-                real[] V={t,1,t,0};
-                real[] T={t,t,1-t,1-t};
-                real N=sign*normal(U[i],V[i]);
-                if(N < M) {
-                  M=N; cut=i+T[i];
-                }
-              }
-            }
-          }
-        }
-      }
-
-      // Split at the worst boundary degeneracy.
-      if(M < 0 && split(p,cut)) return;
-    }
-    
-    // Split arbitrarily to resolve any remaining (internal) degeneracy.
-    checkboundary=false;
-    for(int i=0; i < L; ++i)
-      if(!straight(p,i) && split(p,i+0.5)) return;
-
-    while(true)
-      for(int i=0; i < L; ++i)
-        if(!straight(p,i) && split(p,i+unitrand())) return;
+    for(path g : regularize(p))
+      s.push(patch(coons(g),plane,straight));
   }
 
   void operator init(explicit path[] g, triple plane(pair)=XYplane) {
@@ -715,12 +797,13 @@ struct surface {
   // An optional surface pen color(int i, real j) may be specified
   // to override the color at vertex(i,j).
   void operator init(triple c, path3 g, triple axis, int n=nslice,
-                     real angle1=0, real angle2= 360,
+                     real angle1=0, real angle2=360,
                      pen color(int i, real j)=null) {
     axis=unit(axis);
     real w=(angle2-angle1)/n;
     int L=length(g);
     s=new patch[L*n];
+    index=new int[n][L];
     int m=-1;
     transform3[] T=new transform3[n+1];
     transform3 t=rotate(w,c,c+axis);
@@ -728,6 +811,9 @@ struct surface {
     for(int k=1; k <= n; ++k)
       T[k]=T[k-1]*t;
 
+    typedef pen colorfcn(int i, real j);
+    bool defaultcolors=(colorfcn) color == null;
+    
     for(int i=0; i < L; ++i) {
       path3 h=subpath(g,i,i+1);
       path3 r=reverse(h);
@@ -754,13 +840,16 @@ struct surface {
       for(int k=0; k < n; ++k, j += w) {
         transform3 Tp=T[k+1];
         triple dirp=dir(j+w);
-        path3 G=Tk*h{dirj}..{dirp}Tp*r{-dirp}..{-dirj}cycle;
+        path3 G=reverse(Tk*h{dirj}..{dirp}Tp*r{-dirp}..{-dirj}cycle);
         Tk=Tp;
         dirj=dirp;
-        s[++m]=color == null ? patch(G) :
-          patch(G,new pen[] {color(i,j),color(i+1,j),color(i+1,j+w),
-                             color(i,j+w)});
+        s[++m]=defaultcolors ? patch(G) :
+          patch(G,new pen[] {color(i,j),color(i,j+w),color(i+1,j+w),
+                             color(i+1,j)});
+        index[k][i]=m;
       }
+      ucyclic((angle2-angle1) % 360 == 0);
+      vcyclic(cyclic(g));
     }
   }
 
@@ -784,6 +873,9 @@ surface operator * (transform3 t, surface s)
   S.s=new patch[s.s.length];
   for(int i=0; i < s.s.length; ++i)
     S.s[i]=t*s.s[i];
+  S.index=copy(s.index);
+  S.vcyclic=(bool) s.vcyclic;
+  
   return S;
 }
 
@@ -894,10 +986,11 @@ triple[][] subpatchend(triple[][] P, real u, real v)
           {c4[2],c5[2],c6[2],c7[2]}};
 }
 
-patch subpatch(patch s, real ua, real va, real ub, real vb)
+patch subpatch(patch s, pair a, pair b)
 {
-  assert(ua >= 0 && va >= 0 && ub <= 1 && vb <= 1 && ua < ub && va < vb);
-  return patch(subpatchbegin(subpatchend(s.P,ub,vb),ua/ub,va/vb),
+  assert(a.x >= 0 && a.y >= 0 && b.x <= 1 && b.y <= 1 &&
+         a.x < b.x && a.y < b.y);
+  return patch(subpatchbegin(subpatchend(s.P,b.x,b.y),a.x/b.x,a.y/b.y),
                s.straight,s.planar);
 }
 
@@ -958,8 +1051,17 @@ void draw3D(frame f, patch s, material m, light light=currentlight)
   bool lighton=light.on();
   if(!lighton && !invisible((pen) m))
     m=emissive(m);
+  real PRCshininess;
+  if(prc()) {
+    // Empirical translation table from Phong-Blinn to PRC shininess model:
+    static real[] x={0.015,0.025,0.05,0.07,0.1,0.14,0.23,0.5,0.65,0.75,0.85,
+                     0.875,0.9,1};
+    static real[] y={0.05,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.55,0.6,0.7,0.8,0.9,1};
+    static realfunction s=fspline(x,y,monotonic);
+    PRCshininess=s(m.shininess);
+  }
   real granularity=m.granularity >= 0 ? m.granularity : defaultgranularity;
-  draw(f,s.P,s.straight,m.p,m.opacity,m.shininess,granularity,
+  draw(f,s.P,s.straight,m.p,m.opacity,m.shininess,PRCshininess,granularity,
        s.planar ? s.normal(0.5,0.5) : O,lighton,s.colors);
 }
 
@@ -972,7 +1074,7 @@ void tensorshade(transform t=identity(), frame f, patch s,
 }
 
 restricted pen[] nullpens={nullpen};
-nullpens.cyclic(true);
+nullpens.cyclic=true;
 
 void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
           material[] surfacepen, pen[] meshpen=nullpens,
@@ -1011,7 +1113,7 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
 
     depth=sort(depth);
 
-    light.T=shiftless(P.modelview());
+    light.T=shiftless(P.T.modelview);
 
     // Draw from farthest to nearest
     while(depth.length > 0) {
@@ -1033,8 +1135,8 @@ void draw(transform t=identity(), frame f, surface s, int nu=1, int nv=1,
 {
   material[] surfacepen={surfacepen};
   pen[] meshpen={meshpen};
-  surfacepen.cyclic(true);
-  meshpen.cyclic(true);
+  surfacepen.cyclic=true;
+  meshpen.cyclic=true;
   draw(t,f,s,nu,nv,surfacepen,meshpen,light,meshlight,P);
 }
 
@@ -1082,8 +1184,8 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
 {
   material[] surfacepen={surfacepen};
   pen[] meshpen={meshpen};
-  surfacepen.cyclic(true);
-  meshpen.cyclic(true);
+  surfacepen.cyclic=true;
+  meshpen.cyclic=true;
   draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight);
 }
 
@@ -1092,18 +1194,26 @@ void draw(picture pic=currentpicture, surface s, int nu=1, int nv=1,
           light light=currentlight, light meshlight=light)
 {
   pen[] meshpen={meshpen};
-  meshpen.cyclic(true);
+  meshpen.cyclic=true;
   draw(pic,s,nu,nv,surfacepen,meshpen,light,meshlight);
 }
 
-surface extrude(path p, triple axis=Z)
+surface extrude(path3 p, path3 q)
 {
   static patch[] allocate;
-  path3 G=path3(p);
-  path3 G2=shift(axis)*G;
   return surface(...sequence(new patch(int i) {
-        return patch(subpath(G,i,i+1)--subpath(G2,i+1,i)--cycle);
-      },length(G)));
+        return patch(subpath(p,i,i+1)--subpath(q,i+1,i)--cycle);
+      },length(p)));
+}
+
+surface extrude(path3 p, triple axis=Z)
+{
+  return extrude(p,shift(axis)*p);
+}
+
+surface extrude(path p, triple axis=Z)
+{
+  return extrude(path3(p),axis);
 }
 
 surface extrude(explicit path[] p, triple axis=Z)
@@ -1179,11 +1289,11 @@ void label(frame f, Label L, triple position, align align=NoAlign,
   } else {
     if(L.filltype == NoFill)
       fill(f,path(L,project(position,P.t),P),
-           light.color(L.T3*Z,L.p,shiftless(P.modelview())));
+           color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
     else {
       frame d;
       fill(d,path(L,project(position,P.t),P),
-           light.color(L.T3*Z,L.p,shiftless(P.modelview())));
+           color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
       add(f,d,L.filltype);
     }
   }
@@ -1215,11 +1325,11 @@ void label(picture pic=currentpicture, Label L, triple position,
       if(pic != null) {
         if(L.filltype == NoFill)
           fill(project(v,P.t),pic,path(L,P),
-               light.color(L.T3*Z,L.p,shiftless(P.modelview())));
+               color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
         else {
           picture d;
           fill(project(v,P.t),d,path(L,P),
-               light.color(L.T3*Z,L.p,shiftless(P.modelview())));
+               color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
           add(pic,d,L.filltype);
         }
       }
@@ -1271,30 +1381,74 @@ surface extrude(Label L, triple axis=Z)
 
 restricted surface nullsurface;
 
+surface labelsurface(Label L, surface s, real uoffset, real voffset,
+                     real height=0, bool bottom=false, bool top=true)
+{
+  int nu=s.index.length;
+  if(nu == 0) return nullsurface;
+  int nv=s.index[0].length;
+  if(nv == 0) return nullsurface;
+
+  path[] g=texpath(L);
+  pair m=min(g);
+  pair M=max(g);
+  pair lambda=inverse(L.T*scale(nu-epsilon,nv-epsilon))*(M-m);
+  lambda=(abs(lambda.x),abs(lambda.y));
+  path[] G=bezulate(g);
+
+  path3 transpath(path p, real height) {
+    return path3(unstraighten(p),new triple(pair z) {
+        real u=uoffset+(z.x-m.x)/lambda.x;
+        real v=voffset+(z.y-m.y)/lambda.y;
+        if(((u < 0 || u >= nu) && !s.ucyclic()) ||
+           ((v < 0 || v >= nv) && !s.vcyclic()))
+          abort("cannot fit string to surface");
+        return s.point(u,v)+height*unit(s.normal(u,v));
+      });
+  }
+
+  surface s;
+  for(path p : G) {
+    for(path g : regularize(p)) {
+      path3 b;
+      bool extrude=height > 0;
+      if(bottom || extrude) 
+        b=transpath(g,0);
+      if(bottom) s.s.push(patch(b));
+      if(top || extrude) {
+        path3 h=transpath(g,height);
+        if(top) s.s.push(patch(h));
+        if(extrude) s.append(extrude(b,h));
+      }
+    }
+  }
+  return s;
+}
+
 private real a=4/3*(sqrt(2)-1);
-private transform3 t=rotate(90,O,Z);
-private transform3 t2=t*t;
-private transform3 t3=t2*t;
+private transform3 t1=rotate(90,O,Z);
+private transform3 t2=t1*t1;
+private transform3 t3=t2*t1;
 private transform3 i=xscale3(-1)*zscale3(-1);
 
-restricted patch octant1=patch(X{Z}..{-X}Z..Z{Y}..{-Z}Y{X}..{-Y}cycle,
-                               new triple[] {(1,a,a),(a,a^2,1),(a^2,a,1),
-                                             (a,1,a)});
+restricted patch octant1=patch(X{Y}..{-X}Y{Z}..{-Y}Z..Z{X}..{-Z}cycle,
+                               new triple[] {(1,a,a),(a,1,a),(a^2,a,1),
+                                             (a,a^2,1)});
 
-restricted surface unithemisphere=surface(octant1,t*octant1,t2*octant1,
+restricted surface unithemisphere=surface(octant1,t1*octant1,t2*octant1,
                                           t3*octant1);
-restricted surface unitsphere=surface(octant1,t*octant1,t2*octant1,t3*octant1,
-                                      i*octant1,i*t*octant1,i*t2*octant1,
+restricted surface unitsphere=surface(octant1,t1*octant1,t2*octant1,t3*octant1,
+                                      i*octant1,i*t1*octant1,i*t2*octant1,
                                       i*t3*octant1);
 
 restricted patch unitfrustum(real t1, real t2)
 {
   real s1=interp(t1,t2,1/3);
   real s2=interp(t1,t2,2/3);
-  return patch(interp(Z,X,t2)--interp(Z,X,t1){Y}..{-X}interp(Z,Y,t1)--
-               interp(Z,Y,t2){X}..{-Y}cycle,
-               new triple[] {(s2,s2*a,1-s2),(s1,s1*a,1-s1),(s1*a,s1,1-s1),
-                                          (s2*a,s2,1-s2)});
+  return patch(interp(Z,X,t2){Y}..{-X}interp(Z,Y,t2)--interp(Z,Y,t1){X}..{-Y}
+               interp(Z,X,t1)--cycle,
+               new triple[] {(s2,s2*a,1-s2),(s2*a,s2,1-s2),(s1*a,s1,1-s1),
+                                          (s1,s1*a,1-s1)});
 }
 
 // Return a unitcone constructed from n frusta (the final one being degenerate)
@@ -1306,7 +1460,7 @@ surface unitcone(int n=6)
   for(int i=0; i < n; ++i) {
     patch s=unitfrustum(i < n-1 ? r^(i+1) : 0,r^i);
     unitcone.s[i]=s;
-    unitcone.s[n+i]=t*s;
+    unitcone.s[n+i]=t1*s;
     unitcone.s[2n+i]=t2*s;
     unitcone.s[3n+i]=t3*s;
   }
@@ -1316,9 +1470,9 @@ surface unitcone(int n=6)
 restricted surface unitcone=unitcone();
 restricted surface unitsolidcone=surface(patch(unitcircle3)...unitcone.s);
 
-private patch unitcylinder1=patch(X--X+Z{Y}..{-X}Y+Z--Y{X}..{-Y}cycle);
+private patch unitcylinder1=patch(X{Y}..{-X}Y--Y+Z{X}..{-Y}X+Z--cycle);
 
-restricted surface unitcylinder=surface(unitcylinder1,t*unitcylinder1,
+restricted surface unitcylinder=surface(unitcylinder1,t1*unitcylinder1,
                                         t2*unitcylinder1,t3*unitcylinder1);
 
 private patch unitplane=patch(new triple[] {O,X,X+Y,Y});
diff --git a/Master/texmf/asymptote/three_tube.asy b/Master/texmf/asymptote/three_tube.asy
new file mode 100644
index 00000000000..c79ae3f619f
--- /dev/null
+++ b/Master/texmf/asymptote/three_tube.asy
@@ -0,0 +1,373 @@
+void render(path3 s, real granularity=linegranularity, void f(path3, real))
+{
+  void Split(triple z0, triple c0, triple c1, triple z1, real t0=0, real t1=1,
+             real depth=mantissaBits) {
+    if(depth > 0) {
+      real S=straightness(z0,c0,c1,z1);
+      if(S > 0) {
+        --depth;
+        if(S > max(granularity*max(abs(z0),abs(c0),abs(c1),abs(z1)))) {
+          triple m0=0.5*(z0+c0);
+          triple m1=0.5*(c0+c1);
+          triple m2=0.5*(c1+z1);
+          triple m3=0.5*(m0+m1);
+          triple m4=0.5*(m1+m2);
+          triple m5=0.5*(m3+m4);
+          real tm=0.5*(t0+t1);
+          Split(z0,m0,m3,m5,t0,tm,depth);
+          Split(m5,m4,m2,z1,tm,t1,depth);
+          return;
+        }
+      }
+    }
+    f(z0..controls c0 and c1..z1,t0);
+  }
+  Split(point(s,0),postcontrol(s,0),precontrol(s,1),point(s,1));
+}
+
+struct rmf
+{
+  triple p,r,t,s;
+  void operator init(triple p, triple r, triple t)
+  {
+    this.p=p;
+    this.r=r;
+    this.t=t;
+    s=cross(t,r);
+  }
+}
+
+// Rotation minimizing frame
+// http://www.cs.hku.hk/research/techreps/document/TR-2007-07.pdf
+rmf[] rmf(path3 g, real[] t)
+{
+  rmf[] R=new rmf[t.length];
+  triple d=dir(g,0);
+  R[0]=rmf(point(g,0),perp(d),d);
+  for(int i=1; i < t.length; ++i) {
+    rmf Ri=R[i-1];
+    real t=t[i];
+    triple p=point(g,t);
+    triple v1=p-Ri.p;
+    if(v1 != O) {
+      triple r=Ri.r;
+      triple u1=unit(v1);
+      triple ti=Ri.t;
+      triple tp=ti-2*dot(u1,ti)*u1;
+      ti=dir(g,t);
+      triple rp=r-2*dot(u1,r)*u1;
+      triple u2=unit(ti-tp);
+      rp=rp-2*dot(u2,rp)*u2;
+      R[i]=rmf(p,unit(rp),unit(ti));
+    } else
+      R[i]=R[i-1];
+  }
+  return R;
+}
+
+surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
+                 real[] x, real[] y, bool[][] cond={})
+{ // z[i][j] is the value at (x[i],y[j])
+  // p and q are the first derivatives with respect to x and y, respectively
+  // r is the second derivative ddu/dxdy
+  int n=x.length-1;
+  int m=y.length-1;
+
+  bool all=cond.length == 0;
+
+  int count;
+  if(all)
+    count=n*m;
+  else {
+    count=0;
+    for(int i=0; i < n; ++i) {
+      bool[] condi=cond[i];
+      for(int j=0; j < m; ++j)
+        if(condi[j]) ++count;
+    }
+  }
+
+  surface s=surface(count);
+  s.index=new int[n][m];
+  int k=-1;
+  for(int i=0; i < n; ++i) {
+    bool[] condi=all ? null : cond[i];
+    real xi=x[i];
+    real[] zi=z[i];
+    real[] zp=z[i+1];
+    real[] ri=r[i];
+    real[] rp=r[i+1];
+    real[] pi=p[i];
+    real[] pp=p[i+1];
+    real[] qi=q[i];
+    real[] qp=q[i+1];
+    real xp=x[i+1];
+    real hx=(xp-xi)/3;
+    int[] indexi=s.index[i];
+    for(int j=0; j < m; ++j) {
+      real yj=y[j];
+      real yp=y[j+1];
+      if(all || condi[j]) {
+        triple[][] P=array(4,array(4,O));
+        real hy=(yp-yj)/3;
+        real hxy=hx*hy;
+        // x and y directions
+        for(int k=0; k < 4; ++k) {
+          P[0][k] += xi*X;
+          P[k][0] += yj*Y;
+          P[1][k] += (xp+2*xi)/3*X;
+          P[k][1] += (yp+2*yj)/3*Y;
+          P[2][k] += (2*xp+xi)/3*X;
+          P[k][2] += (2*yp+yj)/3*Y;
+          P[3][k] += xp*X;
+          P[k][3] += yp*Y;
+        }
+        // z: value 
+        P[0][0] += zi[j]*Z;
+        P[3][0] += zp[j]*Z;
+        P[0][3] += zi[j+1]*Z;
+        P[3][3] += zp[j+1]*Z;
+        // z: first derivative
+        P[1][0] += (P[0][0].z+hx*pi[j])*Z;
+        P[1][3] += (P[0][3].z+hx*pi[j+1])*Z;
+        P[2][0] += (P[3][0].z-hx*pp[j])*Z;
+        P[2][3] += (P[3][3].z-hx*pp[j+1])*Z;
+        P[0][1] += (P[0][0].z+hy*qi[j])*Z;
+        P[3][1] += (P[3][0].z+hy*qp[j])*Z;
+        P[0][2] += (P[0][3].z-hy*qi[j+1])*Z;
+        P[3][2] += (P[3][3].z-hy*qp[j+1])*Z;
+        // z: second derivative
+        P[1][1] += (P[0][1].z+P[1][0].z-P[0][0].z+hxy*ri[j])*Z;
+        P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*ri[j+1])*Z;
+        P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*rp[j])*Z;
+        P[2][2] += (P[2][3].z+P[3][2].z-P[3][3].z+hxy*rp[j+1])*Z;
+        s.s[++k]=patch(P);
+        indexi[j]=k;
+      }
+    }
+  }
+  
+  return s;
+}
+
+// return the surface described by a real matrix f, interpolated with
+// xsplinetype and ysplinetype.
+surface surface(real[][] f, real[] x, real[] y,
+                splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype,
+                bool[][] cond={})
+{
+  real epsilon=sqrtEpsilon*norm(y);
+  if(xsplinetype == null)
+    xsplinetype=(abs(x[0]-x[x.length-1]) <= epsilon) ? periodic : notaknot;
+  if(ysplinetype == null)
+    ysplinetype=(abs(y[0]-y[y.length-1]) <= epsilon) ? periodic : notaknot;
+  int n=x.length; int m=y.length;
+  real[][] ft=transpose(f);
+  real[][] tp=new real[m][];
+  for(int j=0; j < m; ++j)
+    tp[j]=xsplinetype(x,ft[j]);
+  real[][] q=new real[n][];
+  for(int i=0; i < n; ++i)
+    q[i]=ysplinetype(y,f[i]);
+  real[][] qt=transpose(q);
+  real[] d1=xsplinetype(x,qt[0]);
+  real[] d2=xsplinetype(x,qt[m-1]);
+  real[][] r=new real[n][];
+  real[][] p=transpose(tp);
+  for(int i=0; i < n; ++i)
+    r[i]=clamped(d1[i],d2[i])(y,p[i]);
+  surface s=bispline(f,p,q,r,x,y,cond);
+  if(xsplinetype == periodic) s.ucyclic(true);
+  if(ysplinetype == periodic) s.vcyclic(true);
+  return s;
+}
+
+bool uperiodic(real[][] a) {
+  int n=a.length;
+  if(n == 0) return false;
+  int m=a[0].length;
+  real[] a0=a[0];
+  real[] a1=a[n-1];
+  real epsilon=sqrtEpsilon*norm(a);
+  for(int j=0; j < m; ++j)
+    if(abs(a0[j]-a1[j]) > epsilon) return false;
+  return true;
+}
+bool vperiodic(real[][] a) {
+  int n=a.length;
+  if(n == 0) return false;
+  int m=a[0].length-1;
+  real epsilon=sqrtEpsilon*norm(a);
+  for(int i=0; i < n; ++i)
+    if(abs(a[i][0]-a[i][m]) > epsilon) return false;
+  return true;
+}
+
+// return the surface described by a parametric function f evaluated at u and v
+// and interpolated with usplinetype and vsplinetype.
+surface surface(triple f(pair z), real[] u, real[] v,
+                splinetype[] usplinetype, splinetype[] vsplinetype=Spline,
+                bool cond(pair z)=null)
+{
+  int nu=u.length-1;
+  int nv=v.length-1;
+  real[] ipt=sequence(u.length);
+  real[] jpt=sequence(v.length);
+  real[][] fx=new real[u.length][v.length];
+  real[][] fy=new real[u.length][v.length];
+  real[][] fz=new real[u.length][v.length];
+
+  bool[][] active;
+  bool all=cond == null;
+  if(!all) active=new bool[u.length][v.length];
+
+  for(int i=0; i <= nu; ++i) {
+    real ui=u[i];
+    real[] fxi=fx[i];
+    real[] fyi=fy[i];
+    real[] fzi=fz[i];
+    bool[] activei=all ? null : active[i];
+    for(int j=0; j <= nv; ++j) {
+      pair z=(ui,v[j]);
+      triple f=(all || (activei[j]=cond(z))) ? f(z) : O;
+      fxi[j]=f.x;
+      fyi[j]=f.y;
+      fzi[j]=f.z;
+    }
+  }
+
+  if(usplinetype.length == 0) {
+    usplinetype=new splinetype[] {uperiodic(fx) ? periodic : notaknot,
+                                  uperiodic(fy) ? periodic : notaknot,
+                                  uperiodic(fz) ? periodic : notaknot};
+  } else if(usplinetype.length != 3) abort("usplinetype must have length 3");
+
+  if(vsplinetype.length == 0) {
+    vsplinetype=new splinetype[] {vperiodic(fx) ? periodic : notaknot,
+                                  vperiodic(fy) ? periodic : notaknot,
+                                  vperiodic(fz) ? periodic : notaknot};
+  } else if(vsplinetype.length != 3) abort("vsplinetype must have length 3");
+  
+  surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
+  surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
+  surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
+
+  surface s=surface(sx.s.length);
+  s.index=new int[nu][nv];
+  int k=-1;
+  for(int i=0; i < nu; ++i) {
+    int[] indexi=s.index[i];
+    for(int j=0; j < nv; ++j)
+      indexi[j]=++k;
+  }
+
+  for(int k=0; k < sx.s.length; ++k) {
+    triple[][] Q=new triple[4][];
+    for(int i=0; i < 4 ; ++i)
+      Q[i]=sequence(new triple(int j) {
+          return (sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z);
+        },4);
+    s.s[k]=patch(Q);
+  }
+
+  if(usplinetype[0] == periodic && usplinetype[1] == periodic &&
+     usplinetype[1] == periodic) s.ucyclic(true);
+
+  if(vsplinetype[0] == periodic && vsplinetype[1] == periodic &&
+     vsplinetype[1] == periodic) s.vcyclic(true);
+  
+  return s;
+}
+
+path3 interp(path3 a, path3 b, real t) 
+{
+  int n=size(a);
+  return path3(sequence(new triple(int i) {return interp(precontrol(a,i),
+                                                         precontrol(b,i),t);},n),
+    sequence(new triple(int i) {return interp(point(a,i),point(b,i),t);},n),
+    sequence(new triple(int i) {return interp(postcontrol(a,i),
+                                              postcontrol(b,i),t);},n),
+    sequence(new bool(int i) {return straight(a,i) && straight(b,i);},n),
+    cyclic(a) && cyclic(b));
+}
+
+struct tube
+{
+  surface s;
+  path3 center;
+
+  void operator init(path3 p, real width, int sectors=4,
+                     real granularity=linegranularity) {
+    sectors += sectors % 2; // Must be even.
+    int h=quotient(sectors,2);
+    real r=0.5*width;
+
+    void generate(path3 p) {
+      int n=length(p);
+      if(piecewisestraight(p)) {
+        for(int i=0; i < n; ++i) {
+          triple v=point(p,i);
+          triple u=point(p,i+1)-v;
+          s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder);
+        }
+        center=center&p;
+      } else {
+        real[] T;
+        path3 G;
+        for(int i=0; i < n; ++i)
+          render(subpath(p,i,i+1),granularity,
+                 new void(path3 g, real s) {
+                   G=G&g;
+                   T.push(i+s);
+                 });
+        T.push(n);
+        T.cyclic=cyclic(p);
+        rmf[] rmf=rmf(p,T);
+        triple f(pair t) {
+          rmf R=rmf[round(t.x)];
+          return point(G,t.x)+r*(R.r*cos(t.y)-R.s*sin(t.y));
+        }
+
+        real[] v=uniform(0,2pi,sectors);
+        static splinetype[] Monotonic={monotonic,monotonic,monotonic};
+        static splinetype[] Periodic={periodic,periodic,periodic};
+        if(T.length > 0) {
+          surface S=surface(f,sequence(T.length),v,Monotonic,Periodic);
+          s.append(S);
+
+          // Compute center of tube:
+          int n=S.index.length;
+          if(T.cyclic) --n;
+          triple[] pre=new triple[n+1];
+          triple[] point=new triple[n+1];
+          triple[] post=new triple[n+1];
+          int[] index=S.index[0];
+          pre[0]=point[0]=0.5*(S.s[index[0]].P[0][0]+S.s[index[h]].P[0][0]);
+          for(int i=0; i < n; ++i) {
+            index=S.index[i];
+            triple [][] P=S.s[index[0]].P;
+            triple [][] Q=S.s[index[h]].P;
+            post[i]=0.5*(P[1][0]+Q[1][0]);
+            pre[i+1]=0.5*(P[2][0]+Q[2][0]);
+            point[i+1]=0.5*(P[3][0]+Q[3][0]);
+          }
+          index=S.index[n-1];
+          post[n]=0.5*(S.s[index[0]].P[3][0]+S.s[index[h]].P[3][0]);
+          center=center&path3(pre,point,post,array(n+1,false),T.cyclic);
+        }
+      }
+    }
+    
+    transform3 t=scale3(r);
+    bool cyclic=cyclic(p);
+    int begin=0;
+    int n=length(p);
+    for(int i=cyclic ? 0 : 1; i < n; ++i)
+     if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) {
+       generate(subpath(p,begin,i));
+       s.append(shift(point(p,i))*t*align(dir(p,i,-1))*unithemisphere);
+       begin=i;
+     }
+    generate(subpath(p,begin,n));
+  }
+}
diff --git a/Master/texmf/asymptote/tree.asy b/Master/texmf/asymptote/tree.asy
index 0f3559112c3..1e603ec4b6b 100644
--- a/Master/texmf/asymptote/tree.asy
+++ b/Master/texmf/asymptote/tree.asy
@@ -7,8 +7,8 @@
 
 struct tree
 {
-  tree left = null;
-  tree right = null;
+  tree left;
+  tree right;
   int key = 0;
   int value = 0;
 }
diff --git a/Master/texmf/asymptote/trembling.asy b/Master/texmf/asymptote/trembling.asy
index fc50ce9b7f5..b01d572d6e2 100644
--- a/Master/texmf/asymptote/trembling.asy
+++ b/Master/texmf/asymptote/trembling.asy
@@ -84,7 +84,7 @@ real atime(pair m, path g, real fuzz=trembleFuzz())
   //   if(t.length <= 0) { // degenerated case
   //     r=initr;
   //     T *= shift(shx,0);
-  //     write("warning: atime needs numerical adjustment.");
+  //     warning("atime","atime needs numerical adjustment.",position=true);
   //   }
   // } while(t.length <= 0);
   if(t.length > 0) ot=t[1];
diff --git a/Master/texmf/asymptote/tube.asy b/Master/texmf/asymptote/tube.asy
index 1f8ef950251..844a39c92f0 100644
--- a/Master/texmf/asymptote/tube.asy
+++ b/Master/texmf/asymptote/tube.asy
@@ -58,20 +58,7 @@ real[] sample(path3 g, real r, real relstep=0)
   return t;
 }
 
-struct Rmf
-{
-  triple p,r,t; // s=cross(t,r);
-  real reltime;
-  void operator init(triple p, triple r, triple t, real reltime)
-  {
-    this.p=p;
-    this.r=r;
-    this.t=t;
-    this.reltime=reltime;
-  }
-}
-
-real degrees(Rmf a, Rmf b)
+real degrees(rmf a, rmf b)
 {
   real d=degrees(acos1(dot(a.r,b.r)));
   real dt=dot(cross(a.r,b.r),a.t);
@@ -79,37 +66,6 @@ real degrees(Rmf a, Rmf b)
   return d%360;
 }
 
-private Rmf[] rmf(path3 g, Rmf U0=Rmf(O,O,O,0), real[] t)
-{
-  if(U0.t == O) {
-    triple d=dir(g,0);
-    U0=Rmf(point(g,0),perp(d),d,0);
-  }
-  int l=length(g);
-  Rmf[] R={U0};
-  triple rp,v1,v2,tp,ti,p;
-  real c;
-
-  for(int i=0; i < t.length-1; ++i) {
-    p=point(g,t[i+1]);
-    v1=p-R[i].p;
-    c=dot(v1,v1);
-    if(c != 0) {
-      rp=R[i].r-2*dot(v1,R[i].r)*v1/c;
-      ti=R[i].t;
-      tp=ti-2*dot(v1,ti)*v1/c;
-      ti=dir(g,t[i+1]);
-      v2=ti-tp;
-      rp=rp-2*dot(v2,rp)*v2/dot(v2,v2);
-      R.push(Rmf(p,unit(rp),unit(ti),t[i+1]/l));
-    } else {
-      write("Warning: path3 has duplicated point in Rmf.");
-      R.push(R[R.length-1]);
-    }
-  }
-  return R;
-}
-
 restricted int coloredNodes=1;
 restricted int coloredSegments=2;
 
@@ -158,8 +114,7 @@ coloredpath operator cast(guide p)
   return coloredpath(p);
 }
 
-private surface surface(Rmf[] R, coloredpath cp,transform T(real)=
-			new transform(real t) {return identity();},
+private surface surface(rmf[] R, real[] t, coloredpath cp, transform T(real),
 			bool cyclic)
 {
   path g=cp.p;
@@ -169,22 +124,22 @@ private surface surface(Rmf[] R, coloredpath cp,transform T(real)=
     planar[i]=straight(g,i);
 
   surface s;
-  path3 sec=path3(T(R[0].reltime)*g);
+  path3 sec=path3(T(t[0]/l)*g);
   real adjust=0;
   if(cyclic) adjust=-degrees(R[0],R[R.length-1])/(R.length-1);
-  path3 sec1=shift(R[0].p)*transform3(R[0].r,cross(R[0].t,R[0].r),R[0].t)*sec,
+  path3 sec1=shift(R[0].p)*transform3(R[0].r,R[0].s,R[0].t)*sec,
     sec2;
 
   for(int i=1; i < R.length; ++i) {
-    sec=path3(T(R[i].reltime)*g);
+    sec=path3(T(t[i]/l)*g);
     sec2=shift(R[i].p)*transform3(R[i].r,cross(R[i].t,R[i].r),R[i].t)*
       rotate(i*adjust,Z)*sec;
     for(int j=0; j < l; ++j) {
       surface st=surface(subpath(sec1,j,j+1)--subpath(sec2,j+1,j)--cycle,
 			 planar=planar[j]);
       if(cp.usepens) {
-        pen[] tp1=cp.pens(R[i-1].reltime), tp2=cp.pens(R[i].reltime);
-        tp1.cyclic(true); tp2.cyclic(true);
+        pen[] tp1=cp.pens(t[i-1]/l), tp2=cp.pens(t[i]/l);
+        tp1.cyclic=true; tp2.cyclic=true;
         if(cp.colortype == coloredSegments) {
           st.colors(new pen[][] {{tp1[j],tp1[j],tp2[j],tp2[j]}});
 	} else {
@@ -206,6 +161,6 @@ surface tube(path3 g, coloredpath section,
   real[] t=sample(g,max(M.x-m.x,M.y-m.y)/max(realEpsilon,abs(corner)),
                   min(abs(relstep),1));
   bool cyclic=cyclic(g);
-  t.cyclic(cyclic);
-  return surface(rmf(g,t),section,T,cyclic);
+  t.cyclic=cyclic;
+  return surface(rmf(g,t),t,section,T,cyclic);
 }
diff --git a/Master/texmf/asymptote/version.asy b/Master/texmf/asymptote/version.asy
index 713dc924f1c..d6c18f1b6de 100644
--- a/Master/texmf/asymptote/version.asy
+++ b/Master/texmf/asymptote/version.asy
@@ -1 +1 @@
-string VERSION="1.82";
+string VERSION="1.83";