summaryrefslogtreecommitdiff
path: root/Master/texmf/asymptote
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2009-08-15 23:57:48 +0000
committerKarl Berry <karl@freefriends.org>2009-08-15 23:57:48 +0000
commit16d128e5e10d541a78654b86409d5a3539f07708 (patch)
tree66de0af63c3811bb3040c16e9b52c11985f70811 /Master/texmf/asymptote
parentb20f78c549859ec0e8610bdd3ad904245e86b489 (diff)
asymptote 1.83
git-svn-id: svn://tug.org/texlive/trunk@14696 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote')
-rwxr-xr-xMaster/texmf/asymptote/GUI/xasyVersion.py2
-rw-r--r--Master/texmf/asymptote/animation.asy71
-rw-r--r--Master/texmf/asymptote/bezulate.asy2
-rw-r--r--Master/texmf/asymptote/geometry.asy56
-rw-r--r--Master/texmf/asymptote/graph.asy4
-rw-r--r--Master/texmf/asymptote/graph3.asy208
-rw-r--r--Master/texmf/asymptote/graph_splinetype.asy251
-rw-r--r--Master/texmf/asymptote/interpolate.asy2
-rw-r--r--Master/texmf/asymptote/labelpath.asy2
-rw-r--r--Master/texmf/asymptote/labelpath3.asy4
-rw-r--r--Master/texmf/asymptote/obj.asy6
-rw-r--r--Master/texmf/asymptote/plain.asy9
-rw-r--r--Master/texmf/asymptote/plain_Label.asy8
-rw-r--r--Master/texmf/asymptote/plain_arrows.asy47
-rw-r--r--Master/texmf/asymptote/plain_boxes.asy4
-rw-r--r--Master/texmf/asymptote/plain_constants.asy7
-rw-r--r--Master/texmf/asymptote/plain_paths.asy2
-rw-r--r--Master/texmf/asymptote/plain_pens.asy4
-rw-r--r--Master/texmf/asymptote/plain_picture.asy106
-rw-r--r--Master/texmf/asymptote/plain_shipout.asy10
-rw-r--r--Master/texmf/asymptote/slide.asy75
-rw-r--r--Master/texmf/asymptote/solids.asy3
-rw-r--r--Master/texmf/asymptote/syzygy.asy10
-rw-r--r--Master/texmf/asymptote/three.asy719
-rw-r--r--Master/texmf/asymptote/three_arrows.asy84
-rw-r--r--Master/texmf/asymptote/three_light.asy113
-rw-r--r--Master/texmf/asymptote/three_surface.asy812
-rw-r--r--Master/texmf/asymptote/three_tube.asy373
-rw-r--r--Master/texmf/asymptote/tree.asy4
-rw-r--r--Master/texmf/asymptote/trembling.asy2
-rw-r--r--Master/texmf/asymptote/tube.asy63
-rw-r--r--Master/texmf/asymptote/version.asy2
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";