asymptote 1.77 sources

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

11 files changed, 1231 insertions, 93 deletions

diff --git a/Build/source/utils/asymptote/base/graph.asy b/Build/source/utils/asymptote/base/graph.asy
index 805cf15f612..2983137b612 100644
--- a/Build/source/utils/asymptote/base/graph.asy
+++ b/Build/source/utils/asymptote/base/graph.asy
@@ -518,16 +518,17 @@ string autoformat(string format="", real norm ... real[] a)
 
   bool signchange=(A.length > 1 && A[0] < 0 && A[A.length-1] >= 0);
 
-  for(int i=0; i < A.length-1; ++i)
+  for(int i=0; i < A.length; ++i)
     if(a[i] < zerotickfuzz*norm) A[i]=a[i]=0;
 
   int n=0;
 
   bool Fixed=find(a >= 1e4-epsilon | (a > 0 & a <= 1e-4-epsilon)) < 0;
+  
   string Format=defaultformat(4,fixed=Fixed);
 
   if(Fixed && n < 4) {
-    for(int i=0; i < A.length-1; ++i) {
+    for(int i=0; i < A.length; ++i) {
       real a=A[i];
       while(format(defaultformat(n,fixed=Fixed),a) != format(Format,a))
         ++n;
diff --git a/Build/source/utils/asymptote/base/graph3.asy b/Build/source/utils/asymptote/base/graph3.asy
index a1c8258638f..8fe6c9ed5c6 100644
--- a/Build/source/utils/asymptote/base/graph3.asy
+++ b/Build/source/utils/asymptote/base/graph3.asy
@@ -1591,15 +1591,11 @@ private surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
       real yj=y[j];
       real yp=y[j+1];
       if(all || condi[j]) {
-        triple[][] P={
-          {O,O,O,O},
-          {O,O,O,O},
-          {O,O,O,O},
-          {O,O,O,O}};
+        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) {
+        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;
@@ -1636,34 +1632,38 @@ private surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
 }
 
 // return the surface described by a real matrix f, interpolated with
-// splinetype.
+// xsplinetype and ysplinetype.
 surface surface(real[][] f, real[] x, real[] y,
-                splinetype splinetype=null, bool[][] cond={})
+                splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype,
+                bool[][] cond={})
 {
-  if(splinetype == null)
-    splinetype=(x[0] == x[x.length-1] && y[0] == y[y.length-1]) ? 
-      periodic : notaknot;
+  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]=splinetype(x,ft[j]);
+  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]=splinetype(y,f[i]);
+  for(int i=0; i < n; ++i)
+    q[i]=ysplinetype(y,f[i]);
   real[][] qt=transpose(q);
-  real[] d1=splinetype(x,qt[0]);
-  real[] d2=splinetype(x,qt[m-1]);
+  real[] d1=xsplinetype(x,qt[0]);
+  real[] d2=xsplinetype(x,qt[m-1]);
   real[][] r=new real[n][];
-  for(int i=0; i < n ; ++i)
-    r[i]=clamped(d1[i],d2[i])(y,f[i]);
-  return bispline(f,transpose(tp),q,r,x,y,cond);
+  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);
 }
 
 // return the surface described by a real matrix f, interpolated with
-// splinetype.
-surface surface(real[][] f, pair a, pair b, splinetype splinetype,
-                bool[][] cond={})
+// xsplinetype and ysplinetype.
+surface surface(real[][] f, pair a, pair b, splinetype xsplinetype,
+                splinetype ysplinetype=xsplinetype, bool[][] cond={})
 {
   if(!rectangular(f)) abort("matrix is not rectangular");
 
@@ -1674,7 +1674,7 @@ surface surface(real[][] f, pair a, pair b, splinetype splinetype,
 
   real[] x=uniform(a.x,b.x,nx);
   real[] y=uniform(a.y,b.y,ny);
-  return surface(f,x,y,splinetype,cond);
+  return surface(f,x,y,xsplinetype,ysplinetype,cond);
 }
 
 // return the surface described by a real matrix f, interpolated linearly.
@@ -1732,6 +1732,82 @@ surface surface(triple f(pair z), pair a, pair b, int nu=nmesh, int nv=nu,
   return surface(v,active);
 }
   
+// return the surface described by a parametric function f over box(a,b),
+// interpolated with usplinetype and vsplinetype.
+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 the surface described by a real function f over box(a,b),
 // interpolated linearly.
 surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
@@ -1741,9 +1817,10 @@ surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
 }
 
 // return the surface described by a real function f over box(a,b),
-// interpolated with splinetype.
+// interpolated with xsplinetype and ysplinetype.
 surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
-                splinetype splinetype, bool cond(pair z)=null)
+                splinetype xsplinetype, splinetype ysplinetype=xsplinetype,
+                bool cond(pair z)=null)
 {
   bool[][] active;
   bool all=cond == null;
@@ -1767,7 +1844,7 @@ surface surface(real f(pair z), pair a, pair b, int nx=nmesh, int ny=nx,
       if(!all) activei[j]=cond(z);
     }
   }
-  return surface(F,x,y,splinetype,active);
+  return surface(F,x,y,xsplinetype,ysplinetype,active);
 }
 
 guide3[][] lift(real f(real x, real y), guide[][] g,
diff --git a/Build/source/utils/asymptote/base/lmfit.asy b/Build/source/utils/asymptote/base/lmfit.asy
new file mode 100644
index 00000000000..d0367547b74
--- /dev/null
+++ b/Build/source/utils/asymptote/base/lmfit.asy
@@ -0,0 +1,838 @@
+/*
+  Copyright (c) 2009 Philipp Stephani
+
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation files
+  (the "Software"), to deal in the Software without restriction,
+  including without limitation the rights to use, copy, modify, merge,
+  publish, distribute, sublicense, and/or sell copies of the Software,
+  and to permit persons to whom the Software is furnished to do so,
+  subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT.  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+  SOFTWARE.
+*/
+
+/*
+  This module provides an implementation of the Levenberg--Marquardt
+  (LM) algorithm, converted from the C lmfit routine by Joachim Wuttke
+  (see http://www.messen-und-deuten.de/lmfit/).
+
+  Implementation strategy: Fortunately, Asymptote's syntax is very
+  similar to C, and the original code cleanly separates the
+  customizable parts (user-provided data, output routines, etc.) from
+  the dirty number crunching.  Thus, mst of the code was just copied
+  and slightly modified from the original source files.  I have
+  amended the lm_data_type structure and the callback routines with a
+  weight array that can be used to provide experimental errors.  I
+  have also created two simple wrapper functions.
+*/
+
+
+// copied from the C code
+private real LM_MACHEP = realEpsilon;
+private real LM_DWARF = realMin;
+private real LM_SQRT_DWARF = sqrt(realMin);
+private real LM_SQRT_GIANT = sqrt(realMax);
+private real LM_USERTOL = 30 * LM_MACHEP;
+
+restricted string lm_infmsg[] = {
+  "improper input parameters",
+  "the relative error in the sum of squares is at most tol",
+  "the relative error between x and the solution is at most tol",
+  "both errors are at most tol",
+  "fvec is orthogonal to the columns of the jacobian to machine precision",
+  "number of calls to fcn has reached or exceeded maxcall*(n+1)",
+  "ftol is too small: no further reduction in the sum of squares is possible",
+  "xtol too small: no further improvement in approximate solution x possible",
+  "gtol too small: no further improvement in approximate solution x possible",
+  "not enough memory",
+  "break requested within function evaluation"
+};
+
+restricted string lm_shortmsg[] = {
+  "invalid input",
+  "success (f)",
+  "success (p)",
+  "success (f,p)",
+  "degenerate",
+  "call limit",
+  "failed (f)",
+  "failed (p)",
+  "failed (o)",
+  "no memory",
+  "user break"
+};
+
+
+// copied from the C code and amended with the weight (user_w) array
+struct lm_data_type {
+  real[] user_t;
+  real[] user_y;
+  real[] user_w;
+  real user_func(real user_t_point, real[] par);  
+};
+
+
+// Asymptote has no pointer support, so we need reference wrappers for
+// the int and real types
+struct lm_int_type {
+  int val;
+  
+  void operator init(int val) {
+    this.val = val;
+  }
+};
+
+
+struct lm_real_type {
+  real val;
+  
+  void operator init(real val) {
+    this.val = val;
+  }
+};
+
+
+// copied from the C code; the lm_initialize_control function turned
+// into a constructor
+struct lm_control_type {
+  real ftol;
+  real xtol;
+  real gtol;
+  real epsilon;
+  real stepbound;
+  real fnorm;
+  int maxcall;
+  lm_int_type nfev;
+  lm_int_type info;
+
+  void operator init() {
+    maxcall = 100;
+    epsilon = LM_USERTOL;
+    stepbound = 100;
+    ftol = LM_USERTOL;
+    xtol = LM_USERTOL;
+    gtol = LM_USERTOL;
+  }
+};
+
+
+// copied from the C code
+typedef void lm_evaluate_ftype(real[] par, int m_dat, real[] fvec, lm_data_type data, lm_int_type info);
+typedef void lm_print_ftype(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev);
+
+
+// copied from the C code
+private real SQR(real x) {
+  return x * x;
+}
+
+
+// Asymptote doesn't support pointers to arbitrary array elements, so
+// we provide an offset parameter.
+private real lm_enorm(int n, real[] x, int offset=0) {
+  real s1 = 0;
+  real s2 = 0;
+  real s3 = 0;
+  real x1max = 0;
+  real x3max = 0;
+  real agiant = LM_SQRT_GIANT / n;
+  real xabs, temp;
+
+  for (int i = 0; i < n; ++i) {
+    xabs = fabs(x[offset + i]);
+    if (xabs > LM_SQRT_DWARF && xabs < agiant) {
+      s2 += SQR(xabs);
+      continue;
+    }
+
+    if (xabs > LM_SQRT_DWARF) {
+      if (xabs > x1max) {
+        temp = x1max / xabs;
+        s1 = 1 + s1 * SQR(temp);
+        x1max = xabs;
+      } else {
+        temp = xabs / x1max;
+        s1 += SQR(temp);
+      }
+      continue;
+    }
+    if (xabs > x3max) {
+      temp = x3max / xabs;
+      s3 = 1 + s3 * SQR(temp);
+      x3max = xabs;
+    } else {
+      if (xabs != 0.0) {
+        temp = xabs / x3max;
+        s3 += SQR(temp);
+      }
+    }
+  }
+
+  if (s1 != 0)
+    return x1max * sqrt(s1 + (s2 / x1max) / x1max);
+  if (s2 != 0) {
+    if (s2 >= x3max)
+      return sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
+    else
+      return sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+  }
+
+  return x3max * sqrt(s3);
+}
+
+
+// This function calculated the vector whose square sum is to be
+// minimized.  We use a slight modification of the original code that
+// includes the weight factor.  The user may provide different
+// customizations.
+void lm_evaluate_default(real[] par, int m_dat, real[] fvec, lm_data_type data, lm_int_type info) {
+  for (int i = 0; i < m_dat; ++i) {
+    fvec[i] = data.user_w[i] * (data.user_y[i] - data.user_func(data.user_t[i], par));
+  }
+}
+
+
+// Helper functions to print padded strings and numbers (until
+// Asymptote provides a real printf function)
+private string pad(string str, int count, string pad=" ") {
+  string res = str;
+  while (length(res) < count)
+    res = pad + res;
+  return res;
+}
+
+
+private string pad(int num, int digits, string pad=" ") {
+  return pad(string(num), digits, pad);
+}
+
+
+private string pad(real num, int digits, string pad=" ") {
+  return pad(string(num), digits, pad);
+}
+
+
+// Similar to the C code, also prints weights
+void lm_print_default(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev) {
+  real f, y, t, w;
+  int i;
+
+  if (iflag == 2) {
+    write("trying step in gradient direction");
+  } else if (iflag == 1) {
+    write(format("determining gradient (iteration %d)", iter));
+  } else if (iflag == 0) {
+    write("starting minimization");
+  } else if (iflag == -1) {
+    write(format("terminated after %d evaluations", nfev));
+  }
+
+  write("  par: ", none);
+  for (i = 0; i < n_par; ++i) {
+    write(" " + pad(par[i], 12), none);
+  }
+  write(" => norm: " + pad(lm_enorm(m_dat, fvec), 12));
+
+  if (iflag == -1) {
+    write("  fitting data as follows:");
+    for (i = 0; i < m_dat; ++i) {
+      t = data.user_t[i];
+      y = data.user_y[i];
+      w = data.user_w[i];
+      f = data.user_func(t, par);
+      write(format("    t[%2d]=", i) + pad(t, 12) + " y=" + pad(y, 12) + " w=" + pad(w, 12) + " fit=" + pad(f, 12) + " residue=" + pad(y - f, 12));
+    }
+  }
+}
+
+
+// Prints nothing
+void lm_print_quiet(int n_par, real[] par, int m_dat, real[] fvec, lm_data_type data, int iflag, int iter, int nfev) {
+}
+
+
+// copied from the C code
+private void lm_qrfac(int m, int n, real[] a, bool pivot, int[] ipvt, real[] rdiag, real[] acnorm, real[] wa) {
+  int i, j, k, kmax, minmn;
+  real ajnorm, sum, temp;
+  static real p05 = 0.05;
+
+  for (j = 0; j < n; ++j) {
+    acnorm[j] = lm_enorm(m, a, j * m);
+    rdiag[j] = acnorm[j];
+    wa[j] = rdiag[j];
+    if (pivot)
+      ipvt[j] = j;
+  }
+
+  minmn = min(m, n);
+  for (j = 0; j < minmn; ++j) {
+    while (pivot) {
+      kmax = j;
+      for (k = j + 1; k < n; ++k)
+        if (rdiag[k] > rdiag[kmax])
+          kmax = k;
+      if (kmax == j)
+        break;
+
+      for (i = 0; i < m; ++i) {
+        temp = a[j * m + i];
+        a[j * m + i] = a[kmax * m + i];
+        a[kmax * m + i] = temp;
+      }
+      rdiag[kmax] = rdiag[j];
+      wa[kmax] = wa[j];
+      k = ipvt[j];
+      ipvt[j] = ipvt[kmax];
+      ipvt[kmax] = k;
+
+      break;
+    }
+
+    ajnorm = lm_enorm(m - j, a, j * m + j);
+    if (ajnorm == 0.0) {
+      rdiag[j] = 0;
+      continue;
+    }
+
+    if (a[j * m + j] < 0.0)
+      ajnorm = -ajnorm;
+    for (i = j; i < m; ++i)
+      a[j * m + i] /= ajnorm;
+    a[j * m + j] += 1;
+
+    for (k = j + 1; k < n; ++k) {
+      sum = 0;
+
+      for (i = j; i < m; ++i)
+        sum += a[j * m + i] * a[k * m + i];
+
+      temp = sum / a[j + m * j];
+
+      for (i = j; i < m; ++i)
+        a[k * m + i] -= temp * a[j * m + i];
+
+      if (pivot && rdiag[k] != 0.0) {
+        temp = a[m * k + j] / rdiag[k];
+        temp = max(0.0, 1 - SQR(temp));
+        rdiag[k] *= sqrt(temp);
+        temp = rdiag[k] / wa[k];
+        if (p05 * SQR(temp) <= LM_MACHEP) {
+          rdiag[k] = lm_enorm(m - j - 1, a, m * k + j + 1);
+          wa[k] = rdiag[k];
+        }
+      }
+    }
+
+    rdiag[j] = -ajnorm;
+  }
+}
+
+
+// copied from the C code
+private void lm_qrsolv(int n, real[] r, int ldr, int[] ipvt, real[] diag, real[] qtb, real[] x, real[] sdiag, real[] wa) {
+  static real p25 = 0.25;
+  static real p5 = 0.5;
+
+  int i, kk, j, k, nsing;
+  real qtbpj, sum, temp;
+  real _sin, _cos, _tan, _cot;
+
+  for (j = 0; j < n; ++j) {
+    for (i = j; i < n; ++i)
+      r[j * ldr + i] = r[i * ldr + j];
+    x[j] = r[j * ldr + j];
+    wa[j] = qtb[j];
+  }
+
+  for (j = 0; j < n; ++j) {
+    while (diag[ipvt[j]] != 0.0) {
+      for (k = j; k < n; ++k)
+        sdiag[k] = 0.0;
+      sdiag[j] = diag[ipvt[j]];
+
+      qtbpj = 0.;
+      for (k = j; k < n; ++k) {
+        if (sdiag[k] == 0.)
+          continue;
+        kk = k + ldr * k;
+        if (fabs(r[kk]) < fabs(sdiag[k])) {
+          _cot = r[kk] / sdiag[k];
+          _sin = p5 / sqrt(p25 + p25 * _cot * _cot);
+          _cos = _sin * _cot;
+        } else {
+          _tan = sdiag[k] / r[kk];
+          _cos = p5 / sqrt(p25 + p25 * _tan * _tan);
+          _sin = _cos * _tan;
+        }
+
+        r[kk] = _cos * r[kk] + _sin * sdiag[k];
+        temp = _cos * wa[k] + _sin * qtbpj;
+        qtbpj = -_sin * wa[k] + _cos * qtbpj;
+        wa[k] = temp;
+
+        for (i = k + 1; i < n; ++i) {
+          temp = _cos * r[k * ldr + i] + _sin * sdiag[i];
+          sdiag[i] = -_sin * r[k * ldr + i] + _cos * sdiag[i];
+          r[k * ldr + i] = temp;
+        }
+      }
+      break;
+    }
+    
+    sdiag[j] = r[j * ldr + j];
+    r[j * ldr + j] = x[j];
+  }
+
+  nsing = n;
+  for (j = 0; j < n; ++j) {
+    if (sdiag[j] == 0.0 && nsing == n)
+      nsing = j;
+    if (nsing < n)
+      wa[j] = 0;
+  }
+
+  for (j = nsing - 1; j >= 0; --j) {
+    sum = 0;
+    for (i = j + 1; i < nsing; ++i)
+      sum += r[j * ldr + i] * wa[i];
+    wa[j] = (wa[j] - sum) / sdiag[j];
+  }
+
+  for (j = 0; j < n; ++j)
+    x[ipvt[j]] = wa[j];
+}
+
+
+// copied from the C code
+private void lm_lmpar(int n, real[] r, int ldr, int[] ipvt, real[] diag, real[] qtb, real delta, lm_real_type par, real[] x, real[] sdiag, real[] wa1, real[] wa2) {
+  static real p1 = 0.1;
+  static real p001 = 0.001;
+
+  int nsing = n;
+  real parl = 0.0;
+
+  int i, iter, j;
+  real dxnorm, fp, fp_old, gnorm, parc, paru;
+  real sum, temp;
+
+  for (j = 0; j < n; ++j) {
+    wa1[j] = qtb[j];
+    if (r[j * ldr + j] == 0 && nsing == n)
+      nsing = j;
+    if (nsing < n)
+      wa1[j] = 0;
+  }
+  for (j = nsing - 1; j >= 0; --j) {
+    wa1[j] = wa1[j] / r[j + ldr * j];
+    temp = wa1[j];
+    for (i = 0; i < j; ++i)
+      wa1[i] -= r[j * ldr + i] * temp;
+  }
+
+  for (j = 0; j < n; ++j)
+    x[ipvt[j]] = wa1[j];
+
+  iter = 0;
+  for (j = 0; j < n; ++j)
+    wa2[j] = diag[j] * x[j];
+  dxnorm = lm_enorm(n, wa2);
+  fp = dxnorm - delta;
+  if (fp <= p1 * delta) {
+    par.val = 0;
+    return;
+  }
+
+  if (nsing >= n) {
+    for (j = 0; j < n; ++j)
+      wa1[j] = diag[ipvt[j]] * wa2[ipvt[j]] / dxnorm;
+
+    for (j = 0; j < n; ++j) {
+      sum = 0.0;
+      for (i = 0; i < j; ++i)
+        sum += r[j * ldr + i] * wa1[i];
+      wa1[j] = (wa1[j] - sum) / r[j + ldr * j];
+    }
+    temp = lm_enorm(n, wa1);
+    parl = fp / delta / temp / temp;
+  }
+
+  for (j = 0; j < n; ++j) {
+    sum = 0;
+    for (i = 0; i <= j; ++i)
+      sum += r[j * ldr + i] * qtb[i];
+    wa1[j] = sum / diag[ipvt[j]];
+  }
+  gnorm = lm_enorm(n, wa1);
+  paru = gnorm / delta;
+  if (paru == 0.0)
+    paru = LM_DWARF / min(delta, p1);
+
+  par.val = max(par.val, parl);
+  par.val = min(par.val, paru);
+  if (par.val == 0.0)
+    par.val = gnorm / dxnorm;
+
+  for (;; ++iter) {
+    if (par.val == 0.0)
+      par.val = max(LM_DWARF, p001 * paru);
+    temp = sqrt(par.val);
+    for (j = 0; j < n; ++j)
+      wa1[j] = temp * diag[j];
+    lm_qrsolv(n, r, ldr, ipvt, wa1, qtb, x, sdiag, wa2);
+    for (j = 0; j < n; ++j)
+      wa2[j] = diag[j] * x[j];
+    dxnorm = lm_enorm(n, wa2);
+    fp_old = fp;
+    fp = dxnorm - delta;
+        
+    if (fabs(fp) <= p1 * delta || (parl == 0.0 && fp <= fp_old && fp_old < 0.0) || iter == 10)
+      break;
+        
+    for (j = 0; j < n; ++j)
+      wa1[j] = diag[ipvt[j]] * wa2[ipvt[j]] / dxnorm;
+
+    for (j = 0; j < n; ++j) {
+      wa1[j] = wa1[j] / sdiag[j];
+      for (i = j + 1; i < n; ++i)
+        wa1[i] -= r[j * ldr + i] * wa1[j];
+    }
+    temp = lm_enorm(n, wa1);
+    parc = fp / delta / temp / temp;
+    
+    if (fp > 0)
+      parl = max(parl, par.val);
+    else if (fp < 0)
+      paru = min(paru, par.val);
+    
+    par.val = max(parl, par.val + parc);
+  }
+}
+
+
+// copied from the C code; the main function
+void lm_lmdif(int m, int n, real[] x, real[] fvec, real ftol, real xtol, real gtol, int maxfev, real epsfcn, real[] diag, int mode, real factor, lm_int_type info, lm_int_type nfev, real[] fjac, int[] ipvt, real[] qtf, real[] wa1, real[] wa2, real[] wa3, real[] wa4, lm_evaluate_ftype evaluate, lm_print_ftype printout, lm_data_type data) {
+  static real p1 = 0.1;
+  static real p5 = 0.5;
+  static real p25 = 0.25;
+  static real p75 = 0.75;
+  static real p0001 = 1.0e-4;
+  
+  nfev.val = 0;
+  int iter = 1;
+  lm_real_type par = lm_real_type(0);
+  real delta = 0;
+  real xnorm = 0;
+  real temp = max(epsfcn, LM_MACHEP);
+  real eps = sqrt(temp);
+  int i, j;
+  real actred, dirder, fnorm, fnorm1, gnorm, pnorm, prered, ratio, step, sum, temp1, temp2, temp3;
+
+  if ((n <= 0) || (m < n) || (ftol < 0.0) || (xtol < 0.0) || (gtol < 0.0) || (maxfev <= 0) || (factor <= 0)) {
+    info.val = 0;
+    return;
+  }
+  if (mode == 2) {
+    for (j = 0; j < n; ++j) {
+      if (diag[j] <= 0.0) {
+        info.val = 0;
+        return;
+      }
+    }
+  }
+  
+  info.val = 0;
+  evaluate(x, m, fvec, data, info);
+  printout(n, x, m, fvec, data, 0, 0, ++nfev.val);
+  if (info.val < 0)
+    return;
+  fnorm = lm_enorm(m, fvec);
+
+  do {
+    for (j = 0; j < n; ++j) {
+      temp = x[j];
+      step = eps * fabs(temp);
+      if (step == 0.0)
+        step = eps;
+      x[j] = temp + step;
+      info.val = 0;
+      evaluate(x, m, wa4, data, info);
+      printout(n, x, m, wa4, data, 1, iter, ++nfev.val);
+      if (info.val < 0)
+        return;
+      for (i = 0; i < m; ++i)
+        fjac[j * m + i] = (wa4[i] - fvec[i]) / (x[j] - temp);
+      x[j] = temp;
+    }
+    
+    lm_qrfac(m, n, fjac, true, ipvt, wa1, wa2, wa3);
+
+    if (iter == 1) {
+      if (mode != 2) {
+        for (j = 0; j < n; ++j) {
+          diag[j] = wa2[j];
+          if (wa2[j] == 0.0)
+            diag[j] = 1.0;
+        }
+      }
+      for (j = 0; j < n; ++j)
+        wa3[j] = diag[j] * x[j];
+      xnorm = lm_enorm(n, wa3);
+      delta = factor * xnorm;
+      if (delta == 0.0)
+        delta = factor;
+    }
+
+    for (i = 0; i < m; ++i)
+      wa4[i] = fvec[i];
+
+    for (j = 0; j < n; ++j) {
+      temp3 = fjac[j * m + j];
+      if (temp3 != 0.0) {
+        sum = 0;
+        for (i = j; i < m; ++i)
+          sum += fjac[j * m + i] * wa4[i];
+        temp = -sum / temp3;
+        for (i = j; i < m; ++i)
+          wa4[i] += fjac[j * m + i] * temp;
+      }
+      fjac[j * m + j] = wa1[j];
+      qtf[j] = wa4[j];
+    }
+
+    gnorm = 0;
+    if (fnorm != 0) {
+      for (j = 0; j < n; ++j) {
+        if (wa2[ipvt[j]] == 0) continue;
+        sum = 0.0;
+        for (i = 0; i <= j; ++i)
+          sum += fjac[j * m + i] * qtf[i] / fnorm;
+        gnorm = max(gnorm, fabs(sum / wa2[ipvt[j]]));
+      }
+    }
+
+    if (gnorm <= gtol) {
+      info.val = 4;
+      return;
+    }
+
+    if (mode != 2) {
+      for (j = 0; j < n; ++j)
+        diag[j] = max(diag[j], wa2[j]);
+    }
+
+    do {
+      lm_lmpar(n, fjac, m, ipvt, diag, qtf, delta, par, wa1, wa2, wa3, wa4);
+
+      for (j = 0; j < n; ++j) {
+        wa1[j] = -wa1[j];
+        wa2[j] = x[j] + wa1[j];
+        wa3[j] = diag[j] * wa1[j];
+      }
+      pnorm = lm_enorm(n, wa3);
+
+      if (nfev.val <= 1 + n)
+        delta = min(delta, pnorm);
+
+      info.val = 0;
+      evaluate(wa2, m, wa4, data, info);
+      printout(n, x, m, wa4, data, 2, iter, ++nfev.val);
+      if (info.val < 0)
+        return;
+
+      fnorm1 = lm_enorm(m, wa4);
+
+      if (p1 * fnorm1 < fnorm)
+        actred = 1 - SQR(fnorm1 / fnorm);
+      else
+        actred = -1;
+
+      for (j = 0; j < n; ++j) {
+        wa3[j] = 0;
+        for (i = 0; i <= j; ++i)
+          wa3[i] += fjac[j * m + i] * wa1[ipvt[j]];
+      }
+      temp1 = lm_enorm(n, wa3) / fnorm;
+      temp2 = sqrt(par.val) * pnorm / fnorm;
+      prered = SQR(temp1) + 2 * SQR(temp2);
+      dirder = -(SQR(temp1) + SQR(temp2));
+
+      ratio = prered != 0 ? actred / prered : 0;
+
+      if (ratio <= p25) {
+        if (actred >= 0.0)
+          temp = p5;
+        else
+          temp = p5 * dirder / (dirder + p5 * actred);
+        if (p1 * fnorm1 >= fnorm || temp < p1)
+          temp = p1;
+        delta = temp * min(delta, pnorm / p1);
+        par.val /= temp;
+      } else if (par.val == 0.0 || ratio >= p75) {
+        delta = pnorm / p5;
+        par.val *= p5;
+      }
+      
+      if (ratio >= p0001) {
+        for (j = 0; j < n; ++j) {
+          x[j] = wa2[j];
+          wa2[j] = diag[j] * x[j];
+        }
+        for (i = 0; i < m; ++i)
+          fvec[i] = wa4[i];
+        xnorm = lm_enorm(n, wa2);
+        fnorm = fnorm1;
+        ++iter;
+      }
+
+      info.val = 0;
+      if (fabs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1)
+        info.val = 1;
+      if (delta <= xtol * xnorm)
+        info.val += 2;
+      if (info.val != 0)
+        return;
+
+      if (nfev.val >= maxfev)
+        info.val = 5;
+      if (fabs(actred) <= LM_MACHEP && prered <= LM_MACHEP && p5 * ratio <= 1)
+        info.val = 6;
+      if (delta <= LM_MACHEP * xnorm)
+        info.val = 7;
+      if (gnorm <= LM_MACHEP)
+        info.val = 8;
+      if (info.val != 0)
+        return;
+    } while (ratio < p0001);
+  } while (true);
+}
+
+
+// copied from the C code; wrapper of lm_lmdif
+void lm_minimize(int m_dat, int n_par, real[] par, lm_evaluate_ftype evaluate, lm_print_ftype printout, lm_data_type data, lm_control_type control) {
+  int n = n_par;
+  int m = m_dat;
+  
+  real[] fvec = new real[m];
+  real[] diag = new real[n];
+  real[] qtf = new real[n];
+  real[] fjac = new real[n * m];
+  real[] wa1 = new real[n];
+  real[] wa2 = new real[n];
+  real[] wa3 = new real[n];
+  real[] wa4 = new real[m];
+  int[] ipvt = new int[n];
+
+  control.info.val = 0;
+  control.nfev.val = 0;
+
+  lm_lmdif(m, n, par, fvec, control.ftol, control.xtol, control.gtol, control.maxcall * (n + 1), control.epsilon, diag, 1, control.stepbound, control.info, control.nfev, fjac, ipvt, qtf, wa1, wa2, wa3, wa4, evaluate, printout, data);
+
+  printout(n, par, m, fvec, data, -1, 0, control.nfev.val);
+  control.fnorm = lm_enorm(m, fvec);
+  if (control.info.val < 0)
+    control.info.val = 10;
+}
+
+
+// convenience functions; wrappers of lm_minimize
+struct FitControl {
+  real squareSumTolerance;
+  real approximationTolerance;
+  real desiredOrthogonality;
+  real epsilon;
+  real stepBound;
+  int maxIterations;
+  bool verbose;
+
+  void operator init(real squareSumTolerance, real approximationTolerance, real desiredOrthogonality, real epsilon, real stepBound, int maxIterations, bool verbose) {
+    this.squareSumTolerance = squareSumTolerance;
+    this.approximationTolerance = approximationTolerance;
+    this.desiredOrthogonality = desiredOrthogonality;
+    this.epsilon = epsilon;
+    this.stepBound = stepBound;
+    this.maxIterations = maxIterations;
+    this.verbose = verbose;
+  }
+
+  FitControl copy() {
+    FitControl result = new FitControl;
+    result.squareSumTolerance = this.squareSumTolerance;
+    result.approximationTolerance = this.approximationTolerance;
+    result.desiredOrthogonality = this.desiredOrthogonality;
+    result.epsilon = this.epsilon;
+    result.stepBound = this.stepBound;
+    result.maxIterations = this.maxIterations;
+    result.verbose = this.verbose;
+    return result;
+  }
+};
+
+FitControl defaultControl;
+defaultControl.squareSumTolerance = LM_USERTOL;
+defaultControl.approximationTolerance = LM_USERTOL;
+defaultControl.desiredOrthogonality = LM_USERTOL;
+defaultControl.epsilon = LM_USERTOL;
+defaultControl.stepBound = 100;
+defaultControl.maxIterations = 100;
+defaultControl.verbose = false;
+
+
+struct FitResult {
+  restricted real norm;
+  restricted int iterations;
+  restricted int status;
+
+  void operator init(real norm, int status, int iterations) {
+    this.norm = norm;
+    this.status = status;
+    this.iterations = iterations;
+  }
+};
+
+
+// Fits data points (xdata, ydata ± errors) to the given function using the given parameters.
+FitResult fit(real[] xdata, real[] ydata, real[] errors, real function(real[], real), real[] parameters, FitControl control=defaultControl) {
+  int m_dat = min(xdata.length, ydata.length);
+  int n_par = parameters.length;
+  lm_evaluate_ftype evaluate = lm_evaluate_default;
+  lm_print_ftype printout = control.verbose ? lm_print_default : lm_print_quiet;
+  
+  lm_data_type data;
+  data.user_t = xdata;
+  data.user_y = ydata;
+  data.user_w = 1 / errors;
+  data.user_func = new real(real x, real[] params) {
+    return function(params, x);
+  };
+
+  lm_control_type ctrl;
+  ctrl.ftol = control.squareSumTolerance;
+  ctrl.xtol = control.approximationTolerance;
+  ctrl.gtol = control.desiredOrthogonality;
+  ctrl.epsilon = control.epsilon;
+  ctrl.stepbound = control.stepBound;
+  ctrl.maxcall = control.maxIterations;
+
+  lm_minimize(m_dat, n_par, parameters, evaluate, printout, data, ctrl);
+  
+  return FitResult(ctrl.fnorm, ctrl.nfev.val, ctrl.info.val);
+}
+
+
+// Fits data points (xdata, ydata) to the given function using the given parameters.
+FitResult fit(real[] xdata, real[] ydata, real function(real[], real), real[] parameters, FitControl control=defaultControl) {
+  return fit(xdata, ydata, array(min(xdata.length, ydata.length), 1.0), function, parameters, control);
+}
+
diff --git a/Build/source/utils/asymptote/base/math.asy b/Build/source/utils/asymptote/base/math.asy
index f204f746349..85d41d86d07 100644
--- a/Build/source/utils/asymptote/base/math.asy
+++ b/Build/source/utils/asymptote/base/math.asy
@@ -148,7 +148,7 @@ real[][] zero(int n, int m)
   return M;
 }
 
-real[][] operator + (real[][] a, real[][] b)
+real[][] operator +(real[][] a, real[][] b)
 {
   int n=a.length;
   real[][] m=new real[n][];
@@ -157,7 +157,7 @@ real[][] operator + (real[][] a, real[][] b)
   return m;
 }
 
-real[][] operator - (real[][] a, real[][] b)
+real[][] operator -(real[][] a, real[][] b)
 {
   int n=a.length;
   real[][] m=new real[n][];
@@ -166,26 +166,7 @@ real[][] operator - (real[][] a, real[][] b)
   return m;
 }
 
-private string incommensurate=
-  "Multiplication of incommensurate matrices is undefined";
-
-real[] operator * (real[] b, real[][] a)
-{
-  int nb=b.length;
-  if(nb != a.length)
-    abort(incommensurate);
-  int na0=a[0].length;
-  real[] m=new real[na0];
-  for(int j=0; j < na0; ++j) {
-    real sum;
-    for(int k=0; k < nb; ++k)
-      sum += b[k]*a[k][j];
-    m[j]=sum;
-  }
-  return m;
-}
-
-real[][] operator * (real[][] a, real b)
+real[][] operator *(real[][] a, real b)
 {
   int n=a.length;
   real[][] m=new real[n][];
@@ -194,12 +175,12 @@ real[][] operator * (real[][] a, real b)
   return m;
 }
 
-real[][] operator * (real b, real[][] a)
+real[][] operator *(real b, real[][] a)
 {
   return a*b;
 }
 
-real[][] operator / (real[][] a, real b)
+real[][] operator /(real[][] a, real b)
 {
   return a*(1/b);
 }
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index 5135dd9443f..643c56feec7 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -1,15 +1,260 @@
-real euler(real y, real f(real x, real y), real a, real b=a, int n=0,
-           real h=0, bool dynamic=false, real tolmin=0, real tolmax=0)
+real stepfactor=2.0; // Maximum dynamic step size adjustment factor.
+
+struct RKTableau
+{
+  int order;
+  real[] steps;
+  real[][] weights;
+  real[] highOrderWeights;
+  real[] lowOrderWeights;
+  real pgrow;
+  real pshrink;
+
+  void operator init(int order, real[][] weights, real[] highOrderWeights,
+                     real[] lowOrderWeights=new real[],
+                     real[] steps=sequence(new real(int i) {
+                         return sum(weights[i]);},weights.length)) {
+    this.order=order;
+    this.steps=steps;
+    this.weights=weights;
+    this.highOrderWeights=highOrderWeights;
+    this.lowOrderWeights=lowOrderWeights;
+    pgrow=(order > 0) ? 1/order : 0;
+    pshrink=(order > 1) ? 1/(order-1) : pgrow;
+  }
+}
+
+// First-Order Euler
+RKTableau Euler=RKTableau(1,new real[][],
+                          new real[] {1});
+
+// Second-Order Runge-Kutta
+RKTableau RK2=RKTableau(2,new real[][] {{1/2}},
+                        new real[] {0,1});
+
+// Second-Order Predictor-Corrector
+RKTableau PC=RKTableau(2,new real[][] {{1}},
+                       new real[] {1/2,1/2});
+
+// Third-Order Classical Runge-Kutta
+RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}},
+                        new real[] {1/6,2/3,1/6});
+
+// Third-Order Bogacki-Shampine Runge-Kutta
+RKTableau RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}},
+                          new real[] {2/9,1/3,4/9}, // 3rd order
+                          new real[] {7/24,1/4,1/3,1/8}); // 2nd order
+
+// Fourth-Order Classical Runge-Kutta
+RKTableau RK4=RKTableau(4,new real[][] {{1/2},{0,1/2},{0,0,1}},
+                        new real[] {1/6,1/3,1/3,1/6});
+
+// Fifth-Order Cash-Karp Runge-Kutta
+RKTableau RK5CK=RKTableau(5,new real[][] {{1/5},
+                                          {3/40,9/40},
+                                          {3/10,-9/10,6/5},
+                                          {-11/54,5/2,-70/27,35/27},
+                                          {1631/55296,175/512,575/13824,
+                                           44275/110592,253/4096}},
+  new real[] {37/378,0,250/621,125/594,
+              0,512/1771},  // 5th order
+  new real[] {2825/27648,0,18575/48384,13525/55296,
+              277/14336,1/4}); // 4th order
+
+// Fifth-Order Fehlberg Runge-Kutta
+RKTableau RK5F=RKTableau(5,new real[][] {{1/4},
+                                         {3/32,9/32},
+                                         {1932/2197,-7200/2197,7296/2197},
+                                         {439/216,-8,3680/513,-845/4104},
+                                         {-8/27,2,-3544/2565,1859/4104,
+                                          -11/40}},
+  new real[] {16/135,0,6656/12825,28561/56430,-9/50,2/55}, // 5th order
+  new real[] {25/216,0,1408/2565,2197/4104,-1/5,0}); // 4th order
+
+// Fifth-Order Dormand-Prince Runge-Kutta
+RKTableau RK5DP=RKTableau(5,new real[][] {{1/5},
+                                          {3/40,9/40},
+                                          {44/45,-56/15,32/9},
+                                          {19372/6561,-25360/2187,64448/6561,
+                                           -212/729},
+                                          {9017/3168,-355/33,46732/5247,49/176,
+                                           -5103/18656}},
+  new real[] {35/384,0,500/1113,125/192,-2187/6784,
+              11/84}, // 5th order
+  new real[] {5179/57600,0,7571/16695,393/640,
+              -92097/339200,187/2100,1/40}); // 4th order
+
+real error(real error, real initial, real norm, real lowOrder, real diff) 
+{
+  if(initial != 0.0 && lowOrder != initial) {
+    static real epsilon=realMin/realEpsilon;
+    real denom=max(abs(norm),abs(initial))+epsilon;
+    return max(error,max(abs(diff)/denom));
+  }
+  return error;
+}
+
+real adjust(real h, real error, real t, real tolmin, real tolmax,
+            real dtmin, real dtmax, RKTableau tableau, bool verbose=true) 
+{
+  real dt=h;
+  void report(real t) {
+    if(h != dt)
+      write("Time step changed from "+(string) dt+" to "+(string) h+" at t="+
+            (string) t+".");
+  }
+  if(error > tolmax) {
+    h=max(h*max((tolmin/error)^tableau.pshrink,1/stepfactor),dtmin);
+    if(verbose) report(t);
+    return h;
+  }
+  if(error > 0 && error < tolmin) {
+    h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax);
+    if(verbose) report(t+dt);
+  }
+  return h;
+}
+
+// Integrate dy/dt=f(t,y) from a to b using initial conditions y,
+// specifying either the step size h or the number of steps n.
+real integrate(real y, real f(real t, real y), real a, real b=a, real h=0,
+               int n=0, bool dynamic=false, real tolmin=0, real tolmax=0,
+               real dtmin=0, real dtmax=realMax, RKTableau tableau,
+               bool verbose=false)
+{
+  if(h == 0) {
+    if(b == a) return y;
+    if(n == 0) abort("Either n or h must be specified");
+    else h=(b-a)/n;
+  }
+  real t=a;
+  real f0;
+  bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+  if(fsal) f0=f(t,y);
+  if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+  while(t < b) {
+    real[] predictions={fsal ? f0 : f(t,y)};
+    for(int i=0; i < tableau.steps.length; ++i)
+      predictions.push(f(t+h*tableau.steps[i],
+                         y+h*dot(tableau.weights[i],predictions)));
+
+    real highOrder=h*dot(tableau.highOrderWeights,predictions);
+    if(dynamic) {
+      real f1;
+      if(fsal) {
+        f1=f(t+h,y+highOrder);
+        predictions.push(f1);
+      }
+      real lowOrder=h*dot(tableau.lowOrderWeights,predictions);
+      real error;
+      error=error(error,y,y+highOrder,y+lowOrder,highOrder-lowOrder);
+      real dt=h;
+      h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+      if(h >= dt) {
+        t += dt;
+        y += highOrder;
+        f0=f1;
+      }
+    } else {
+      t += h;
+      y += highOrder;
+    }
+    h=min(h,b-t);
+    if(t >= b || t+h == t) break;
+  }
+  return y;
+}
+
+// Integrate a set of equations, dy/dt=f(t,y), from a to b using initial
+// conditions y, specifying either the step size h or the number of steps n.
+real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
+                 real h=0, int n=0, bool dynamic=false,
+                 real tolmin=0, real tolmax=0, real dtmin=0, real dtmax=realMax,
+                 RKTableau tableau, bool verbose=false)
 {
   if(h == 0) {
     if(b == a) return y;
     if(n == 0) abort("Either n or h must be specified");
     else h=(b-a)/n;
   }
-  real x=a;
-  for(int i=0; i < n; ++i) {
-    y += h*f(x,y);
-    x += h;
+  real[] y=copy(y);
+  real t=a;
+  real[] f0;
+  bool fsal=tableau.lowOrderWeights.length > tableau.highOrderWeights.length;
+  if(fsal) f0=f(t,y);
+  if(tableau.lowOrderWeights.length == 0) dynamic=false;
+
+  while(t < b) {
+    real[][] predictions={fsal ? f0 : f(t,y)};
+    for(int i=0; i < tableau.steps.length; ++i)
+      predictions.push(f(t+h*tableau.steps[i],
+                         y+h*tableau.weights[i]*predictions));
+
+    real[] highOrder=h*tableau.highOrderWeights*predictions;
+    if(dynamic) {
+      real[] f1;
+      if(fsal) {
+        f1=f(t+h,y+highOrder);
+        predictions.push(f1);
+      }
+      real[] lowOrder=h*tableau.lowOrderWeights*predictions;
+      real error;
+      for(int i=0; i < y.length; ++i)
+        error=error(error,y[i],y[i]+highOrder[i],y[i]+lowOrder[i],
+                    highOrder[i]-lowOrder[i]);
+      real dt=h;
+      h=adjust(h,error,t,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+      if(h >= dt) {
+        t += dt;
+        y += highOrder;
+        f0=f1;
+      }
+    } else {
+      t += h;
+      y += highOrder;
+    }
+    h=min(h,b-t);
+    if(t >= b || t+h == t) break;
   }
   return y;
 }
+
+real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t,
+                                  real[] h=sqrtEpsilon*abs(t))
+{
+  real[] ft=f(t);
+  real[][] J=new real[t.length][ft.length];
+  real[] ti=copy(t);
+  real tlast=ti[0];
+  ti[0] += h[0];
+  J[0]=(f(ti)-ft)/h[0];
+  for(int i=1; i < t.length; ++i) {
+    ti[i-1]=tlast;
+    tlast=ti[i];
+    ti[i] += h[i];
+    J[i]=(f(ti)-ft)/h[i];
+  }
+  return transpose(J);
+}
+
+// Solve simultaneous nonlinear system by Newton's method.
+real[] newton(int iterations=100, real[] f(real[]), real[][] jacobian(real[]),
+              real[] t)
+{
+  real[] t=copy(t);
+  for(int i=0; i < iterations; ++i)
+    t += solve(jacobian(t),-f(t));
+  return t;
+}
+
+real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0,
+                real[] initial(real[]), real[] discrepancy(real[]),
+                real[] guess, RKTableau tableau, int iterations=100)
+{
+  real[] g(real[] t) {
+    return discrepancy(integrate(initial(t),f,a,b,h,n,tableau));
+  }
+  real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);}
+  return initial(newton(iterations,g,jacobian,guess));
+}
diff --git a/Build/source/utils/asymptote/base/plain.asy b/Build/source/utils/asymptote/base/plain.asy
index a918a6bdc09..8424c767ffa 100644
--- a/Build/source/utils/asymptote/base/plain.asy
+++ b/Build/source/utils/asymptote/base/plain.asy
@@ -38,6 +38,8 @@ include plain_markers;
 include plain_arrows;
 include plain_debugger;
 
+pair viewportmargin=(0,0);   // Horizontal and vertical 3D viewport margins.
+
 typedef void exitfcn();
 
 bool needshipout() {
diff --git a/Build/source/utils/asymptote/base/plain_strings.asy b/Build/source/utils/asymptote/base/plain_strings.asy
index e3d0a797d86..0fe59ba5a63 100644
--- a/Build/source/utils/asymptote/base/plain_strings.asy
+++ b/Build/source/utils/asymptote/base/plain_strings.asy
@@ -181,9 +181,9 @@ string math(real x)
   return math((string) x);
 }
 
-string format(real x)
+string format(real x, string locale="")
 {
-  return format(defaultformat,x);
+  return format(defaultformat,x,locale);
 }
 
 string phantom(string s)
diff --git a/Build/source/utils/asymptote/base/slide.asy b/Build/source/utils/asymptote/base/slide.asy
index 5b3596d02e6..b154c1d4746 100644
--- a/Build/source/utils/asymptote/base/slide.asy
+++ b/Build/source/utils/asymptote/base/slide.asy
@@ -519,7 +519,6 @@ void bibliographystyle(string name)
 {
   settings.twice=true;
   settings.keepaux=true;
-  delete(outprefix()+"_.aux");
   texpreamble("\bibliographystyle{"+name+"}");
 }
 
diff --git a/Build/source/utils/asymptote/base/splinetype.asy b/Build/source/utils/asymptote/base/splinetype.asy
index 366f91ac370..08263ddf4ca 100644
--- a/Build/source/utils/asymptote/base/splinetype.asy
+++ b/Build/source/utils/asymptote/base/splinetype.asy
@@ -2,6 +2,7 @@ typedef real[] splinetype(real[], real[]);
 
 restricted real[] defaultspline(real[] x, real[] y);
 restricted real[] Spline(real[] x, real[] y);
+restricted splinetype[] Spline;
 
 string morepoints="interpolation requires at least 2 points";
 string differentlengths="arrays have different lengths";
@@ -63,7 +64,8 @@ real[] periodic(real[] x, real[] y)
 {
   int n=x.length;
   checklengths(n,y.length);
-  if(y[n-1] != y[0]) abort("function values are not periodic");
+  if(abs(y[n-1]-y[0]) > sqrtEpsilon*max(abs(y)))
+    abort("function values are not periodic");
   real[] d;
   if(n > 2) {
     real[] a=new real[n-1];
diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy
index 0af0d2c3bdf..d454d032b23 100644
--- a/Build/source/utils/asymptote/base/three.asy
+++ b/Build/source/utils/asymptote/base/three.asy
@@ -16,11 +16,10 @@ real defaultgranularity=0;
 real linegranularity=0.01;
 real tubegranularity=0.003;
 real dotgranularity=0.0001;
-pair viewportmargin=0;       // Horizontal and vertical viewport margins.
 real viewportfactor=1.002;   // Factor used to expand orthographic viewport.
-real viewportpadding=1.2;    // Offset used to expand PRC viewport.
+real viewportpadding=1;      // Offset used to expand PRC viewport.
 real angleprecision=1e-3;    // Precision for centering perspective projections.
-real anglefactor=max(1.005,1+angleprecision);
+real anglefactor=max(1.01,1+angleprecision);
 // Factor used to expand perspective viewport.
 
 string defaultembed3Doptions;
@@ -2178,7 +2177,7 @@ projection perspective(string s)
   return P;
 }
 
-private string format(real x)
+private string Format(real x)
 {
   // Work around movie15.sty division by zero bug;
   // e.g. u=unit((1e-10,1e-10,0.9));
@@ -2187,14 +2186,14 @@ private string format(real x)
   return format("%.18f",x,"C");
 }
 
-private string format(triple v, string sep=" ")
+private string Format(triple v, string sep=" ")
 {
-  return format(v.x)+sep+format(v.y)+sep+format(v.z);
+  return Format(v.x)+sep+Format(v.y)+sep+Format(v.z);
 }
 
-private string format(real[] c)
+private string Format(real[] c)
 {
-  return format((c[0],c[1],c[2]));
+  return Format((c[0],c[1],c[2]));
 }
 
 private string[] file3;
@@ -2227,8 +2226,8 @@ string lightscript(light light) {
     string Li="L"+string(i);
     real[] diffuse=light.diffuse[i];
     script += Li+"=scene.createLight();"+'\n'+
-      Li+".direction.set("+format(-light.position[i],",")+");"+'\n'+
-      Li+".color.set("+format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n';
+      Li+".direction.set("+Format(-light.position[i],",")+");"+'\n'+
+      Li+".color.set("+Format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n';
   }
   // Work around initialization bug in Adobe Reader 8.0:
   return script +"
@@ -2255,14 +2254,9 @@ void writeJavaScript(string name, string preamble, string script)
     file3.push(name);
 }
 
-pair viewportmargin(real width, real height) 
+pair viewportmargin(pair lambda)
 {
-  pair viewportmargin=viewportmargin;
-  real xmargin=viewportmargin.x;
-  real ymargin=viewportmargin.y;
-  if(xmargin <= 0) xmargin=max(0.5*(viewportsize.x-width),0);
-  if(ymargin <= 0) ymargin=max(0.5*(viewportsize.y-height),0);
-  return (xmargin,ymargin);
+  return maxbound(0.5*(viewportsize-lambda),viewportmargin);
 }
 
 string embed3D(string label="", string text=label, string prefix,
@@ -2284,7 +2278,7 @@ string embed3D(string label="", string text=label, string prefix,
   real viewplanesize;
   if(P.infinity) {
     triple lambda=max3(f)-min3(f);
-    pair margin=viewportpadding*(1,1)+viewportmargin(lambda.x,lambda.y);
+    pair margin=viewportpadding*(1,1)+viewportmargin((lambda.x,lambda.y));
     viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/cm;
   } else
     if(!P.absolute) angle=2*aTan(Tan(0.5*angle)-viewportpadding/P.target.z);
@@ -2314,12 +2308,12 @@ 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)+
-    ",3Dc2c="+format(u)+
-    ",3Dcoo="+format(P.target/cm)+
-    ",3Droll="+format(roll)+
-    ",3Droo="+format(abs(v))+
-    ",3Dbg="+format(light.background());
+    ",3Daac="+Format(P.absolute ? P.angle : angle)+
+    ",3Dc2c="+Format(u)+
+    ",3Dcoo="+Format(P.target/cm)+
+    ",3Droll="+Format(roll)+
+    ",3Droo="+Format(abs(v))+
+    ",3Dbg="+Format(light.background());
   if(options != "") options3 += ","+options;
   if(name != "") options3 += ",3Djscript="+stripdirectory(name);
 
@@ -2406,7 +2400,7 @@ object embed(string label="", string text=label,
     pair m2=pic2.min(s);
     pair M2=pic2.max(s);
     pair lambda=M2-m2;
-    pair viewportmargin=viewportmargin(lambda.x,lambda.y);
+    pair viewportmargin=viewportmargin(lambda);
     real width=ceil(lambda.x+2*viewportmargin.x);
     real height=ceil(lambda.y+2*viewportmargin.y);
 
@@ -2440,7 +2434,7 @@ object embed(string label="", string text=label,
         triple m=min3(f);
         triple M=max3(f);
         triple lambda=M-m;
-        viewportmargin=viewportmargin(lambda.x,lambda.y);
+        viewportmargin=viewportmargin((lambda.x,lambda.y));
         width=lambda.x+2*viewportmargin.x;
         height=lambda.y+2*viewportmargin.y;
 
@@ -2542,7 +2536,6 @@ object embed(string label="", string text=label,
       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");
diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy
index 3c49b382133..597104ad202 100644
--- a/Build/source/utils/asymptote/base/three_surface.asy
+++ b/Build/source/utils/asymptote/base/three_surface.asy
@@ -1191,11 +1191,11 @@ void label(picture pic=currentpicture, Label L, triple position,
   pic.add(new void(frame f, transform3 t, picture pic, projection P) {
       // Handle relative projected 3D alignments.
       Label L=L.copy();
+      triple v=t*position;
       if(!align.is3D && L.align.relative && L.align.dir3 != O &&
          determinant(P.t) != 0)
-        L.align(L.align.dir*unit(project(L.align.dir3,P.t))); 
+          L.align(L.align.dir*unit(project(v+L.align.dir3,P.t)-project(v,P.t)));
       
-      triple v=t*position;
       if(L.defaulttransform3)
         L.T3=transform3(P);
       if(is3D())
@@ -1234,7 +1234,7 @@ void label(picture pic=currentpicture, Label L, path3 g, align align=NoAlign,
   if(L.align.default) {
     align a;
     a.init(-I*(position <= sqrtEpsilon ? S :
-                   position >= length(g)-sqrtEpsilon ? N : E),relative=true);
+              position >= length(g)-sqrtEpsilon ? N : E),relative=true);
     a.dir3=dir(g,position); // Pass 3D direction via unused field.
     L.align(a);             
   }