summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2010-04-19 23:59:33 +0000
committerKarl Berry <karl@freefriends.org>2010-04-19 23:59:33 +0000
commitdddce3148a2bc785603576b18ffcf72d39adbe47 (patch)
tree14f4c79ef11bcf0820a8d4ff2ab0c5ac009c365c /Build/source/utils/asymptote/base
parenta01e51b01f5819b6091af48cdca581e9f2a9282e (diff)
asy 1.93
git-svn-id: svn://tug.org/texlive/trunk@17934 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r--Build/source/utils/asymptote/base/CAD.asy10
-rw-r--r--Build/source/utils/asymptote/base/animation.asy3
-rw-r--r--Build/source/utils/asymptote/base/contour.asy1706
-rw-r--r--Build/source/utils/asymptote/base/contour3.asy2
-rw-r--r--Build/source/utils/asymptote/base/embed.asy3
-rw-r--r--Build/source/utils/asymptote/base/external.asy7
-rw-r--r--Build/source/utils/asymptote/base/flowchart.asy202
-rw-r--r--Build/source/utils/asymptote/base/geometry.asy13
-rw-r--r--Build/source/utils/asymptote/base/graph.asy73
-rw-r--r--Build/source/utils/asymptote/base/graph3.asy30
-rw-r--r--Build/source/utils/asymptote/base/graph_splinetype.asy13
-rw-r--r--Build/source/utils/asymptote/base/math.asy50
-rw-r--r--Build/source/utils/asymptote/base/ode.asy73
-rw-r--r--Build/source/utils/asymptote/base/palette.asy8
-rw-r--r--Build/source/utils/asymptote/base/plain.asy2
-rw-r--r--Build/source/utils/asymptote/base/plain_Label.asy17
-rw-r--r--Build/source/utils/asymptote/base/plain_arcs.asy5
-rw-r--r--Build/source/utils/asymptote/base/plain_boxes.asy7
-rw-r--r--Build/source/utils/asymptote/base/plain_constants.asy4
-rw-r--r--Build/source/utils/asymptote/base/plain_paths.asy46
-rw-r--r--Build/source/utils/asymptote/base/plain_pens.asy17
-rw-r--r--Build/source/utils/asymptote/base/plain_picture.asy96
-rw-r--r--Build/source/utils/asymptote/base/plain_shipout.asy5
-rw-r--r--Build/source/utils/asymptote/base/plain_strings.asy43
-rw-r--r--Build/source/utils/asymptote/base/slide.asy4
-rw-r--r--Build/source/utils/asymptote/base/solids.asy2
-rw-r--r--Build/source/utils/asymptote/base/three.asy92
-rw-r--r--Build/source/utils/asymptote/base/three_arrows.asy21
-rw-r--r--Build/source/utils/asymptote/base/three_surface.asy110
29 files changed, 1932 insertions, 732 deletions
diff --git a/Build/source/utils/asymptote/base/CAD.asy b/Build/source/utils/asymptote/base/CAD.asy
index 7c333318f5c..b65e6478838 100644
--- a/Build/source/utils/asymptote/base/CAD.asy
+++ b/Build/source/utils/asymptote/base/CAD.asy
@@ -106,12 +106,12 @@ struct sCAD
// E
cad.pE =
cad.pSurfaceTreatmentAllowed =
- pFullWidth + linetype("10 2.5");
+ pFullWidth + linetype(new real[] {10,2.5});
// F
cad.pF =
cad.pInvisibleEdge =
cad.pInvisibleContour =
- pHalfWidth + linetype("20 5");
+ pHalfWidth + linetype(new real[] {20,5});
// G
cad.pG =
cad.pMiddleLine =
@@ -120,7 +120,7 @@ struct sCAD
cad.pCircularHole =
cad.pDivisionPlane =
cad.pTransferLine =
- pHalfWidth + linetype("40 5 5 5");
+ pHalfWidth + linetype(new real[] {40,5,5,5});
// H
// see J
// I
@@ -129,7 +129,7 @@ struct sCAD
cad.pJ =
cad.pCuttingPlane =
cad.pSurfaceTreatmentRequested =
- pFullWidth + linetype("20 2.5 2.5 2.5");
+ pFullWidth + linetype(new real[] {20,2.5,2.5,2.5});
// K
cad.pK =
cad.pContourBeforeDeformation =
@@ -137,7 +137,7 @@ struct sCAD
cad.pEndShapeRawMaterial =
cad.pContourEligibleType =
cad.pPartInFrontOfCuttingPlane =
- pHalfWidth + linetype("40 5 5 5 5 5");
+ pHalfWidth + linetype(new real[] {40,5,5,5,5,5});
return cad;
} // end of Create
diff --git a/Build/source/utils/asymptote/base/animation.asy b/Build/source/utils/asymptote/base/animation.asy
index 27f7171a639..582b32bced4 100644
--- a/Build/source/utils/asymptote/base/animation.asy
+++ b/Build/source/utils/asymptote/base/animation.asy
@@ -33,8 +33,7 @@ struct animation {
// been generated.
void operator init(string prefix="", bool global=true) {
- prefix=replace((prefix == "") ? outprefix() : stripdirectory(prefix),
- " ","_");
+ prefix=replace(stripdirectory(outprefix(prefix))," ","_");
this.prefix=prefix;
this.global=global;
}
diff --git a/Build/source/utils/asymptote/base/contour.asy b/Build/source/utils/asymptote/base/contour.asy
index fbb4cd1c71a..3c01640b2cc 100644
--- a/Build/source/utils/asymptote/base/contour.asy
+++ b/Build/source/utils/asymptote/base/contour.asy
@@ -1,307 +1,923 @@
-// Contour routines written by Radoslav Marinov and John Bowman.
-
+/*
+ Contour routines written by Radoslav Marinov, John Bowman, and Chris Savage.
+
+ [2009/10/15: C Savage] generate oriented contours
+ [2009/10/19: C Savage] use boxes instead of triangles
+*/
+
+/*
+ Contours lines/guides are oriented throughout. By convention,
+ for a single contour, higher values are to the left and/or lower
+ values are to the right along the direction of the lines/guide.
+*/
+
import graph_settings;
-real eps=10000*realEpsilon;
-
-// 1
-// 6 +-------------------+ 5
-// | \ / |
-// | \ / |
-// | \ / |
-// | \ / |
-// 2 | X | 0
-// | / \ |
-// | / \ |
-// | / \ |
-// | / \ |
-// 7 +-------------------+ 4 or 8
-// 3
+real eps=sqrtEpsilon;
-private struct segment
-{
- bool active;
- pair a,b; // Endpoints; a is always an edge point if one exists.
- int c; // Contour value.
- int edge; // -1: interior, 0 to 3: edge,
- // 4-8: single-vertex edge, 9: double-vertex edge.
+/*
+ GRID CONTOURS
+
+ Contours on a grid of points are determined as follows:
+ for each grid square, the function is approximated as the unique
+ paraboloid passing through the function values at the four
+ corners. The intersection of a paraboloid with the f(x,y)=c
+ plane is a line or hyperbola.
+
+ Grid data structures:
+
+ boxcontour:
+ Describes a particular contour segment in a grid square.
+
+ boxdata:
+ Describes contours in a grid square (holds boxcontours).
+
+ segment:
+ Describes a contour line. Usually a closed (interior) contour,
+ a line that terminates on the border, or a border segment used
+ to enclose a region.
+
+ Segment:
+ Describes a contour line.
+
+ Main grid routines:
+
+ setcontour:
+ Determines the contours in a grid square.
+
+ contouredges:
+ Determines the contour segments over a grid of function values.
+
+ connect:
+ Converts contours into guides
+
+*/
+
+typedef int boxtype;
+boxtype exterior=-1;
+boxtype edge = 0;
+boxtype interior=+1;
+
+typedef int contourshape;
+contourshape none =0;
+contourshape line =1;
+contourshape hyperbola=2;
+
+// Describe position by grid square and position in square
+private struct gridpoint {
+ int i,j;
+ pair z;
+ void operator init(int i, int j, pair z) {
+ this.i=i;
+ this.j=j;
+ this.z=z;
+ }
+ void operator init(gridpoint gp) {
+ this.i=gp.i;
+ this.j=gp.j;
+ this.z=gp.z;
+ }
}
-// Case 1: line passes through two vertices of a triangle
-private segment case1(pair p0, pair p1, int edge)
+private bool same(gridpoint gp1, gridpoint gp2)
{
- // Will cause a duplicate guide; luckily case1 is rare
- segment rtrn;
- rtrn.active=true;
- rtrn.a=p0;
- rtrn.b=p1;
- rtrn.edge=edge;
- return rtrn;
+ return abs(gp2.z-gp1.z+(gp2.i-gp1.i,gp2.j-gp1.j)) < eps;
}
-// Case 2: line passes through a vertex and a side of a triangle
-// (the first vertex passed and the side between the other two)
-private segment case2(pair p0, pair p1, pair p2,
- real v0, real v1, real v2, int edge)
-{
- segment rtrn;
- pair val=interp(p1,p2,abs(v1/(v2-v1)));
- rtrn.active=true;
- if(edge < 4) {
- rtrn.a=val;
- rtrn.b=p0;
- } else {
- rtrn.a=p0;
- rtrn.b=val;
+
+// Describe contour in unit square(scaling to be done later).
+private struct boxcontour {
+ bool active;
+ contourshape type; // Shape of contour segment(line or hyperbola)
+ pair a,b; // Start/end point of contour segment.
+ // Higher values to left along a--b.
+ real x0,y0,m; // For hyperbola: (x-x0)*(y-y0)=m
+ int signx,signy; // Sign of x-x0&y-y0 for hyperbola piece;
+ // identifies which direction it opens
+ int i,j; // Indices of lower left corner in position or
+ // data array.
+ int index; // Contour index
+
+ void operator init(contourshape type, pair a, pair b,
+ real x0, real y0, real m, int signx, int signy,
+ int i, int j, int index) {
+ this.active=true;
+ this.type=type;
+ this.a=a;
+ this.b=b;
+
+ this.x0=x0;
+ this.y0=y0;
+ this.m=m;
+ this.signx=signx;
+ this.signy=signy;
+
+ this.i=i;
+ this.j=j;
+ this.index=index;
+ }
+ // Generate list of points along the line/hyperbola segment
+ // representing the contour in the box
+ gridpoint[] points(int subsample=1, bool first=true, bool last=true) {
+ gridpoint[] gp;
+ if(first)
+ gp.push(gridpoint(i,j,a));
+ if(subsample > 0) {
+ // Linear case
+ if(type == line) {
+ for(int k=1; k <= subsample; ++k) {
+ pair z=interp(a,b,k/(subsample+1));
+ gp.push(gridpoint(i,j,z));
+ }
+ } else if(type == hyperbola) {
+ // Special hyperbolic case of m=0
+ // The contours here are infinite lines at x=x0 and y=y0,
+ // but handedness always connects a semi-infinite
+ // horizontal segment with a semi-infinite vertical segment
+ // connected at (x0,y0).
+ // If (x0,y0) is outside the unit box, there is only one
+ // line segment to include; otherwise, there are both
+ // a horizontal and a vertical line segment to include.
+ if(m == 0) {
+ // Single line
+ if(a.x == b.x || a.y == b.y) {
+ for(int k=1; k <= subsample; ++k) {
+ pair z=interp(a,b,k/(subsample+1));
+ gp.push(gridpoint(i,j,z));
+ }
+ // Two lines(may get one extra point here)
+ } else {
+ int nsub=quotient(subsample,2);
+ pair mid=(x0,y0);
+ for(int k=1; k <= nsub; ++k) {
+ pair z=interp(a,mid,k/(nsub+1));
+ gp.push(gridpoint(i,j,z));
+ }
+ gp.push(gridpoint(i,j,mid));
+ for(int k=1; k <= nsub; ++k) {
+ pair z=interp(mid,b,k/(nsub+1));
+ gp.push(gridpoint(i,j,z));
+ }
+ }
+ // General hyperbolic case (m != 0).
+ // Parametric equations(m > 0):
+ // x(t)=x0 +/- sqrt(m)*exp(t)
+ // y(t)=y0 +/- sqrt(m)*exp(-t)
+ // Parametric equations (m < 0):
+ // x(t)=x0 +/- sqrt(-m)*exp(t)
+ // y(t)=y0 -/+ sqrt(-m)*exp(-t)
+ // Points will be taken equally spaced in parameter t.
+ } else {
+ real sqrtm=sqrt(abs(m));
+ real ta=log(signx*(a.x-x0)/sqrtm);
+ real tb=log(signx*(b.x-x0)/sqrtm);
+ real[] t=uniform(ta,tb,subsample+1);
+ for(int k=1; k <= subsample; ++k) {
+ pair z=(x0+signx*sqrtm*exp(t[k]),
+ y0+signy*sqrtm*exp(-t[k]));
+ gp.push(gridpoint(i,j,z));
+ }
+ }
+ }
+ }
+ if(last)
+ gp.push(gridpoint(i,j,b));
+
+ return gp;
}
- rtrn.edge=edge;
- return rtrn;
}
-// Case 3: line passes through two sides of a triangle
-// (through the sides formed by the first & second, and second & third
-// vertices)
-private segment case3(pair p0, pair p1, pair p2,
- real v0, real v1, real v2, int edge=-1)
-{
- segment rtrn;
- rtrn.active=true;
- rtrn.a=interp(p1,p0,abs(v1/(v0-v1)));
- rtrn.b=interp(p1,p2,abs(v1/(v2-v1)));
- rtrn.edge=edge;
- return rtrn;
+// Hold data for a single grid square
+private struct boxdata {
+ boxtype type; // Does box contain a contour segment (edge of
+ // contour region) or is it entirely interior/
+ // exterior to contour region ?
+ real min,max; // Smallest/largest corner value
+ real max2; // Second-largest corner value
+ boxcontour[] data; // Stores actual contour segment data
+
+ int count() {return data.length;}
+ void operator init(real f00, real f10, real f01, real f11) {
+ real[] X={f00,f10,f01,f11};
+ min=min(X);
+ max=max(X);
+ X.delete(find(X == max));
+ max2=max(X);
+ }
+ void settype(real c) {
+ // Interior case(f >= c)
+ if(min > c-eps) {
+ type=interior;
+ // Exterior case(f < c)
+ } else if(max < c-eps) {
+ type=exterior;
+ // Special case: only one corner at f=c, f < c elsewhere
+ //(no segment in this case)
+ } else if((max < c+eps) && (max2 < c-eps)) {
+ type=exterior;
+ // Edge of contour passes through box
+ } else {
+ type=edge;
+ }
+ }
}
-// Check if a line passes through a triangle, and draw the required line.
-private segment checktriangle(pair p0, pair p1, pair p2,
- real v0, real v1, real v2, int edge=-1)
-{
- // default null return
- static segment dflt;
- real eps=eps*max(abs(v0),abs(v1),abs(v2));
+/*
+ Determine contours within a unit square box.
- if(v0 < -eps) {
- if(v1 < -eps) {
- if(v2 < -eps) return dflt; // nothing to do
- else if(v2 <= eps) return dflt; // nothing to do
- else return case3(p0,p2,p1,v0,v2,v1);
- } else if(v1 <= eps) {
- if(v2 < -eps) return dflt; // nothing to do
- else if(v2 <= eps) return case1(p1,p2,5+edge);
- else return case2(p1,p0,p2,v1,v0,v2,5+edge);
- } else {
- if(v2 < -eps) return case3(p0,p1,p2,v0,v1,v2,edge);
- else if(v2 <= eps)
- return case2(p2,p0,p1,v2,v0,v1,edge);
- else return case3(p1,p0,p2,v1,v0,v2,edge);
- }
- } else if(v0 <= eps) {
- if(v1 < -eps) {
- if(v2 < -eps) return dflt; // nothing to do
- else if(v2 <= eps) return case1(p0,p2,4+edge);
- else return case2(p0,p1,p2,v0,v1,v2,4+edge);
- } else if(v1 <= eps) {
- if(v2 < -eps) return case1(p0,p1,9);
- else if(v2 <= eps) return dflt; // use finer partitioning.
- else return case1(p0,p1,9);
+ Here, we approximate the function on the unit square to be a quadric
+ surface passing through the specified values at the four corners:
+ f(x,y)=(1-x)(1-y) f00+x(1-y) f10+(1-x)y f01+xy f11
+ =a0+ax x+ay y+axy xy
+ where f00, f10, f01 and f11 are the function values at the four
+ corners of the unit square 0 < x < 1&0 < y < 1 and:
+ a0 =f00
+ ax =f10-f00
+ ay =f01-f00
+ axy=f00+f11-f10-f01
+ This can also be expressed in paraboloid form as:
+ f(x,y)=alpha [(x+y-cp)^2-(x-y-cn)^2]+d
+ where:
+ alpha=axy/4
+ cp =-(ax+ay)/a11
+ cn =-(ax-ay)/a11
+ d =(a0 axy-ax ay)/axy
+ In the procedure below, we take f00 - > f00-c etc. for a contour
+ level c and we search for f=0.
+
+ For this surface, there are two possible contour shapes:
+ linear: (y-y0)/(x-x0)=m
+ hyperbolic: (x-x0)*(y-y0)=m
+ The linear case has a single line. The hyperbolic case may have
+ zero, one or two segments within the box (there are two sides of
+ a hyperbola, each of which may or may not pass through the unit
+ square). A hyperbola with m=0 is a special case that is handled
+ separately below.
+
+ If c0 is the desired contour level, we effectively find the
+ contours at c0-epsilon for arbitrarily small epsilon. Flat
+ regions equal to c0 are considered to be interior to the
+ contour curves. Regions that lie at the contour level are
+ considered to be interior to the contour curves. As a result,
+ contours are only constructed if they are immediately adjacent
+ to some region interior to the square that falls below the
+ contour value; in other words, if an edge falls on the contour
+ value, but a point within the square arbitrarily close to the
+ edge falls above the contour value, that edge (or applicable
+ portion) is not included. This requirement gives the following:
+ *) ensures contours on an edge are unique (do not appear in
+ an adjacent square with the same orientation)
+ *) no three line vertices (four line vertices are possible, but
+ are not usually an issue)
+ *) all segments can be joined into closed curves or curves that
+ terminate on the boundary (no unclosed curves terminate in
+ the interior region of the grid)
+
+ Note the logic below skips cases that have been filtered out
+ by the boxdata.settype() routine.
+*/
+private void setcontour(real f00, real f10, real f01, real f11, real epsf,
+ boxdata bd, int i, int j, int index) {
+ // SPECIAL CASE: two diagonal corners at the contour level with
+ // the other two below does not yield any contours within the
+ // unit box, but may have been previously misidentified as an
+ // edge containing region.
+ if(((f00*f11 == 0) && (f10*f01 > 0)) || ((f01*f10 == 0) && (f00*f11 > 0))) {
+ bd.type=exterior;
+ return;
+ }
+
+ // NOTE: From this point on, we can assume at least one contour
+ // segment exists in the square. This allows several cases to
+ // be ignored or simplified below, particularly edge cases.
+
+ // Form used to approximate function on unit square
+ real F(real x, real y) {
+ return interp(interp(f00,f10,x),interp(f01,f11,x),y);
+ }
+
+ // Write contour as a0+ax*x+ay*y +axy*x*y=0
+ real a0 =f00;
+ real ax =f10-f00;
+ if(abs(ax) < epsf) ax=0;
+ real ay =f01-f00;
+ if(abs(ay) < epsf) ay=0;
+ real axy=f00+f11-f01 -f10;
+ if(abs(axy) < epsf) axy=0;
+
+ // Linear contour(s)
+ if(axy == 0) {
+ pair a,b;
+ // Horizontal
+ if(ax == 0) {
+ if(ay == 0) return; // Contour is at most an isolated point; ignore.
+ real y0=-a0/ay;
+ if(abs(y0-1) < eps) y0=1;
+ if((f00 > 0) || (f01 < 0)) {
+ a=(1,y0);
+ b=(0,y0);
+ } else {
+ a=(0,y0);
+ b=(1,y0);
+ }
+ // Vertical
+ } else if(ay == 0) {
+ real x0=-a0/ax;
+ if(abs(x0-1) < eps) x0=1;
+ if((f00 > 0) || (f10 < 0)) {
+ a=(x0,0);
+ b=(x0,1);
+ } else {
+ a=(x0,1);
+ b=(x0,0);
+ }
+ // Angled line
} else {
- if(v2 < -eps) return case2(p0,p1,p2,v0,v1,v2,4+edge);
- else if(v2 <= eps) return case1(p0,p2,4+edge);
- else return dflt; // nothing to do
- }
- } else {
- if(v1 < -eps) {
- if(v2 < -eps) return case3(p1,p0,p2,v1,v0,v2,edge);
- else if(v2 <= eps)
- return case2(p2,p0,p1,v2,v0,v1,edge);
- else return case3(p0,p1,p2,v0,v1,v2,edge);
- } else if(v1 <= eps) {
- if(v2 < -eps) return case2(p1,p0,p2,v1,v0,v2,5+edge);
- else if(v2 <= eps) return case1(p1,p2,5+edge);
- else return dflt; // nothing to do
+ real x0=-a0/ax;
+ if(abs(x0-1) < eps) x0=1;
+ real y0=-a0/ay;
+ if(abs(y0-1) < eps) y0=1;
+ int count=0;
+ real[] farr={f00,f10,f11,f01};
+ farr.cyclic=true;
+ pair[] corners={(0,0),(1,0),(1,1),(0,1)};
+ pair[] sidedir={(1,0),(0,1),(-1,0),(0,-1)};
+
+ int count=0;
+ for(int i=0; i < farr.length; ++i) {
+ // Corner
+ if(farr[i] == 0) {
+ ++count;
+ if(farr[i-1] > 0) {
+ a=corners[i];
+ } else {
+ b=corners[i];
+ }
+ // Side
+ } else if(farr[i]*farr[i+1] < 0) {
+ ++count;
+ if(farr[i] > 0) {
+ a=corners[i]-(farr[i]/(farr[i+1]-farr[i]))*sidedir[i];
+ } else {
+ b=corners[i]-(farr[i]/(farr[i+1]-farr[i]))*sidedir[i];
+ }
+ }
+ }
+ // Check(if logic is correct above, this will not happen)
+ if(count != 2) {
+ abort("Unexpected error in setcontour routine: odd number of"
+ +" crossings (linear case)");
+ }
+ }
+ boxcontour bc=boxcontour(line,a,b,0,0,0,1,1,i,j,index);
+ bd.data.push(bc);
+ return;
+ }
+
+ // Hyperbolic contour(s)
+ // Described in form: (x-x0)*(y-y0)=m
+ real x0=-ay/axy;
+ if(abs(x0-1) < eps) x0=1;
+ real y0=-ax/axy;
+ if(abs(y0-1) < eps) y0=1;
+ real m =ay*ax-a0*axy;
+ m=(abs(m) < eps) ? 0 : m/axy^2;
+
+ // Special case here: straight segments (possibly crossing)
+ if(m == 0) {
+ pair a,b;
+ int signx,signy;
+ // Assuming at least one corner is below contour level here
+ if(x0 == 0) {
+ signx=+1;
+ if(y0 == 0) {
+ a=(1,0);
+ b=(0,1);
+ signy=+1;
+ } else if(y0 == 1) {
+ a=(0,0);
+ b=(1,1);
+ signy=-1;
+ } else if(y0 < 0 || y0 > 1) {
+ a=(0,0);
+ b=(0,1);
+ signy=y0 > 0 ? -1 : +1;
+ } else {
+ if(f10 > 0) {
+ a=(1,y0);
+ b=(0,1);
+ signy=+1;
+ } else {
+ a=(0,0);
+ b=(1,y0);
+ signy=-1;
+ }
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else if(x0 == 1) {
+ signx=-1;
+ if(y0 == 0) {
+ a=(1,1);
+ b=(0,0);
+ signy=+1;
+ } else if(y0 == 1) {
+ a=(0,1);
+ b=(1,0);
+ signy=-1;
+ } else if(y0 < 0 || y0 > 1) {
+ a=(1,1);
+ b=(1,0);
+ signy=y0 > 0 ? -1 : +1;
+ } else {
+ if(f01 > 0) {
+ a=(0,y0);
+ b=(1,0);
+ signy=-1;
+ } else {
+ a=(1,1);
+ b=(0,y0);
+ signy=+1;
+ }
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else if(y0 == 0) {
+ signy=+1;
+ if(x0 < 0 || x0 > 1) {
+ a=(1,0);
+ b=(0,0);
+ signx=x0 > 0 ? -1 : +1;
+ } else {
+ if(f11 > 0) {
+ a=(x0,1);
+ b=(0,0);
+ signx=-1;
+ } else {
+ a=(1,0);
+ b=(x0,1);
+ signx=+1;
+ }
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else if(y0 == 1) {
+ signy=-1;
+ if(x0 < 0 || x0 > 1) {
+ a=(0,1);
+ b=(1,1);
+ signx=x0 > 0 ? -1 : +1;
+ } else {
+ if(f00 > 0) {
+ a=(x0,0);
+ b=(1,1);
+ signx=+1;
+ } else {
+ a=(0,1);
+ b=(x0,0);
+ signx=-1;
+ }
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else if(x0 < 0 || x0 > 1) {
+ signx=x0 > 0 ? -1 : +1;
+ if(f00 > 0) {
+ a=(1,y0);
+ b=(0,y0);
+ signy=+1;
+ } else {
+ a=(0,y0);
+ b=(1,y0);
+ signy=-1;
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else if(y0 < 0 || y0 > 1) {
+ signy=y0 > 0 ? -1 : +1;
+ if(f00 > 0) {
+ a=(x0,0);
+ b=(x0,1);
+ signx=+1;
+ } else {
+ a=(x0,1);
+ b=(x0,0);
+ signx=-1;
+ }
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,signx,signy,i,j,index);
+ bd.data.push(bc);
+ return;
} else {
- if(v2 < -eps) return case3(p0,p2,p1,v0,v2,v1);
- else if(v2 <= eps) return dflt; // nothing to do
- else return dflt; // nothing to do
- }
- }
-}
+ if(f10 > 0) {
+ a=(0,y0);
+ b=(x0,0);
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,-1,-1,i,j,index);
+ bd.data.push(bc);
+ a=(1,y0);
+ b=(x0,1);
+ bc=boxcontour(hyperbola,a,b,x0,y0,m,+1,+1,i,j,index);
+ bd.data.push(bc);
+ return;
+ } else {
+ a=(x0,0);
+ b=(1,y0);
+ boxcontour bc=boxcontour(hyperbola,a,b,x0,y0,m,+1,-1,i,j,index);
+ bd.data.push(bc);
+ a=(x0,1);
+ b=(0,y0);
+ bc=boxcontour(hyperbola,a,b,x0,y0,m,-1,+1,i,j,index);
+ bd.data.push(bc);
+ return;
+ }
+ }
+ }
+
+ // General hyperbola case
+ int signc=(F(x0,y0) > 0) ? +1 : -1;
+
+ pair[] points;
+
+ real xB=(y0 == 0) ? infinity : x0-m/y0;
+ if(abs(xB) < eps) xB=0;
+ if(xB >= 0 && xB <= 1-eps) points.push((xB,0));
-// Collect connecting path segments.
-private void collect(pair[][][] points, real[] c)
-{
- // use to reverse an array, omitting the first point
- int[] reverseF(int n) {return sequence(new int(int x){return n-1-x;},n-1);}
- // use to reverse an array, omitting the last point
- int[] reverseL(int n) {return sequence(new int(int x){return n-2-x;},n-1);}
+ real xT=(y0 == 1) ? infinity : x0+m/(1-y0);
+ if(abs(xT-1) < eps) xT=1;
+ if(xT >= eps && xT <= 1) points.push((xT,1));
+
+ real yL=(x0 == 0) ? infinity : y0-m/x0;
+ if(abs(yL-1) < eps) yL=1;
- for(int cnt=0; cnt < c.length; ++cnt) {
- pair[][] gdscnt=points[cnt];
- for(int i=0; i < gdscnt.length; ++i) {
- pair[] gig=gdscnt[i];
- int Li=gig.length;
- for(int j=i+1; j < gdscnt.length; ++j) {
- pair[] gjg=gdscnt[j];
- int Lj=gjg.length;
- if(abs(gig[0]-gjg[0]) < eps) {
- gdscnt[j]=gjg[reverseF(Lj)];
- gdscnt[j].append(gig);
- gdscnt.delete(i);
- --i;
- break;
- } else if(abs(gig[0]-gjg[Lj-1]) < eps) {
- gig.delete(0);
- gdscnt[j].append(gig);
- gdscnt.delete(i);
- --i;
- break;
- } else if(abs(gig[Li-1]-gjg[0]) < eps) {
- gjg.delete(0);
- gig.append(gjg);
- gdscnt[j]=gig;
- gdscnt.delete(i);
- --i;
- break;
- } else if(abs(gig[Li-1]-gjg[Lj-1]) < eps) {
- gig.append(gjg[reverseL(Lj)]);
- gdscnt[j]=gig;
- gdscnt.delete(i);
- --i;
- break;
- }
+ if(yL > eps && yL <= 1) points.push((0,yL));
+
+ real yR=(x0 == 1) ? infinity : y0+m/(1-x0);
+ if(abs(yR) < eps) yR=0;
+ if(yR >= 0 && yR <= 1-eps) points.push((1,yR));
+
+ // Check (if logic is correct above, this will not happen)
+ if(!(points.length == 2 || points.length == 4)) {
+ abort("Unexpected error in setcontour routine: odd number of"
+ +" crossings (hyperbolic case)");
+ }
+
+ // Lower left side
+ if((x0 > 0) && (y0 > 0) && (f00*signc < 0)) {
+ pair[] pts0;
+ for(int i=0; i < points.length; ++i) {
+ if((points[i].x < x0) && (points[i].y < y0)) {
+ pts0.push(points[i]);
+ }
+ }
+ if(pts0.length == 2) {
+ pair a0,b0;
+ if((f00 > 0) ^(pts0[0].x < pts0[1].x)) {
+ a0=pts0[0];
+ b0=pts0[1];
+ } else {
+ a0=pts0[1];
+ b0=pts0[0];
+ }
+ boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,-1,-1,i,j,index);
+ bd.data.push(bc);
+ }
+ }
+
+ // Lower right side
+ if((x0 < 1) && (y0 > 0) && (f10*signc < 0)) {
+ pair[] pts0;
+ for(int i=0; i < points.length; ++i) {
+ if((points[i].x > x0) && (points[i].y < y0)) {
+ pts0.push(points[i]);
+ }
+ }
+ if(pts0.length == 2) {
+ pair a0,b0;
+ if((f10 > 0) ^(pts0[0].x < pts0[1].x)) {
+ a0=pts0[0];
+ b0=pts0[1];
+ } else {
+ a0=pts0[1];
+ b0=pts0[0];
+ }
+ boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,+1,-1,i,j,index);
+ bd.data.push(bc);
+ }
+ }
+
+ // Upper right side
+ if((x0 < 1) && (y0 < 1) && (f11*signc < 0)) {
+ pair[] pts0;
+ for(int i=0; i < points.length; ++i) {
+ if((points[i].x > x0) && (points[i].y > y0)) {
+ pts0.push(points[i]);
+ }
+ }
+ if(pts0.length == 2) {
+ pair a0,b0;
+ if((f11 > 0) ^(pts0[0].x > pts0[1].x)) {
+ a0=pts0[0];
+ b0=pts0[1];
+ } else {
+ a0=pts0[1];
+ b0=pts0[0];
}
+ boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,+1,+1,i,j,index);
+ bd.data.push(bc);
}
}
+
+ // Upper left side
+ if((x0 > 0) && (y0 < 1) && (f01*signc < 0)) {
+ pair[] pts0;
+ for(int i=0; i < points.length; ++i) {
+ if((points[i].x < x0) && (points[i].y > y0)) {
+ pts0.push(points[i]);
+ }
+ }
+ if(pts0.length == 2) {
+ pair a0,b0;
+ if((f01 > 0) ^(pts0[0].x > pts0[1].x)) {
+ a0=pts0[0];
+ b0=pts0[1];
+ } else {
+ a0=pts0[1];
+ b0=pts0[0];
+ }
+ boxcontour bc=boxcontour(hyperbola,a0,b0,x0,y0,m,-1,+1,i,j,index);
+ bd.data.push(bc);
+ }
+ }
+ return;
}
-// Join path segments.
-private guide[][] connect(pair[][][] points, real[] c, interpolate join)
-{
- // set up return value
- guide[][] result=new guide[c.length][];
- for(int cnt=0; cnt < c.length; ++cnt) {
- pair[][] pointscnt=points[cnt];
- guide[] resultcnt=result[cnt]=new guide[pointscnt.length];
- for(int i=0; i < pointscnt.length; ++i) {
- pair[] pts=pointscnt[i];
- guide gd;
- if(pts.length > 0) {
- if(pts.length > 1 && abs(pts[0]-pts[pts.length-1]) < eps) {
- guide[] g=sequence(new guide(int i) {
- return pts[i];
- },pts.length-1);
- g.push(cycle);
- gd=join(...g);
- } else
- gd=join(...sequence(new guide(int i) {
- return pts[i];
- },pts.length));
+
+// Checks if end of first contour segment matches the beginning of
+// the second.
+private bool connected(boxcontour bc1, boxcontour bc2) {
+ return abs(bc2.a-bc1.b+(bc2.i-bc1.i,bc2.j-bc1.j)) < eps;
+}
+
+// Returns index of first active element in bca that with beginning
+// that connects to the end of bc, or -1 if no such element.
+private int connectedindex(boxcontour bc, boxcontour[] bca,
+ bool activeonly=true) {
+ for(int i=0; i < bca.length; ++i) {
+ if(!bca[i].active) continue;
+ if(connected(bc,bca[i])) {
+ return i;
+ }
+ }
+ return -1;
+}
+
+// Returns index of first active element in bca with end that connects
+// to the start of bc, or -1 if no such element.
+private int connectedindex(boxcontour[] bca, boxcontour bc,
+ bool activeonly=true) {
+ for(int i=0; i < bca.length; ++i) {
+ if(!bca[i].active) continue;
+ if(connected(bca[i],bc)) {
+ return i;
+ }
+ }
+ return -1;
+}
+
+
+// Processes indices for grid regions touching the
+// end/start (forward=true/false) of the contour segment
+private void searchindex(boxcontour bc, bool forward, void f(int i, int j)) {
+ pair z=forward ? bc.b : bc.a;
+
+ int i=bc.i;
+ int j=bc.j;
+
+ if(z == (0,0)) f(i-1,j-1);
+ if(z.y == 0) f(i,j-1);
+ if(z == (1,0)) f(i+1,j-1);
+ if(z.x == 1) f(i+1,j);
+ if(z == (1,1)) f(i+1,j+1);
+ if(z.y == 1) f(i,j+1);
+ if(z == (0,1)) f(i-1,j+1);
+ if(z.x == 0) f(i-1,j);
+}
+
+// Contour segment
+private struct segment {
+ gridpoint[] data;
+ void operator init() {
+ }
+ void operator init(boxcontour bc, int subsample=1) {
+ bc.active=false;
+ this.data.append(bc.points(subsample,first=true,last=true));
+ }
+ void operator init(int i, int j, pair z) {
+ gridpoint gp=gridpoint(i,j,z);
+ data.push(gp);
+ }
+ void operator init(gridpoint[] gp) {
+ this.data.append(gp);
+ }
+ gridpoint start() {
+ if(data.length == 0) {
+ return gridpoint(-1,-1,(-infinity,-infinity));
+ }
+ gridpoint gp=data[0];
+ return gridpoint(gp.i,gp.j,gp.z);
+ }
+ gridpoint end() {
+ if(data.length == 0) {
+ return gridpoint(-1,-1,(-infinity,-infinity));
+ }
+ gridpoint gp=data[data.length-1];
+ return gridpoint(gp.i,gp.j,gp.z);
+ }
+ bool closed() {
+ return same(this.start(),this.end());
+ }
+ void append(boxcontour bc, int subsample=1) {
+ bc.active=false;
+ data.append(bc.points(subsample,first=false,last=true));
+ }
+ void prepend(boxcontour bc, int subsample=1) {
+ bc.active=false;
+ data.insert(0 ... bc.points(subsample,first=true,last=false));
+ }
+ void append(int i, int j, pair z) {
+ gridpoint gp=gridpoint(i,j,z);
+ data.push(gp);
+ }
+ void prepend(int i, int j, pair z) {
+ gridpoint gp=gridpoint(i,j,z);
+ data.insert(0,gp);
+ }
+ segment copy() {
+ segment seg=new segment;
+ seg.data=new gridpoint[data.length];
+ for(int i=0; i < data.length; ++i) {
+ seg.data[i]=gridpoint(data[i].i,data[i].j,data[i].z);
+ }
+ return seg;
+ }
+ segment reversecopy() {
+ segment seg=new segment;
+ seg.data=new gridpoint[data.length];
+ for(int i=0; i < data.length; ++i) {
+ seg.data[data.length-i-1]=gridpoint(data[i].i,data[i].j,data[i].z);
+ }
+ return seg;
+ }
+}
+
+// Container to hold edge and border segments that form one continuous line
+private struct Segment {
+ segment[] edges;
+ segment[] borders;
+ void operator init() {
+ }
+ void operator init(segment seg) {
+ edges.push(seg);
+ }
+ void operator init(gridpoint[] gp) {
+ segment seg=segment(gp);
+ edges.push(seg);
+ }
+ gridpoint start() {
+ if(edges.length == 0) {
+ if(borders.length > 0) {
+ return borders[0].start();
}
- resultcnt[i]=gd;
+ return gridpoint(-1,-1,(-infinity,-infinity));
}
+ return edges[0].start();
+ }
+ gridpoint end() {
+ if(edges.length == 0 && borders.length == 0) {
+ return gridpoint(-1,-1,(-infinity,-infinity));
+ }
+ if(edges.length > borders.length) {
+ return edges[edges.length-1].end();
+ } else {
+ return borders[borders.length-1].end();
+ }
+ }
+ bool closed() {
+ return same(this.start(),this.end());
+ }
+ void addedge(segment seg) {
+ edges.push(seg);
+ }
+ void addedge(gridpoint[] gp) {
+ segment seg=segment(gp);
+ edges.push(seg);
+ }
+ void addborder(segment seg) {
+ borders.push(seg);
+ }
+ void addborder(gridpoint[] gp) {
+ segment seg=segment(gp);
+ borders.push(seg);
+ }
+ void append(Segment S) {
+ edges.append(S.edges);
+ borders.append(S.borders);
}
- return result;
}
+private Segment[] Segment(segment[] s)
+{
+ return sequence(new Segment(int i) {return Segment(s[i]);},s.length);
+}
-// Return contour guides for a 2D data array.
-// z: two-dimensional array of nonoverlapping mesh points
+private Segment[][] Segment(segment[][] s)
+{
+ Segment[][] S=new Segment[s.length][];
+ for(int i=0; i < s.length; ++i)
+ S[i]=Segment(s[i]);
+ return S;
+}
+
+// Return contour points for a 2D data array.
// f: two-dimensional array of corresponding f(z) data values
-// midpoint: optional array containing values of f at cell midpoints
// c: array of contour values
-// join: interpolation operator (e.g. operator -- or operator ..)
-guide[][] contour(pair[][] z, real[][] f,
- real[][] midpoint=new real[][], real[] c,
- interpolate join=operator --)
+// subsample: number of points to use in each box in addition to endpoints
+segment[][] contouredges(real[][] f, real[] c, int subsample=1)
{
- int nx=z.length-1;
- if(nx == 0)
- abort("array z must have length >= 2");
- int ny=z[0].length-1;
- if(ny == 0)
- abort("array z[0] must have length >= 2");
+ int nx=f.length-1;
+ if(nx <= 0)
+ abort("array f must have length >= 2");
+ int ny=f[0].length-1;
+ if(ny <= 0)
+ abort("array f[0] must have length >= 2");
c=sort(c);
- bool midpoints=midpoint.length > 0;
+ boxdata[][] bd=new boxdata[nx][ny];
+
+ segment[][] result=new segment[c.length][];
- segment segments[][][]=new segment[nx][ny][];
-
- // go over region a rectangle at a time
for(int i=0; i < nx; ++i) {
- pair[] zi=z[i];
- pair[] zp=z[i+1];
+ boxdata[] bdi=bd[i];
real[] fi=f[i];
real[] fp=f[i+1];
- real[] midpointi;
- if(midpoints) midpointi=midpoint[i];
- segment[][] segmentsi=segments[i];
+
for(int j=0; j < ny; ++j) {
- segment[] segmentsij=segmentsi[j];
-
- // define points
- pair bleft=zi[j];
- pair bright=zp[j];
- pair tleft=zi[j+1];
- pair tright=zp[j+1];
- pair middle=0.25*(bleft+bright+tleft+tright);
-
- real f00=fi[j];
- real f01=fi[j+1];
- real f10=fp[j];
- real f11=fp[j+1];
- real fmm=midpoints ? midpoint[i][j] : 0.25*(f00+f01+f10+f11);
-
- // optimization: we make sure we don't work with empty rectangles
+ boxdata bdij=bdi[j]=boxdata(fi[j],fp[j],fi[j+1],fp[j+1]);
+
int checkcell(int cnt) {
real C=c[cnt];
- real vertdat0=f00-C; // bottom-left vertex
- real vertdat1=f10-C; // bottom-right vertex
- real vertdat2=f01-C; // top-left vertex
- real vertdat3=f11-C; // top-right vertex
+
+ real f00=fi[j];
+ real f10=fp[j];
+ real f01=fi[j+1];
+ real f11=fp[j+1];
+
+ real epsf=eps*max(abs(f00),abs(f10),abs(f01),abs(f11),abs(C));
+
+ f00=f00-C;
+ f10=f10-C;
+ f01=f01-C;
+ f11=f11-C;
+
+ if(abs(f00) < epsf) f00=0;
+ if(abs(f10) < epsf) f10=0;
+ if(abs(f01) < epsf) f01=0;
+ if(abs(f11) < epsf) f11=0;
+
- // optimization: we make sure we don't work with empty rectangles
int countm=0;
int countz=0;
int countp=0;
-
+
void check(real vertdat) {
- if(vertdat < -eps) ++countm;
+ if(vertdat < -eps)++countm;
else {
- if(vertdat <= eps) ++countz;
- else ++countp;
+ if(vertdat <= eps)++countz;
+ else++countp;
}
}
- check(vertdat0);
- check(vertdat1);
- check(vertdat2);
- check(vertdat3);
+ check(f00);
+ check(f10);
+ check(f01);
+ check(f11);
if(countm == 4) return 1; // nothing to do
if(countp == 4) return -1; // nothing to do
if((countm == 3 || countp == 3) && countz == 1) return 0;
- // go through the triangles
-
- void addseg(segment seg) {
- if(seg.active) {
- seg.c=cnt;
- segmentsij.push(seg);
- }
- }
- real vertdat4=fmm-C;
- addseg(checktriangle(bright,tright,middle,
- vertdat1,vertdat3,vertdat4,0));
- addseg(checktriangle(tright,tleft,middle,
- vertdat3,vertdat2,vertdat4,1));
- addseg(checktriangle(tleft,bleft,middle,
- vertdat2,vertdat0,vertdat4,2));
- addseg(checktriangle(bleft,bright,middle,
- vertdat0,vertdat1,vertdat4,3));
+ // Calculate individual box contours
+ bdij.settype(C);
+ if(bdij.type == edge)
+ setcontour(f00,f10,f01,f11,epsf,bdij,i,j,cnt);
return 0;
}
-
+
void process(int l, int u) {
if(l >= u) return;
int i=quotient(l+u,2);
@@ -313,150 +929,171 @@ guide[][] contour(pair[][] z, real[][] f,
process(i+1,u);
}
}
-
+
process(0,c.length);
}
}
-
- // set up return value
- pair[][][] points=new pair[c.length][][];
-
+
+ // Find contours and follow them
for(int i=0; i < nx; ++i) {
- segment[][] segmentsi=segments[i];
+ boxdata[] bdi=bd[i];
for(int j=0; j < ny; ++j) {
- segment[] segmentsij=segmentsi[j];
- for(int k=0; k < segmentsij.length; ++k) {
- segment C=segmentsij[k];
-
- if(!C.active) continue;
-
- pair[] g=new pair[] {C.a,C.b};
- segmentsij[k].active=false;
-
- int forward(int I, int J, bool first=true) {
- if(I >= 0 && I < nx && J >= 0 && J < ny) {
- segment[] segmentsIJ=segments[I][J];
- for(int l=0; l < segmentsIJ.length; ++l) {
- segment D=segmentsIJ[l];
- if(!D.active) continue;
- if(abs(D.a-g[g.length-1]) < eps) {
- g.push(D.b);
- segmentsIJ[l].active=false;
- if(D.edge >= 0 && !first) return D.edge;
- first=false;
- l=-1;
- } else if(abs(D.b-g[g.length-1]) < eps) {
- g.push(D.a);
- segmentsIJ[l].active=false;
- if(D.edge >= 0 && !first) return D.edge;
- first=false;
- l=-1;
- }
- }
- }
- return -1;
- }
-
- int backward(int I, int J, bool first=true) {
- if(I >= 0 && I < nx && J >= 0 && J < ny) {
- segment[] segmentsIJ=segments[I][J];
- for(int l=0; l < segmentsIJ.length; ++l) {
- segment D=segmentsIJ[l];
- if(!D.active) continue;
- if(abs(D.a-g[0]) < eps) {
- g.insert(0,D.b);
- segmentsIJ[l].active=false;
- if(D.edge >= 0 && !first) return D.edge;
- first=false;
- l=-1;
- } else if(abs(D.b-g[0]) < eps) {
- g.insert(0,D.a);
- segmentsIJ[l].active=false;
- if(D.edge >= 0 && !first) return D.edge;
- first=false;
- l=-1;
+ boxdata bd0=bdi[j];
+ if(bd0.count() == 0) continue;
+ for(int k=0; k < bd0.count(); ++k) {
+ boxcontour bc0=bd0.data[k];
+
+ if(!bc0.active) continue;
+
+ // Note: boxcontour set inactive when added to segment
+ segment seg=segment(bc0,subsample);
+
+ // Forward direction
+ bool foundnext=true;
+ while(foundnext) {
+ foundnext=false;
+ searchindex(bc0,true,new void(int i, int j) {
+ if((i >= 0) && (i < nx) && (j >= 0) && (j < ny)) {
+ boxcontour[] data=bd[i][j].data;
+ int k0=connectedindex(bc0,data);
+ if(k0 >= 0) {
+ bc0=data[k0];
+ seg.append(bc0,subsample);
+ foundnext=true;
+ }
}
- }
- }
- return -1;
+ });
}
-
- void follow(int f(int, int, bool first=true), int edge) {
- int I=i;
- int J=j;
- while(true) {
- static int ix[]={1,0,-1,0};
- static int iy[]={0,1,0,-1};
- if(edge >= 0 && edge < 4) {
- I += ix[edge];
- J += iy[edge];
- edge=f(I,J);
- } else {
- if(edge == -1) break;
- if(edge < 9) {
- int edge0=(edge-5) % 4;
- int edge1=(edge-4) % 4;
- int ix0=ix[edge0];
- int iy0=iy[edge0];
- I += ix0;
- J += iy0;
- // Search all 3 corner cells
- if((edge=f(I,J)) == -1) {
- I += ix[edge1];
- J += iy[edge1];
- if((edge=f(I,J)) == -1) {
- I -= ix0;
- J -= iy0;
- edge=f(I,J);
- }
- }
- } else {
- // Double-vertex edge: search all 8 surrounding cells
- void search() {
- for(int i=-1; i <= 1; ++i) {
- for(int j=-1; j <= 1; ++j) {
- if((edge=f(I+i,J+j,false)) >= 0) {
- I += i;
- J += j;
- return;
- }
- }
- }
+
+ // Backward direction
+ bc0=bd0.data[k];
+ bool foundprev=true;
+ while(foundprev) {
+ foundprev=false;
+ searchindex(bc0,false,new void(int i, int j) {
+ if((i >= 0) && (i < nx) && (j >= 0) && (j < ny)) {
+ boxcontour[] data=bd[i][j].data;
+ int k0=connectedindex(data,bc0);
+ if(k0 >= 0) {
+ bc0=data[k0];
+ seg.prepend(bc0,subsample);
+ foundprev=true;
}
- search();
}
- }
- }
+ });
}
- // Follow contour in cell
- int edge=forward(i,j,first=false);
-
- // Follow contour forward outside of cell
- follow(forward,edge);
-
- // Follow contour backward outside of cell
- follow(backward,C.edge);
-
- points[C.c].push(g);
+ result[bc0.index].push(seg);
}
}
}
+
+ // Note: every segment here _should_ be cyclic or terminate on the
+ // boundary
+ return result;
+}
+
+// Connect contours into guides.
+// Same initial/final points indicates a closed path.
+// Borders are always joined using--.
+private guide connect(Segment S, pair[][] z, interpolate join)
+{
+ pair loc(gridpoint gp) {
+ pair offset=z[gp.i][gp.j];
+ pair size=z[gp.i+1][gp.j+1]-z[gp.i][gp.j];
+ return offset+(size.x*gp.z.x,size.y*gp.z.y);
+ }
+ pair[] loc(gridpoint[] gp) {
+ pair[] result=new pair[gp.length];
+ for(int i; i < gp.length; ++i) {
+ result[i]=loc(gp[i]);
+ }
+ return result;
+ }
+
+ bool closed=S.closed();
+
+ pair[][] edges=new pair[S.edges.length][];
+ for(int i; i < S.edges.length; ++i) {
+ edges[i]=loc(S.edges[i].data);
+ }
+ pair[][] borders=new pair[S.borders.length][];
+ for(int i; i < S.borders.length; ++i) {
+ borders[i]=loc(S.borders[i].data);
+ }
+
+ if(edges.length == 0 && borders.length == 1) {
+ guide g=operator--(...borders[0]);
+ if(closed) g=g--cycle;
+ return g;
+ }
+
+ if(edges.length == 1 && borders.length == 0) {
+ pair[] pts=edges[0];
+ if(closed) pts.delete(pts.length-1);
+ guide g=join(...pts);
+ if(closed) g=join(g,cycle);
+ return g;
+ }
+
+ guide[] ge=new guide[edges.length];
+ for(int i=0; i < ge.length; ++i)
+ ge[i]=join(...edges[i]);
+
+ guide[] gb=new guide[borders.length];
+ for(int i=0; i < gb.length; ++i)
+ gb[i]=operator--(...borders[i]);
+
+ guide g=ge[0];
+ if(0 < gb.length) g=g&gb[0];
+ for(int i=1; i < ge.length; ++i) {
+ g=g&ge[i];
+ if(i < gb.length) g=g&gb[i];
+ }
+ if(closed) g=g&cycle;
+ return g;
+}
- collect(points,c); // Required to join remaining case1 cycles.
+// Connect contours into guides.
+private guide[] connect(Segment[] S, pair[][] z, interpolate join)
+{
+ return sequence(new guide(int i) {return connect(S[i],z,join);},S.length);
+}
- return connect(points,c,join);
+// Connect contours into guides.
+private guide[][] connect(Segment[][] S, pair[][] z, interpolate join)
+{
+ guide[][] result=new guide[S.length][];
+ for(int i=0; i < S.length; ++i) {
+ result[i]=connect(S[i],z,join);
+ }
+ return result;
+}
+
+// Return contour guides for a 2D data array.
+// z: two-dimensional array of nonoverlapping mesh points
+// f: two-dimensional array of corresponding f(z) data values
+// c: array of contour values
+// join: interpolation operator (e.g. operator--or operator ..)
+// subsample: number of interior points to include in each grid square
+// (in addition to points on edge)
+guide[][] contour(pair[][] z, real[][] f, real[] c,
+ interpolate join=operator--, int subsample=1)
+{
+ segment[][] seg=contouredges(f,c,subsample);
+ Segment[][] Seg=Segment(seg);
+ return connect(Seg,z,join);
}
// Return contour guides for a 2D data array on a uniform lattice
// f: two-dimensional array of real data values
-// midpoint: optional array containing data values at cell midpoints
// a,b: diagonally opposite vertices of rectangular domain
// c: array of contour values
-// join: interpolation operator (e.g. operator -- or operator ..)
-guide[][] contour(real[][] f, real[][] midpoint=new real[][],
- pair a, pair b, real[] c,
- interpolate join=operator --)
+// join: interpolation operator (e.g. operator--or operator ..)
+// subsample: number of interior points to include in each grid square
+// (in addition to points on edge)
+guide[][] contour(real[][] f, pair a, pair b, real[] c,
+ interpolate join=operator--, int subsample=1)
{
int nx=f.length-1;
if(nx == 0)
@@ -473,37 +1110,43 @@ guide[][] contour(real[][] f, real[][] midpoint=new real[][],
zi[j]=(xi,interp(a.y,b.y,j/ny));
}
}
- return contour(z,f,midpoint,c,join);
+ return contour(z,f,c,join,subsample);
}
// return contour guides for a real-valued function
-// f: real-valued function of two real variables
-// a,b: diagonally opposite vertices of rectangular domain
-// c: array of contour values
-// nx,ny: number of subdivisions in x and y directions (determines accuracy)
-// join: interpolation operator (e.g. operator -- or operator ..)
+// f: real-valued function of two real variables
+// a,b: diagonally opposite vertices of rectangular domain
+// c: array of contour values
+// nx,ny: number of subdivisions in x and y directions(determines accuracy)
+// join: interpolation operator (e.g. operator--or operator ..)
+// subsample: number of interior points to include in each grid square
+// (in addition to points on edge)
guide[][] contour(real f(real, real), pair a, pair b,
real[] c, int nx=ngraph, int ny=nx,
- interpolate join=operator --)
+ interpolate join=operator--, int subsample=1)
{
- // evaluate function at points and midpoints
+ // evaluate function at points and subsample
real[][] dat=new real[nx+1][ny+1];
- real[][] midpoint=new real[nx+1][ny+1];
for(int i=0; i <= nx; ++i) {
real x=interp(a.x,b.x,i/nx);
- real x2=interp(a.x,b.x,(i+0.5)/nx);
real[] dati=dat[i];
- real[] midpointi=midpoint[i];
for(int j=0; j <= ny; ++j) {
dati[j]=f(x,interp(a.y,b.y,j/ny));
- midpointi[j]=f(x2,interp(a.y,b.y,(j+0.5)/ny));
}
}
- return contour(dat,midpoint,a,b,c,join);
+ return contour(dat,a,b,c,join,subsample);
}
-
+
+guide[][] contour(real f(pair), pair a, pair b,
+ real[] c, int nx=ngraph, int ny=nx,
+ interpolate join=operator--, int subsample=1)
+{
+ return contour(new real(real x, real y) {return f((x,y));},
+ a,b,c,nx,ny,join,subsample);
+}
+
void draw(picture pic=currentpicture, Label[] L=new Label[],
guide[][] g, pen[] p)
{
@@ -530,6 +1173,239 @@ void draw(picture pic=currentpicture, Label[] L=new Label[],
draw(pic,L,g,sequence(new pen(int) {return p;},g.length));
}
+// Draw the contour
+void draw(picture pic=currentpicture, Label L,
+ guide[] g, pen p=currentpen)
+{
+ draw(pic,g,p);
+ for(int i=0; i < g.length; ++i) {
+ if(L.s != "" && size(g[i]) > 1) {
+ label(pic,L,g[i],p);
+ }
+ }
+}
+
+/* CONTOURS FOR IRREGULARLY SPACED POINTS */
+//
+// +---------+
+// |\ /|
+// | \ / |
+// | \ / |
+// | \ / |
+// | X |
+// | / \ |
+// | / \ |
+// | / \ |
+// |/ \|
+// +---------+
+//
+
+// Is triangle p0--p1--p2--cycle counterclockwise ?
+private bool isCCW(pair p0, pair p1, pair p2) {return side(p0,p1,p2) < 0;}
+
+private struct segment
+{
+ bool active;
+ bool reversed; // True if lower values are to the left along line a--b.
+ pair a,b; // Endpoints; a is always an edge point if one exists.
+ int c; // Contour value.
+}
+
+// Case 1: line passes through two vertices of a triangle
+private segment case1(pair p0, pair p1, pair p2,
+ real v0, real v1, real v2)
+{
+ // Will cause a duplicate guide; luckily case1 is rare
+ segment rtrn;
+ rtrn.active=true;
+ rtrn.a=p0;
+ rtrn.b=p1;
+ rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0));
+ return rtrn;
+}
+
+// Cases 2 and 3: line passes through a vertex and a side of a triangle
+//(the first vertex passed and the side between the other two)
+private segment case2(pair p0, pair p1, pair p2,
+ real v0, real v1, real v2)
+{
+ segment rtrn;
+ rtrn.active=true;
+ pair val=interp(p1,p2,abs(v1/(v2-v1)));
+ rtrn.a=val;
+ rtrn.b=p0;
+ rtrn.reversed=!(isCCW(p0,p1,p2) ^(v2 > 0));
+ return rtrn;
+}
+
+private segment case3(pair p0, pair p1, pair p2,
+ real v0, real v1, real v2)
+{
+ segment rtrn;
+ rtrn.active=true;
+ pair val=interp(p1,p2,abs(v1/(v2-v1)));
+ rtrn.a=p0;
+ rtrn.b=val;
+ rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0));
+ return rtrn;
+}
+
+// Case 4: line passes through two sides of a triangle
+//(through the sides formed by the first&second, and second&third vertices)
+private segment case4(pair p0, pair p1, pair p2,
+ real v0, real v1, real v2)
+{
+ segment rtrn;
+ rtrn.active=true;
+ rtrn.a=interp(p1,p0,abs(v1/(v0-v1)));
+ rtrn.b=interp(p1,p2,abs(v1/(v2-v1)));
+ rtrn.reversed=(isCCW(p0,p1,p2) ^(v2 > 0));
+ return rtrn;
+}
+
+// Check if a line passes through a triangle, and draw the required line.
+private segment checktriangle(pair p0, pair p1, pair p2,
+ real v0, real v1, real v2)
+{
+ // default null return
+ static segment dflt;
+
+ real eps=eps*max(abs(v0),abs(v1),abs(v2),1);
+
+ if(v0 < -eps) {
+ if(v1 < -eps) {
+ if(v2 < -eps) return dflt; // nothing to do
+ else if(v2 <= eps) return dflt; // nothing to do
+ else return case4(p0,p2,p1,v0,v2,v1);
+ } else if(v1 <= eps) {
+ if(v2 < -eps) return dflt; // nothing to do
+ else if(v2 <= eps) return case1(p1,p2,p0,v1,v2,v0);
+ else return case3(p1,p0,p2,v1,v0,v2);
+ } else {
+ if(v2 < -eps) return case4(p0,p1,p2,v0,v1,v2);
+ else if(v2 <= eps)
+ return case2(p2,p0,p1,v2,v0,v1);
+ else return case4(p1,p0,p2,v1,v0,v2);
+ }
+ } else if(v0 <= eps) {
+ if(v1 < -eps) {
+ if(v2 < -eps) return dflt; // nothing to do
+ else if(v2 <= eps) return case1(p0,p2,p1,v0,v2,v1);
+ else return case2(p0,p1,p2,v0,v1,v2);
+ } else if(v1 <= eps) {
+ if(v2 < -eps) return case1(p0,p1,p2,v0,v1,v2);
+ else if(v2 <= eps) return dflt; // use finer partitioning.
+ else return case1(p0,p1,p2,v0,v1,v2);
+ } else {
+ if(v2 < -eps) return case2(p0,p1,p2,v0,v1,v2);
+ else if(v2 <= eps) return case1(p0,p2,p1,v0,v2,v1);
+ else return dflt; // nothing to do
+ }
+ } else {
+ if(v1 < -eps) {
+ if(v2 < -eps) return case4(p1,p0,p2,v1,v0,v2);
+ else if(v2 <= eps)
+ return case2(p2,p0,p1,v2,v0,v1);
+ else return case4(p0,p1,p2,v0,v1,v2);
+ } else if(v1 <= eps) {
+ if(v2 < -eps) return case3(p1,p0,p2,v1,v0,v2);
+ else if(v2 <= eps) return case1(p1,p2,p0,v1,v2,v0);
+ else return dflt; // nothing to do
+ } else {
+ if(v2 < -eps) return case4(p0,p2,p1,v0,v2,v1);
+ else if(v2 <= eps) return dflt; // nothing to do
+ else return dflt; // nothing to do
+ }
+ }
+}
+
+// Collect connecting path segments.
+private void collect(pair[][][] points, real[] c)
+{
+ for(int cnt=0; cnt < c.length; ++cnt) {
+ pair[][] gdscnt=points[cnt];
+ for(int i=0; i < gdscnt.length; ++i) {
+ pair[] gig=gdscnt[i];
+ int Li=gig.length;
+ for(int j=i+1; j < gdscnt.length; ++j) {
+ pair[] gjg=gdscnt[j];
+ int Lj=gjg.length;
+ if(abs(gig[0]-gjg[Lj-1]) < eps) {
+ gig.delete(0);
+ gdscnt[j].append(gig);
+ gdscnt.delete(i);
+ --i;
+ break;
+ } else if(abs(gig[Li-1]-gjg[0]) < eps) {
+ gjg.delete(0);
+ gig.append(gjg);
+ gdscnt[j]=gig;
+ gdscnt.delete(i);
+ --i;
+ break;
+ }
+ }
+ }
+ }
+}
+
+// Join path segments.
+private guide[][] connect(pair[][][] points, real[] c, interpolate join)
+{
+ // set up return value
+ guide[][] result=new guide[c.length][];
+ for(int cnt=0; cnt < c.length; ++cnt) {
+ pair[][] pointscnt=points[cnt];
+ guide[] resultcnt=result[cnt]=new guide[pointscnt.length];
+ for(int i=0; i < pointscnt.length; ++i) {
+ pair[] pts=pointscnt[i];
+ guide gd;
+ if(pts.length > 0) {
+ if(pts.length > 1 && abs(pts[0]-pts[pts.length-1]) < eps) {
+ guide[] g=sequence(new guide(int i) {
+ return pts[i];
+ },pts.length-1);
+ g.push(cycle);
+ gd=join(...g);
+ } else
+ gd=join(...sequence(new guide(int i) {
+ return pts[i];
+ },pts.length));
+ }
+ resultcnt[i]=gd;
+ }
+ }
+ return result;
+}
+
+guide[][] contour(pair[] z, real[] f, real[] c, interpolate join=operator--)
+{
+ if(z.length != f.length)
+ abort("z and f arrays have different lengths");
+
+ int[][] trn=triangulate(z);
+
+ // array to store guides found so far
+ pair[][][] points=new pair[c.length][][];
+
+ for(int cnt=0; cnt < c.length; ++cnt) {
+ pair[][] pointscnt=points[cnt];
+ real C=c[cnt];
+ for(int i=0; i < trn.length; ++i) {
+ int[] trni=trn[i];
+ int i0=trni[0], i1=trni[1], i2=trni[2];
+ segment seg=checktriangle(z[i0],z[i1],z[i2],f[i0]-C,f[i1]-C,f[i2]-C);
+ if(seg.active)
+ pointscnt.push(seg.reversed ? new pair[] {seg.b,seg.a} :
+ new pair[] {seg.a,seg.b});
+ }
+ }
+
+ collect(points,c);
+
+ return connect(points,c,join);
+}
+
// Extend palette by the colors below and above at each end.
pen[] extend(pen[] palette, pen below, pen above) {
pen[] p=copy(palette);
@@ -614,69 +1490,3 @@ void fill(picture pic=currentpicture, guide[][] g, pen[][] palette)
}
}
}
-
-// routines for irregularly spaced points:
-
-// check existing guides and adds new segment to them if possible,
-// or otherwise store segment as a new guide
-private void addseg(pair[][] gds, segment seg)
-{
- if(!seg.active) return;
- // search for a path to extend
- for(int i=0; i < gds.length; ++i) {
- pair[] gd=gds[i];
- if(abs(gd[0]-seg.b) < eps) {
- gd.insert(0,seg.a);
- return;
- } else if(abs(gd[gd.length-1]-seg.b) < eps) {
- gd.push(seg.a);
- return;
- } else if(abs(gd[0]-seg.a) < eps) {
- gd.insert(0,seg.b);
- return;
- } else if(abs(gd[gd.length-1]-seg.a) < eps) {
- gd.push(seg.b);
- return;
- }
- }
-
- // in case nothing is found
- pair[] segm;
- segm=new pair[] {seg.a,seg.b};
- gds.push(segm);
-
- return;
-}
-
-guide[][] contour(real f(pair), pair a, pair b,
- real[] c, int nx=ngraph, int ny=nx,
- interpolate join=operator --)
-{
- return contour(new real(real x, real y) {return f((x,y));},a,b,c,nx,ny,join);
-}
-
-guide[][] contour(pair[] z, real[] f, real[] c, interpolate join=operator --)
-{
- if(z.length != f.length)
- abort("z and f arrays have different lengths");
-
- int[][] trn=triangulate(z);
-
- // array to store guides found so far
- pair[][][] points=new pair[c.length][][];
-
- for(int cnt=0; cnt < c.length; ++cnt) {
- pair[][] pointscnt=points[cnt];
- real C=c[cnt];
- for(int i=0; i < trn.length; ++i) {
- int[] trni=trn[i];
- int i0=trni[0], i1=trni[1], i2=trni[2];
- addseg(pointscnt,checktriangle(z[i0],z[i1],z[i2],
- f[i0]-C,f[i1]-C,f[i2]-C));
- }
- }
-
- collect(points,c);
-
- return connect(points,c,join);
-}
diff --git a/Build/source/utils/asymptote/base/contour3.asy b/Build/source/utils/asymptote/base/contour3.asy
index e9d8d928991..3d925be684e 100644
--- a/Build/source/utils/asymptote/base/contour3.asy
+++ b/Build/source/utils/asymptote/base/contour3.asy
@@ -181,7 +181,7 @@ vertex[][] contour3(triple[][][] v, real[][][] f,
addval(w.kpb0,w.kpb1,w.kpb2,val2,w.v);
}
- triple dir=P.vector();
+ triple dir=P.normal;
void addnormals(weighted[] pts) {
triple vec2=pts[1].v-pts[0].v;
diff --git a/Build/source/utils/asymptote/base/embed.asy b/Build/source/utils/asymptote/base/embed.asy
index ac8a09acb83..504b3af66b0 100644
--- a/Build/source/utils/asymptote/base/embed.asy
+++ b/Build/source/utils/asymptote/base/embed.asy
@@ -1,11 +1,12 @@
if(latex()) {
+usepackage("hyperref");
+texpreamble("\hypersetup{"+settings.hyperrefOptions+"}");
texpreamble("
\ifx\pdfhorigin\undefined%
\usepackage[3D,dvipdfmx]{movie15}
\else%
\usepackage[3D]{movie15}
\fi%");
-usepackage("hyperref",settings.hyperrefOptions);
}
// See http://www.ctan.org/tex-archive/macros/latex/contrib/movie15/README
diff --git a/Build/source/utils/asymptote/base/external.asy b/Build/source/utils/asymptote/base/external.asy
index 9560e319299..053f1eff4c1 100644
--- a/Build/source/utils/asymptote/base/external.asy
+++ b/Build/source/utils/asymptote/base/external.asy
@@ -1,4 +1,5 @@
-usepackage("hyperref",settings.hyperrefOptions);
+usepackage("hyperref");
+texpreamble("\hypersetup{"+settings.hyperrefOptions+"}");
// Embed object to be run in an external window. An image file name can be
// specified; if not given one will be automatically generated.
@@ -19,8 +20,8 @@ string embed(string name, string options="", real width=0, real height=0,
atexit(exitfunction);
}
}
- if(width != 0) options += ", width="+(string) (width*pt)+"pt";
- if(height != 0) options +=", height="+(string) (height*pt)+"pt";
+ if(width != 0) options += ", width="+(string) (width/pt)+"pt";
+ if(height != 0) options +=", height="+(string) (height/pt)+"pt";
return "\href{run:"+name+"}{"+graphic(image,options)+"}";
}
diff --git a/Build/source/utils/asymptote/base/flowchart.asy b/Build/source/utils/asymptote/base/flowchart.asy
index 5fe96ae8508..bb8fe21a4d6 100644
--- a/Build/source/utils/asymptote/base/flowchart.asy
+++ b/Build/source/utils/asymptote/base/flowchart.asy
@@ -10,28 +10,17 @@ restricted flowdir Vertical;
real minblockwidth=0;
real minblockheight=0;
real mincirclediameter=0;
+real defaultexcursion=0.1;
struct block {
// The absolute center of the block in user coordinates.
pair center;
- // The relative center of the block.
- pair f_center;
-
// The size of the block
pair size;
- // Returns the relative position along the boundary of the block.
- pair f_position(real x);
-
- pair shift(transform t=identity()) {
- return t*center-f_center;
- }
-
- // Returns the absolute position along the boundary of the block.
- pair position(real x, transform t=identity()) {
- return shift(t)+f_position(x);
- }
+ // The relative center of the block.
+ pair f_center;
// These eight variables return the appropriate location on the block
// in relative coordinates, where the lower left corner of the block is (0,0).
@@ -44,6 +33,26 @@ struct block {
pair f_bottomleft;
pair f_bottomright;
+ void operator init(pair z) {
+ center=z;
+ }
+
+ void operator init(real x, real y) {
+ center=(x,y);
+ }
+
+ pair shift(transform t=identity()) {
+ return t*center-f_center;
+ }
+
+ // Returns the relative position along the boundary of the block.
+ pair f_position(real x);
+
+ // Returns the absolute position along the boundary of the block.
+ pair position(real x, transform t=identity()) {
+ return shift(t)+f_position(x);
+ }
+
// These eight functions return the appropriate location on the block
// in absolute coordinates.
pair top(transform t=identity()) {
@@ -73,6 +82,15 @@ struct block {
// Return a frame representing the block.
frame draw(pen p=currentpen);
+
+ // Store optional label on outgoing edge.
+ Label label;
+
+ // Store rectilinear path directions.
+ pair[] dirs;
+
+ // Store optional arrow.
+ arrowbar arrow=None;
};
// Construct a rectangular block with header and body objects.
@@ -160,6 +178,47 @@ block rectangle(object body, pair center=(0,0),
return block;
}
+block parallelogram(object body, pair center=(0,0),
+ pen fillpen=invisible, pen drawpen=currentpen,
+ real dx=3, real slope=2,
+ real minwidth=minblockwidth,
+ real minheight=minblockheight)
+{
+ frame f=body.f;
+ pair m=min(f);
+ pair M=max(f);
+ pair bound=maxbound(M-m+dx*(0,2),(minwidth,minheight));
+
+ real skew=bound.y/slope;
+ real a=bound.x+skew;
+ real b=bound.y;
+
+ path shape=(0,0)--(a,0)--(a+skew,b)--(skew,b)--cycle;
+
+ block block;
+ block.draw=new frame(pen p) {
+ frame block;
+ filldraw(block,shape,fillpen,drawpen);
+ add(block,shift(-0.5*(M+m))*f,((a+skew)/2,b/2));
+ return block;
+ };
+ block.f_position=new pair(real x) {
+ return point(shape,x);
+ };
+ block.f_center=((a+skew)/2,b/2);
+ block.center=center;
+ block.size=(a+skew,b);
+ block.f_bottomleft=(0,0);
+ block.f_bottom=((a+skew)/2,0);
+ block.f_bottomright=(a,0);
+ block.f_right=(a+skew/2,b/2);
+ block.f_topright=(a+skew,b);
+ block.f_top=((a+skew)/2,b);
+ block.f_topleft=(skew,b);
+ block.f_left=(skew/2,b/2);
+ return block;
+}
+
block diamond(object body, pair center=(0,0),
pen fillpen=invisible, pen drawpen=currentpen,
real ds=5, real dw=1,
@@ -327,7 +386,7 @@ path path(pair point[] ... flowdir dir[])
path line=point[0];
pair current, prev=point[0];
for(int i=1; i < point.length; ++i) {
- if(dir[i-1] == Horizontal)
+ if(i-1 >= dir.length || dir[i-1] == Horizontal)
current=(point[i].x,point[i-1].y);
else
current=(point[i-1].x,point[i].y);
@@ -348,9 +407,120 @@ path path(pair point[] ... flowdir dir[])
void draw(picture pic=currentpicture, block block, pen p=currentpen)
{
- pic.add(new void (frame f, transform t) {
+ pic.add(new void(frame f, transform t) {
add(f,shift(block.shift(t))*block.draw(p));
},true);
pic.addBox(block.center,block.center,
-0.5*block.size+min(p),0.5*block.size+max(p));
}
+
+typedef block blockconnector(block, block);
+
+blockconnector blockconnector(picture pic, transform t, pen p=currentpen,
+ margin margin=PenMargin)
+{
+ return new block(block b1, block b2) {
+ if(b1.dirs.length == 0) {
+ if(abs(b1.center.y-b2.center.y) < sqrtEpsilon) {
+ // horizontally aligned
+ b1.dirs[0]=b1.center.x < b2.center.x ? right : left;
+ blockconnector(pic,t)(b1,b2);
+ } else if(abs(b1.center.x-b2.center.x) < sqrtEpsilon) {
+ // vertically aligned
+ b1.dirs[0]=b1.center.y < b2.center.y ? up : down;
+ blockconnector(pic,t)(b1,b2);
+ } else {
+ if(abs(b1.center.y-b2.center.y) < abs(b1.center.x-b2.center.x)) {
+ b1.dirs[0]=b1.center.x < b2.center.x ? right : left;
+ b1.dirs[1]=b1.center.y < b2.center.y ? up : down;
+ blockconnector(pic,t)(b1,b2);
+ } else {
+ b1.dirs[0]=b1.center.y < b2.center.y ? up : down;
+ b1.dirs[1]=b1.center.x < b2.center.x ? right : left;
+ blockconnector(pic,t)(b1,b2);
+ }
+ }
+ return b2;
+ }
+
+ // compute the link for given directions (and label if any)
+ pair[] dirs=copy(b1.dirs); // deep copy
+ pair current,prev;
+ pair dir=dirs[0];
+ if(dir == up) prev=b1.top(t);
+ if(dir == down) prev=b1.bottom(t);
+ if(dir == left) prev=b1.left(t);
+ if(dir == right) prev=b1.right(t);
+ path line=prev;
+ arrowbar arrow=b1.arrow;
+
+ int i;
+ for(i=1; i < dirs.length-1; ++i) {
+ if(abs(length(dirs[i-1])-1) < sqrtEpsilon)
+ current=prev+t*dirs[i-1]*defaultexcursion;
+ else
+ current=prev+t*dirs[i-1];
+
+ if(current != prev) {
+ line=line--current;
+ prev=current;
+ }
+ }
+ dir=dirs[dirs.length-1];
+ current=0;
+ if(dir == up) current=b2.bottom(t);
+ if(dir == down) current=b2.top(t);
+ if(dir == left) current=b2.right(t);
+ if(dir == right) current=b2.left(t);
+ if(abs(dirs[i-1].y) < sqrtEpsilon &&
+ abs(prev.x-current.x) > sqrtEpsilon) {
+ prev=(current.x,prev.y);
+ line=line--prev; // horizontal
+ } else if(abs(dirs[i-1].x) < sqrtEpsilon &&
+ abs(prev.y-current.y) > sqrtEpsilon) {
+ prev=(prev.x,current.y);
+ line=line--prev;
+ }
+ if(current != prev)
+ line=line--current;
+
+ draw(pic,b1.label,line,p,arrow,margin);
+
+ b1.label="";
+ b1.dirs.delete();
+ b1.arrow=None;
+ return b2;
+ };
+}
+
+struct Dir
+{
+ pair z;
+ void operator init(pair z) {this.z=z;}
+}
+
+Dir Right=Dir(right);
+Dir Left=Dir(left);
+Dir Up=Dir(up);
+Dir Down=Dir(down);
+
+// Add a label to the current link
+block operator --(block b1, Label label)
+{
+ b1.label=label;
+ return b1;
+}
+
+// Add a direction to the current link
+block operator --(block b1, Dir dir)
+{
+ b1.dirs.push(dir.z);
+ return b1;
+}
+
+// Add an arrowbar to the current link
+block operator --(block b, arrowbar arrowbar)
+{
+ b.arrow=arrowbar;
+ return b;
+}
diff --git a/Build/source/utils/asymptote/base/geometry.asy b/Build/source/utils/asymptote/base/geometry.asy
index 3d9d84563b7..8679ada4c75 100644
--- a/Build/source/utils/asymptote/base/geometry.asy
+++ b/Build/source/utils/asymptote/base/geometry.asy
@@ -1163,7 +1163,7 @@ private void Drawline(picture pic=currentpicture, Label L="",pair P, bool dirP=t
// Calculate the points and direction vector in the transformed space.
pair z=t*P;
pair q=t*Q;
- pair v=t*Q-z;
+ pair v=q-z;
// path g;
pair ptp,ptq;
real cp = dirP ? 1:0;
@@ -1171,7 +1171,7 @@ private void Drawline(picture pic=currentpicture, Label L="",pair P, bool dirP=t
// Handle horizontal and vertical lines.
if(v.x == 0) {
if(m.x <= z.x && z.x <= M.x)
- if (dot(v,(z.x,m.y)) < 0) {
+ if (dot(v,m-z) < 0) {
ptp=(z.x,z.y+cp*(m.y-z.y));
ptq=(z.x,q.y+cq*(M.y-q.y));
} else {
@@ -1179,7 +1179,7 @@ private void Drawline(picture pic=currentpicture, Label L="",pair P, bool dirP=t
ptq=(z.x,z.y+cp*(M.y-z.y));
}
} else if(v.y == 0) {
- if (dot(v,(m.x,z.y)) < 0) {
+ if (dot(v,m-z) < 0) {
ptp=(z.x+cp*(m.x-z.x),z.y);
ptq=(q.x+cq*(M.x-q.x),z.y);
} else {
@@ -1210,7 +1210,7 @@ private void Drawline(picture pic=currentpicture, Label L="",pair P, bool dirP=t
lL.out(opic,g);
}
g=pathModifier(g);
- if(linetype(p) == ""){
+ if(linetype(p).length == 0){
pair m=midpoint(g);
pen tp;
tp=dirP ? p : addpenline(p);
@@ -7016,7 +7016,7 @@ arc arc(ellipse el, explicit abscissa x1, explicit abscissa x2, bool direction=C
{/*<asyxml></code><documentation>Return the arc from 'point(c,x1)' to 'point(c,x2)' in the direction 'direction'.</documentation></function></asyxml>*/
real a=degrees(point(el,x1)-el.C);
real b=degrees(point(el,x2)-el.C);
- arc oa=arc(el,a,b,fromCenter,direction);
+ arc oa=arc(el,a-el.angle,b-el.angle,fromCenter,direction);
return oa;
}
@@ -7135,8 +7135,7 @@ void perpendicular(picture pic=currentpicture, pair z, pair align, path g,
// Return an interior arc BAC of triangle ABC, given a radius r > 0.
// If r < 0, return the corresponding exterior arc of radius |r|.
-path arc(explicit pair B, explicit pair A, explicit pair C,
- real r=arrowfactor)
+path arc(explicit pair B, explicit pair A, explicit pair C, real r)
{
return arc(A,r,degrees(B-A),degrees(C-A));
}
diff --git a/Build/source/utils/asymptote/base/graph.asy b/Build/source/utils/asymptote/base/graph.asy
index e9f01880d67..9ac9983ddef 100644
--- a/Build/source/utils/asymptote/base/graph.asy
+++ b/Build/source/utils/asymptote/base/graph.asy
@@ -1,5 +1,4 @@
private import math;
-
import graph_splinetype;
import graph_settings;
@@ -826,14 +825,19 @@ ticks Ticks(int sign, Label F="", ticklabel ticklabel=null,
typedef tickvalues tickmodifier(tickvalues);
tickvalues None(tickvalues v) {return v;}
-tickmodifier OmitTick(... real[] x) {
+// Tickmodifier that removes all ticks in the intervals [a[i],b[i]].
+tickmodifier OmitTickIntervals(real[] a, real[] b) {
return new tickvalues(tickvalues v) {
- void omit(real[] a) {
- if(a.length != 0) {
- real norm=max(abs(a));
- for(int i=0; i < x.length; ++i) {
- int j=find(abs(a-x[i]) < zerotickfuzz*norm);
- if(j >= 0) a.delete(j);
+ if(a.length != b.length) abort(differentlengths);
+ void omit(real[] A) {
+ if(A.length != 0) {
+ real norm=max(abs(A));
+ for(int i=0; i < a.length; ++i) {
+ int j;
+ while((j=find(A > a[i]-zerotickfuzz*norm
+ & A < b[i]+zerotickfuzz*norm)) >= 0) {
+ A.delete(j);
+ }
}
}
}
@@ -843,20 +847,20 @@ tickmodifier OmitTick(... real[] x) {
};
}
-tickmodifier NoZero=OmitTick(0);
+// Tickmodifier that removes all ticks in the interval [a,b].
+tickmodifier OmitTickInterval(real a, real b) {
+ return OmitTickIntervals(new real[] {a}, new real[] {b});
+}
-// Tickmodifier that removes all major ticks in the interval [a,b].
-tickmodifier Break(real a, real b) {
- return new tickvalues(tickvalues v) {
- real[] V=v.major;
- real[] major;
- for(int i=0; i < V.length; ++i)
- if(V[i] < a-epsilon || V[i] > b+epsilon) major.push(V[i]);
- v.major=major;
- return v;
- };
+// Tickmodifier that removes the specified ticks.
+tickmodifier OmitTick(... real[] x) {
+ return OmitTickIntervals(x,x);
}
+tickmodifier NoZero=OmitTick(0);
+
+tickmodifier Break(real, real)=OmitTickInterval;
+
// Automatic tick construction routine.
ticks Ticks(int sign, Label F="", ticklabel ticklabel=null,
bool beginlabel=true, bool endlabel=true,
@@ -1536,7 +1540,7 @@ void xaxis(picture pic=currentpicture, Label L="", axis axis=YZero,
void yaxis(picture pic=currentpicture, Label L="", axis axis=XZero,
real ymin=-infinity, real ymax=infinity, pen p=currentpen,
ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin,
- bool above=false)
+ bool above=false, bool autorotate=true)
{
if(ymin > ymax) return;
@@ -1590,7 +1594,7 @@ void yaxis(picture pic=currentpicture, Label L="", axis axis=XZero,
if(L.defaultposition) L.position(axis.position);
L.align(L.align,axis.align);
- if(L.defaulttransform) {
+ if(autorotate && L.defaulttransform) {
frame f;
add(f,Label(L.s,(0,0),L.p));
if(length(max(f)-min(f)) > ylabelwidth*fontsize(L.p))
@@ -1977,15 +1981,8 @@ guide[] graph(picture pic=currentpicture, pair z(real), real a, real b,
},a,b,n);
}
-string differentlengths="attempt to graph arrays of different lengths";
string conditionlength="condition array has different length than data";
-void checklengths(int x, int y, string text=differentlengths)
-{
- if(x != y)
- abort(text+": "+string(x)+" != "+string(y));
-}
-
void checkconditionlength(int x, int y)
{
checklengths(x,y,conditionlength);
@@ -2070,6 +2067,19 @@ guide polargraph(picture pic=currentpicture, real r(real), real a, real b,
},a,b,n);
}
+guide polargraph(picture pic=currentpicture, real[] r, real[] theta,
+ interpolate join=operator--)
+{
+ int n=r.length;
+ checklengths(n,theta.length);
+ int i=0;
+ return graph(join)(new pair(real) {
+ pair w=Scale(pic,polar(r[i],theta[i]));
+ ++i;
+ return w;
+ },0,0,n-1);
+}
+
void errorbar(picture pic, pair z, pair dp, pair dm, pen p=currentpen,
real size=0)
{
@@ -2166,7 +2176,7 @@ picture vectorfield(path vector(pair), pair a, pair b,
path g=vector(z);
return abs(point(g,size(g)-1)-point(g,0));
}
- real max=size((0,0));
+ real max=size(a);
for(int i=0; i <= nx; ++i) {
real x=interp(a.x,b.x,i*dx);
for(int j=0; j <= ny; ++j)
@@ -2198,8 +2208,9 @@ path Arc(pair c, real r, real angle1, real angle2, bool direction,
{
angle1=radians(angle1);
angle2=radians(angle2);
- if(angle1 >= angle2 && direction) angle1 -= 2pi;
- if(angle2 >= angle1 && !direction) angle2 -= 2pi;
+ if(direction) {
+ if(angle1 >= angle2) angle1 -= 2pi;
+ } else if(angle2 >= angle1) angle2 -= 2pi;
return shift(c)*polargraph(new real(real t){return r;},angle1,angle2,n,
operator ..);
}
diff --git a/Build/source/utils/asymptote/base/graph3.asy b/Build/source/utils/asymptote/base/graph3.asy
index 1357f3e7c87..32f66851e75 100644
--- a/Build/source/utils/asymptote/base/graph3.asy
+++ b/Build/source/utils/asymptote/base/graph3.asy
@@ -126,7 +126,7 @@ void labelaxis(picture pic, transform3 T, Label L, path3 g,
locate1.dir(T,g,locate,t);
triple pathdir=locate1.pathdir;
- triple perp=cross(pathdir,P.vector());
+ triple perp=cross(pathdir,P.normal);
if(align == O)
align=unit(sgn(dot(sign*locate1.dir,perp))*perp);
path[] g=project(box(T*m,T*M),P);
@@ -149,9 +149,11 @@ void labelaxis(picture pic, transform3 T, Label L, path3 g,
},exact=false);
path3[] G=path3(texpath(L));
- G=L.align.is3D ? align(G,O,align,L.p) : L.T3*G;
- triple v=point(g,relative(L,g));
- pic.addBox(v,v,min(G),max(G));
+ if(G.length > 0) {
+ G=L.align.is3D ? align(G,O,align,L.p) : L.T3*G;
+ triple v=point(g,relative(L,g));
+ pic.addBox(v,v,min(G),max(G));
+ }
}
// Tick construction routine for a user-specified array of tick values.
@@ -509,7 +511,7 @@ real ztrans(transform3 t, real z)
private triple defaultdir(triple X, triple Y, triple Z, bool opposite=false,
projection P) {
- triple u=cross(P.vector(),Z);
+ triple u=cross(P.normal,Z);
return abs(dot(u,X)) > abs(dot(u,Y)) ? -X : (opposite ? Y : -Y);
}
@@ -964,7 +966,7 @@ void xaxis3(picture pic=currentpicture, Label L="", axis axis=YZZero,
bool back=false;
if(axis.type == Both) {
- triple v=currentprojection.vector();
+ triple v=currentprojection.normal;
back=dot((0,pic.userMax.y-pic.userMin.y,0),v)*sgn(v.z) > 0;
}
@@ -1039,7 +1041,7 @@ void yaxis3(picture pic=currentpicture, Label L="", axis axis=XZZero,
bool back=false;
if(axis.type == Both) {
- triple v=currentprojection.vector();
+ triple v=currentprojection.normal;
back=dot((pic.userMax.x-pic.userMin.x,0,0),v)*sgn(v.z) > 0;
}
@@ -1842,10 +1844,11 @@ path3 Arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW,
real phi1=radians(longitude(v1,warn=false));
real phi2=radians(longitude(v2,warn=false));
- if(phi1 >= phi2 && direction) phi1 -= 2pi;
- if(phi2 >= phi1 && !direction) phi2 -= 2pi;
+ if(direction) {
+ if(phi1 >= phi2) phi1 -= 2pi;
+ } else if(phi2 >= phi1) phi2 -= 2pi;
- real piby2=pi/2;
+ static real piby2=pi/2;
return shift(c)*T*polargraph(new real(real theta, real phi) {return r;},
new real(real t) {return piby2;},
new real(real t) {return interp(phi1,phi2,t);},
@@ -1869,5 +1872,10 @@ path3 Arc(triple c, real r, real theta1, real phi1, real theta2, real phi2,
// True circle
path3 Circle(triple c, real r, triple normal=Z, int n=nCircle)
{
- return Arc(c,r,90,0,90,360,normal,n)&cycle;
+ static real piby2=pi/2;
+ return shift(c)*align(unit(normal))*
+ polargraph(new real(real theta, real phi) {return r;},
+ new real(real t) {return piby2;},
+ new real(real t) {return interp(0,2pi,t);},n,operator ..);
+
}
diff --git a/Build/source/utils/asymptote/base/graph_splinetype.asy b/Build/source/utils/asymptote/base/graph_splinetype.asy
index fbcb1b00c75..02f780ca02e 100644
--- a/Build/source/utils/asymptote/base/graph_splinetype.asy
+++ b/Build/source/utils/asymptote/base/graph_splinetype.asy
@@ -1,3 +1,5 @@
+private import math;
+
typedef real[] splinetype(real[], real[]);
restricted real[] Spline(real[] x, real[] y);
@@ -11,6 +13,12 @@ void checklengths(int x, int y, string text=differentlengths)
abort(text+": "+string(x)+" != "+string(y));
}
+void checkincreasing(real[] x)
+{
+ if(!increasing(x,true))
+ abort("strictly increasing array expected");
+}
+
// Linear interpolation
real[] linear(real[] x, real[] y)
{
@@ -32,6 +40,7 @@ real[] notaknot(real[] x, real[] y)
{
int n=x.length;
checklengths(n,y.length);
+ checkincreasing(x);
real[] d;
if(n > 3) {
real[] a=new real[n];
@@ -75,6 +84,7 @@ real[] periodic(real[] x, real[] y)
{
int n=x.length;
checklengths(n,y.length);
+ checkincreasing(x);
if(abs(y[n-1]-y[0]) > sqrtEpsilon*norm(y))
abort("function values are not periodic");
real[] d;
@@ -110,6 +120,7 @@ real[] natural(real[] x, real[] y)
{
int n=x.length;
checklengths(n,y.length);
+ checkincreasing(x);
real[] d;
if(n > 2) {
real[] a=new real[n];
@@ -144,6 +155,7 @@ splinetype clamped(real slopea, real slopeb)
return new real[] (real[] x, real[] y) {
int n=x.length;
checklengths(n,y.length);
+ checkincreasing(x);
real[] d;
if(n > 2) {
real[] a=new real[n];
@@ -183,6 +195,7 @@ real[] monotonic(real[] x, real[] y)
{
int n=x.length;
checklengths(n,y.length);
+ checkincreasing(x);
real[] d=new real[n];
if(n > 2) {
real[] h=new real[n-1];
diff --git a/Build/source/utils/asymptote/base/math.asy b/Build/source/utils/asymptote/base/math.asy
index 85d41d86d07..eef703cb7ee 100644
--- a/Build/source/utils/asymptote/base/math.asy
+++ b/Build/source/utils/asymptote/base/math.asy
@@ -38,6 +38,15 @@ bool polygon(path p)
return cyclic(p) && piecewisestraight(p);
}
+// Return the intersection time of the point on the line through p and q
+// that is closest to z.
+real intersect(pair p, pair q, pair z)
+{
+ pair u=q-p;
+ real denom=dot(u,u);
+ return denom == 0 ? infinity : dot(z-p,u)/denom;
+}
+
// Return the intersection time of the extension of the line segment PQ
// with the plane perpendicular to n and passing through Z.
real intersect(triple P, triple Q, triple n, triple Z)
@@ -185,6 +194,30 @@ real[][] operator /(real[][] a, real b)
return a*(1/b);
}
+private string incommensurate=
+ "Multiplication of incommensurate matrices is undefined";
+
+pair[][] operator * (pair[][] a, pair[][] b)
+{
+ int n=a.length;
+ int nb=b.length;
+ int nb0=b[0].length;
+ pair[][] m=new pair[n][nb0];
+ for(int i=0; i < n; ++i) {
+ pair[] ai=a[i];
+ pair[] mi=m[i];
+ if(ai.length != nb)
+ abort(incommensurate);
+ for(int j=0; j < nb0; ++j) {
+ pair sum;
+ for(int k=0; k < nb; ++k)
+ sum += ai[k]*b[k][j];
+ mi[j]=sum;
+ }
+ }
+ return m;
+}
+
bool square(real[][] m)
{
int n=m.length;
@@ -371,3 +404,20 @@ pair[] quarticroots(real a, real b, real c, real d, real e)
return roots;
}
+
+pair[][] fft(pair[][] a, int sign=1)
+{
+ pair[][] A=new pair[a.length][];
+ int k=0;
+ for(pair[] v : a) {
+ A[k]=fft(v,sign);
+ ++k;
+ }
+ a=transpose(A);
+ k=0;
+ for(pair[] v : a) {
+ A[k]=fft(v,sign);
+ ++k;
+ }
+ return transpose(A);
+}
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index 77ab8d6e589..0b88060e081 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -244,17 +244,31 @@ real adjust(real h, real error, real tolmin, real tolmax, RKTableau tableau)
return h;
}
+struct solution
+{
+ real[] t;
+ real[] y;
+}
+
+void write(solution S)
+{
+ for(int i=0; i < S.t.length; ++i)
+ write(S.t[i],S.y[i]);
+}
+
// Integrate dy/dt+cy=f(t,y) from a to b using initial conditions y,
// specifying either the step size h or the number of steps n.
-real[] integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
- real h=0, int n=0, bool dynamic=false, real tolmin=0,
- real tolmax=0, real dtmin=0, real dtmax=realMax,
- RKTableau tableau, bool verbose=false)
+solution integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
+ real h=0, int n=0, bool dynamic=false, real tolmin=0,
+ real tolmax=0, real dtmin=0, real dtmax=realMax,
+ RKTableau tableau, bool verbose=false)
{
- real[] Y={y};
+ solution S;
+ S.t=new real[] {a};
+ S.y=new real[]{y};
if(h == 0) {
- if(b == a) return Y;
+ if(b == a) return S;
if(n == 0) abort("Either n or h must be specified");
else h=(b-a)/n;
}
@@ -302,30 +316,49 @@ real[] integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
if(h >= dt) {
t += dt;
y=y0+highOrder;
- Y.push(y);
+ S.t.push(t);
+ S.y.push(y);
f0=f1;
}
h=min(max(h,dtmin),dtmax);
} else {
t += h;
y=y0+highOrder;
- Y.push(y);
+ S.y.push(y);
}
}
- return Y;
+ return S;
+}
+
+struct Solution
+{
+ real[] t;
+ real[][] y;
+}
+
+void write(Solution S)
+{
+ for(int i=0; i < S.t.length; ++i) {
+ write(S.t[i],tab);
+ for(real y : S.y[i])
+ write(y,tab);
+ write();
+ }
}
// Integrate a set of equations, dy/dt=f(t,y), from a to b using initial
// conditions y, specifying either the step size h or the number of steps n.
-real[][] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
+Solution integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
real h=0, int n=0, bool dynamic=false,
real tolmin=0, real tolmax=0, real dtmin=0,
real dtmax=realMax, RKTableau tableau, bool verbose=false)
{
- real[][] Y={copy(y)};
+ Solution S;
+ S.t=new real[] {a};
+ S.y=new real[][] {copy(y)};
if(h == 0) {
- if(b == a) return Y;
+ if(b == a) return S;
if(n == 0) abort("Either n or h must be specified");
else h=(b-a)/n;
}
@@ -366,17 +399,19 @@ real[][] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
if(h >= dt) {
t += dt;
y += highOrder;
- Y.push(y);
+ S.t.push(t);
+ S.y.push(y);
f0=f1;
}
h=min(max(h,dtmin),dtmax);
} else {
t += h;
y += highOrder;
- Y.push(y);
+ S.t.push(t);
+ S.y.push(y);
}
}
- return Y;
+ return S;
}
real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t,
@@ -408,12 +443,14 @@ real[] newton(int iterations=100, real[] f(real[]), real[][] jacobian(real[]),
}
real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0,
+ bool dynamic=false, real tolmin=0, real tolmax=0, real dtmin=0,
+ real dtmax=realMax, RKTableau tableau, bool verbose=false,
real[] initial(real[]), real[] discrepancy(real[]),
- real[] guess, RKTableau tableau, int iterations=100)
+ real[] guess, int iterations=100)
{
real[] g(real[] t) {
- real[][] y=integrate(initial(t),f,a,b,h,n,tableau);
- return discrepancy(y[y.length-1]);
+ real[][] y=integrate(initial(t),f,a,b,h,n,dynamic,tolmin,tolmax,dtmin,dtmax,
+ tableau,verbose).y;return discrepancy(y[y.length-1]);
}
real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);}
return initial(newton(iterations,g,jacobian,guess));
diff --git a/Build/source/utils/asymptote/base/palette.asy b/Build/source/utils/asymptote/base/palette.asy
index fabbdc2cd5b..41a40a04a85 100644
--- a/Build/source/utils/asymptote/base/palette.asy
+++ b/Build/source/utils/asymptote/base/palette.asy
@@ -305,7 +305,7 @@ pen[] Wheel(int NColors=32766)
if(settings.gray) return Grayscale(NColors);
int nintervals=6;
- int n=quotient(NColors,nintervals);
+ int n=-quotient(NColors,-nintervals);
pen[] Palette;
if(n == 0) return Palette;
@@ -333,7 +333,7 @@ pen[] Rainbow(int NColors=32766)
int offset=1;
int nintervals=5;
- int n=quotient(NColors-1,nintervals);
+ int n=-quotient(NColors-1,-nintervals);
pen[] Palette;
if(n == 0) return Palette;
@@ -366,7 +366,7 @@ private pen[] BWRainbow(int NColors, bool two)
if(two) nintervals += 6;
int num=NColors-offset;
- int n=quotient(num,nintervals*divisor)*divisor;
+ int n=-quotient(num,-nintervals*divisor)*divisor;
NColors=n*nintervals+offset;
pen[] Palette;
@@ -395,7 +395,7 @@ private pen[] BWRainbow(int NColors, bool two)
for(int i=0; i < n; ++i)
Palette[k+i]=rgb(1.0-i*ninv,0.0,1.0);
else {
- int n3=quotient(n,3);
+ int n3=-quotient(n,-3);
int n23=2*n3;
real third=n3*ninv;
real twothirds=n23*ninv;
diff --git a/Build/source/utils/asymptote/base/plain.asy b/Build/source/utils/asymptote/base/plain.asy
index f9e18b883ea..cc1a5064038 100644
--- a/Build/source/utils/asymptote/base/plain.asy
+++ b/Build/source/utils/asymptote/base/plain.asy
@@ -24,6 +24,7 @@ if(version.VERSION != VERSION) {
version.VERSION+" of plain.asy"+'\n');
}
+include plain_strings;
include plain_pens;
include plain_paths;
include plain_filldraw;
@@ -33,7 +34,6 @@ include plain_Label;
include plain_shipout;
include plain_arcs;
include plain_boxes;
-include plain_strings;
include plain_markers;
include plain_arrows;
include plain_debugger;
diff --git a/Build/source/utils/asymptote/base/plain_Label.asy b/Build/source/utils/asymptote/base/plain_Label.asy
index b75d0de3967..6e6f231375a 100644
--- a/Build/source/utils/asymptote/base/plain_Label.asy
+++ b/Build/source/utils/asymptote/base/plain_Label.asy
@@ -75,11 +75,19 @@ transform scaleless(transform t)
A=conj(U)*A*U;
real D=abs(A[0][0]);
- if(D != 0) A[0][0] /= D;
+ if(D != 0) {
+ A[0][0] /= D;
+ A[0][1] /= D;
+ }
+
D=abs(A[1][1]);
- if(D != 0) A[1][1] /= D;
+ if(D != 0) {
+ A[1][0] /= D;
+ A[1][1] /= D;
+ }
A=U*A*conj(U);
+
return (0,0,A[0][0].x,A[0][1].x,A[1][0].x,A[1][1].x);
}
@@ -209,7 +217,8 @@ transform Slant(transform t) {return scaleless(t);}
transform Scale(transform t) {return t;}
embed Rotate(pair z) {
- return new transform(transform t) {return rotate(degrees(shiftless(t)*z));};
+ return new transform(transform t) {return rotate(degrees(shiftless(t)*z,
+ warn=false));};
}
struct Label {
@@ -288,7 +297,7 @@ struct Label {
void label(frame f, transform t=identity(), pair position, pair align) {
pen p0=p == nullpen ? currentpen : p;
- align=length(align)*unit(scaleless(shiftless(t))*align);
+ align=length(align)*unit(shiftless(t)*align);
label(f,s,size,embed(t)*shiftless(T),
t*position+align*labelmargin(p0)+shift(T)*0,align,p0);
}
diff --git a/Build/source/utils/asymptote/base/plain_arcs.asy b/Build/source/utils/asymptote/base/plain_arcs.asy
index 01e750f42cd..140bc3cff19 100644
--- a/Build/source/utils/asymptote/base/plain_arcs.asy
+++ b/Build/source/utils/asymptote/base/plain_arcs.asy
@@ -23,8 +23,9 @@ path arc(pair c, explicit pair z1, explicit pair z2, bool direction=CCW)
real t1=intersect(unitcircle,(0,0)--2*z1)[0];
real t2=intersect(unitcircle,(0,0)--2*z2)[0];
static int n=length(unitcircle);
- if(t1 >= t2 && direction) t1 -= n;
- if(t2 >= t1 && !direction) t2 -= n;
+ if(direction) {
+ if (t1 >= t2) t1 -= n;
+ } else if(t2 >= t1) t2 -= n;
return shift(c)*scale(r)*subpath(unitcircle,t1,t2);
}
diff --git a/Build/source/utils/asymptote/base/plain_boxes.asy b/Build/source/utils/asymptote/base/plain_boxes.asy
index ff82a0b0b32..0c6401311d8 100644
--- a/Build/source/utils/asymptote/base/plain_boxes.asy
+++ b/Build/source/utils/asymptote/base/plain_boxes.asy
@@ -4,7 +4,8 @@ path box(frame dest, frame src=dest, real xmargin=0, real ymargin=xmargin,
{
pair z=(xmargin,ymargin);
int sign=filltype == NoFill ? 1 : -1;
- path g=box(min(src)+0.5*sign*min(p)-z,max(src)+0.5*sign*max(p)+z);
+ pair h=0.5*sign*(max(p)-min(p));
+ path g=box(min(src)-h-z,max(src)+h+z);
frame F;
if(above == false) {
filltype.fill(F,g,p);
@@ -45,8 +46,8 @@ path ellipse(frame dest, frame src=dest, real xmargin=0, real ymargin=xmargin,
pair D=M-m;
static real factor=0.5*sqrt(2);
int sign=filltype == NoFill ? 1 : -1;
- path g=ellipse(0.5*(M+m),factor*D.x+0.5*sign*max(p).x+xmargin,
- factor*D.y+0.5*sign*max(p).y+ymargin);
+ pair h=0.5*sign*(max(p)-min(p));
+ path g=ellipse(0.5*(M+m),factor*D.x+h.x+xmargin,factor*D.y+h.y+ymargin);
frame F;
if(above == false) {
filltype.fill(F,g,p);
diff --git a/Build/source/utils/asymptote/base/plain_constants.asy b/Build/source/utils/asymptote/base/plain_constants.asy
index 32138756891..ba8947c91c2 100644
--- a/Build/source/utils/asymptote/base/plain_constants.asy
+++ b/Build/source/utils/asymptote/base/plain_constants.asy
@@ -8,10 +8,10 @@ restricted real bp=1; // A PostScript point.
restricted real pt=72.0/72.27; // A TeX pt; smaller than a PostScript bp.
restricted pair I=(0,1);
-restricted pair up=(0,1);
-restricted pair down=(0,-1);
restricted pair right=(1,0);
restricted pair left=(-1,0);
+restricted pair up=(0,1);
+restricted pair down=(0,-1);
restricted pair E=(1,0);
restricted pair N=(0,1);
diff --git a/Build/source/utils/asymptote/base/plain_paths.asy b/Build/source/utils/asymptote/base/plain_paths.asy
index bf3f2790c1d..3dab68fe32e 100644
--- a/Build/source/utils/asymptote/base/plain_paths.asy
+++ b/Build/source/utils/asymptote/base/plain_paths.asy
@@ -15,6 +15,7 @@ tensionSpecifier operator tension(real t, bool atLeast)
{
return operator tension(t,t,atLeast);
}
+
guide operator controls(pair z)
{
return operator controls(z,z);
@@ -92,19 +93,26 @@ void write(file file, string s="", explicit path[] x, suffix suffix=none)
write(file,suffix);
}
-void write(file file, string s="", explicit guide[] x, suffix suffix=none)
+void write(string s="", explicit path[] x, suffix suffix=endl)
{
- write(file,s,(path[]) x,suffix);
+ write(stdout,s,x,suffix);
}
-void write(string s="", explicit path[] x, suffix suffix=endl)
+void write(file file, string s="", explicit guide[] x, suffix suffix=none)
{
- write(stdout,s,x,suffix);
+ write(file,s);
+ if(x.length > 0) write(file,x[0]);
+ for(int i=1; i < x.length; ++i) {
+ write(file,endl);
+ write(file," ^^");
+ write(file,x[i]);
+ }
+ write(file,suffix);
}
void write(string s="", explicit guide[] x, suffix suffix=endl)
{
- write(stdout,s,(path[]) x,suffix);
+ write(stdout,s,x,suffix);
}
private string nopoints="nullpath has no points";
@@ -127,23 +135,19 @@ pair max(explicit path[] p)
return maxp;
}
-guide operator ::(... guide[] a)
+interpolate operator ..(tensionSpecifier t)
{
- if(a.length == 0) return nullpath;
- guide g=a[0];
- for(int i=1; i < a.length; ++i)
- g=g..operator tension(1,true)..a[i];
- return g;
+ return new guide(... guide[] a) {
+ if(a.length == 0) return nullpath;
+ guide g=a[0];
+ for(int i=1; i < a.length; ++i)
+ g=g..t..a[i];
+ return g;
+ };
}
-guide operator ---(... guide[] a)
-{
- if(a.length == 0) return nullpath;
- guide g=a[0];
- for(int i=1; i < a.length; ++i)
- g=g..operator tension(infinity,true)..a[i];
- return g;
-}
+interpolate operator ::=operator ..(operator tension(1,true));
+interpolate operator ---=operator ..(operator tension(infinity,true));
// return an arbitrary intersection point of paths p and q
pair intersectionpoint(path p, path q, real fuzz=-1)
@@ -283,7 +287,7 @@ path buildcycle(... path[] p)
for(int i=0; i < n; ++i) {
real[][] t=intersections(p[i],p[j]);
if(t.length == 0)
- return nullpath;
+ return nullpath;
ta[i]=t[0][0]; tb[j]=t[0][1];
j=i;
}
@@ -302,7 +306,7 @@ path buildcycle(... path[] p)
int L=length(p[i]);
real t=Tb-L;
if(abs(c-point(p[i],0.5(Ta+t))) <
- abs(c-point(p[i],0.5(Ta+Tb)))) Tb=t;
+ abs(c-point(p[i],0.5(Ta+Tb)))) Tb=t;
while(Tb < Ta) Tb += L;
}
G=G&subpath(p[i],Ta,Tb);
diff --git a/Build/source/utils/asymptote/base/plain_pens.asy b/Build/source/utils/asymptote/base/plain_pens.asy
index a8a0d0c7461..5b0bdbd5879 100644
--- a/Build/source/utils/asymptote/base/plain_pens.asy
+++ b/Build/source/utils/asymptote/base/plain_pens.asy
@@ -1,12 +1,17 @@
real labelmargin=0.3;
real dotfactor=6;
-pen solid=linetype("");
-pen dotted=linetype("0 4");
-pen dashed=linetype("8 8");
-pen longdashed=linetype("24 8");
-pen dashdotted=linetype("8 8 0 8");
-pen longdashdotted=linetype("24 8 0 8");
+pen solid=linetype(new real[]);
+pen dotted=linetype(new real[] {0,4});
+pen dashed=linetype(new real[] {8,8});
+pen longdashed=linetype(new real[] {24,8});
+pen dashdotted=linetype(new real[] {8,8,0,8});
+pen longdashdotted=linetype(new real[] {24,8,0,8});
+
+pen linetype(string pattern, real offset=0, bool scale=true, bool adjust=true)
+{
+ return linetype((real[]) split(pattern),offset,scale,adjust);
+}
void defaultpen(real w) {defaultpen(linewidth(w));}
pen operator +(pen p, real w) {return p+linewidth(w);}
diff --git a/Build/source/utils/asymptote/base/plain_picture.asy b/Build/source/utils/asymptote/base/plain_picture.asy
index e13120a198c..56ebf9f4bbe 100644
--- a/Build/source/utils/asymptote/base/plain_picture.asy
+++ b/Build/source/utils/asymptote/base/plain_picture.asy
@@ -4,11 +4,6 @@ pair viewportsize=0; // Horizontal and vertical viewport limits.
restricted bool Aspect=true;
restricted bool IgnoreAspect=false;
-pair size(frame f)
-{
- return max(f)-min(f);
-}
-
typedef real[][] transform3;
restricted transform3 identity4=identity(4);
@@ -360,6 +355,7 @@ pair point(frame f, pair dir)
path[] align(path[] g, transform t=identity(), pair position,
pair align, pen p=currentpen)
{
+ if(g.length == 0) return g;
pair m=min(g);
pair M=max(g);
pair dir=rectify(inverse(t)*-align);
@@ -385,18 +381,15 @@ struct transformation {
transform3 modelview; // For orientation and positioning
transform3 projection; // For 3D to 2D projection
bool infinity;
- bool oblique;
- void operator init(transform3 modelview, bool oblique=false) {
+ void operator init(transform3 modelview) {
this.modelview=modelview;
this.projection=identity4;
infinity=true;
- this.oblique=oblique;
}
void operator init(transform3 modelview, transform3 projection) {
this.modelview=modelview;
this.projection=projection;
infinity=false;
- oblique=false;
}
transform3 compute() {
return infinity ? modelview : projection*modelview;
@@ -406,7 +399,6 @@ struct transformation {
T.modelview=copy(modelview);
T.projection=copy(projection);
T.infinity=infinity;
- T.oblique=oblique;
return T;
}
}
@@ -414,11 +406,11 @@ struct transformation {
struct projection {
transform3 t; // projection*modelview (cached)
bool infinity;
- bool oblique;
bool absolute=false;
triple camera; // Position of camera.
triple up; // A vector that should be projected to direction (0,1).
triple target; // Point where camera is looking at.
+ triple normal; // Normal vector from target to projection plane.
pair viewportshift; // Fractional viewport shift.
real zoom=1; // Zoom factor.
real angle; // Lens angle (for perspective projection).
@@ -434,7 +426,6 @@ struct projection {
T=projector(camera,up,target);
t=T.compute();
infinity=T.infinity;
- oblique=T.oblique;
ninterpolate=infinity ? 1 : 16;
}
@@ -443,12 +434,14 @@ struct projection {
}
void operator init(triple camera, triple up=(0,0,1), triple target=(0,0,0),
+ triple normal=camera-target,
real zoom=1, real angle=0, pair viewportshift=0,
bool showtarget=true, bool autoadjust=true,
bool center=false, projector projector) {
this.camera=camera;
this.up=up;
this.target=target;
+ this.normal=normal;
this.zoom=zoom;
this.angle=angle;
this.viewportshift=viewportshift;
@@ -464,10 +457,10 @@ struct projection {
P.t=t;
P.infinity=infinity;
P.absolute=absolute;
- P.oblique=oblique;
P.camera=camera;
P.up=up;
P.target=target;
+ P.normal=normal;
P.zoom=zoom;
P.angle=angle;
P.viewportshift=viewportshift;
@@ -1281,7 +1274,7 @@ struct picture {
}
- // Calculate additional scaling required if only an approximate thisture
+ // Calculate additional scaling required if only an approximate picture
// size estimate is available.
transform3 scale3(frame f, real xsize3=this.xsize3,
real ysize3=this.ysize3, real zsize3=this.zsize3,
@@ -1567,12 +1560,24 @@ void add(picture pic=currentpicture, drawer d, bool exact=false)
pic.add(d,exact);
}
+typedef void drawer3(frame f, transform3 t, picture pic, projection P);
+void add(picture pic=currentpicture, drawer3 d, bool exact=false)
+{
+ pic.add(d,exact);
+}
+
void add(picture pic=currentpicture, void d(picture,transform),
bool exact=false)
{
pic.add(d,exact);
}
+void add(picture pic=currentpicture, void d(picture,transform3),
+ bool exact=false)
+{
+ pic.add(d,exact);
+}
+
void begingroup(picture pic=currentpicture)
{
pic.add(new void(frame f, transform) {
@@ -1628,7 +1633,7 @@ void latticeshade(picture pic=currentpicture, path[] g, bool stroke=false,
p=copy(p);
}
pic.add(new void(frame f, transform t) {
- latticeshade(f,t*g,stroke,fillrule,p,false);
+ latticeshade(f,t*g,stroke,fillrule,p,t,false);
},true);
pic.addPath(g);
}
@@ -1738,7 +1743,8 @@ void tensorshade(picture pic=currentpicture, path[] g, bool stroke=false,
// Smoothly shade the regions between consecutive paths of a sequence using a
// given array of pens:
-void draw(picture pic=currentpicture, path[] g, pen[] p)
+void draw(picture pic=currentpicture, path[] g, pen fillrule=currentpen,
+ pen[] p)
{
path[] G;
pen[][] P;
@@ -1753,7 +1759,7 @@ void draw(picture pic=currentpicture, path[] g, pen[] p)
P.push(new pen[] {p[i],p[i],p[i+1],p[i+1]});
}
}
- tensorshade(pic,G,P);
+ tensorshade(pic,G,fillrule,P);
}
void functionshade(picture pic=currentpicture, path[] g, bool stroke=false,
@@ -1883,7 +1889,14 @@ void add(picture dest=currentpicture, frame src, pair position, pair align,
add(dest,align(src,align),position,group,filltype,above);
}
-// Like attach(picture,frame,pair) but extend picture to accommodate frame;
+// Like add(frame,frame,pair) but align frame in direction align.
+void add(frame dest, frame src, pair position, pair align,
+ bool group=true, filltype filltype=NoFill, bool above=true)
+{
+ add(dest,align(src,align),position,group,filltype,above);
+}
+
+// Like add(picture,frame,pair,pair) but extend picture to accommodate frame;
void attach(picture dest=currentpicture, frame src, pair position,
pair align, bool group=true, filltype filltype=NoFill,
bool above=true)
@@ -1935,22 +1948,29 @@ void postscript(picture pic=currentpicture, string s)
},true);
}
-void tex(picture pic=currentpicture, string s)
+void postscript(picture pic=currentpicture, string s, pair min, pair max)
{
- pic.add(new void(frame f, transform) {
- tex(f,s);
+ pic.add(new void(frame f, transform t) {
+ postscript(f,s,t*min,t*max);
},true);
}
-void postscript(picture pic=currentpicture, string s, pair min, pair max)
+void tex(picture pic=currentpicture, string s)
{
- pic.add(new void(frame f, transform t) {
- postscript(f,s,t*min,t*max);
+ // Force TeX string s to be evaluated immediately (in case it is a macro).
+ frame g;
+ tex(g,s);
+ size(g);
+ pic.add(new void(frame f, transform) {
+ tex(f,s);
},true);
}
void tex(picture pic=currentpicture, string s, pair min, pair max)
{
+ frame g;
+ tex(g,s);
+ size(g);
pic.add(new void(frame f, transform t) {
tex(f,s,t*min,t*max);
},true);
@@ -1968,3 +1988,31 @@ void erase(picture pic=currentpicture)
pic.uptodate=false;
pic.erase();
}
+
+void begin(picture pic=currentpicture, string name, string id="",
+ bool visible=true)
+{
+ if(!latex() || !pdf()) return;
+ settings.twice=true;
+ if(id == "") id=string(++ocgindex);
+ tex(pic,"\begin{ocg}{"+name+"}{"+id+"}{"+(visible ? "1" : "0")+"}");
+ layer(pic);
+}
+
+void end(picture pic=currentpicture)
+{
+ if(!latex() || !pdf()) return;
+ tex(pic,"\end{ocg}");
+ layer(pic);
+}
+
+// For users of the LaTeX babel package.
+void deactivatequote(picture pic=currentpicture)
+{
+ tex(pic,"\catcode`\"=12");
+}
+
+void activatequote(picture pic=currentpicture)
+{
+ tex(pic,"\catcode`\"=13");
+}
diff --git a/Build/source/utils/asymptote/base/plain_shipout.asy b/Build/source/utils/asymptote/base/plain_shipout.asy
index 9691a8e2fcf..e96a3a8acb1 100644
--- a/Build/source/utils/asymptote/base/plain_shipout.asy
+++ b/Build/source/utils/asymptote/base/plain_shipout.asy
@@ -2,9 +2,7 @@
string defaultfilename;
string outprefix(string prefix=defaultfilename) {
- string s=prefix != "" ? prefix :
- (settings.outname == "" && interactive()) ? "out" : settings.outname;
- return stripextension(s);
+ return stripextension(prefix != "" ? prefix : outname());
}
string outformat(string format="")
@@ -14,7 +12,6 @@ string outformat(string format="")
return format;
}
-
bool shipped; // Was a picture or frame already shipped out?
frame currentpatterns;
diff --git a/Build/source/utils/asymptote/base/plain_strings.asy b/Build/source/utils/asymptote/base/plain_strings.asy
index 598d61e4249..138ab05ca05 100644
--- a/Build/source/utils/asymptote/base/plain_strings.asy
+++ b/Build/source/utils/asymptote/base/plain_strings.asy
@@ -83,19 +83,28 @@ string verbatim(string s)
}
// Split a string into an array of substrings delimited by delimiter
-string[] split(string s, string delimiter=" ")
+// If delimiter is an empty string, use space delimiter, discarding duplicates.
+string[] split(string s, string delimiter="")
{
+ bool prune=false;
+ if(delimiter == "") {
+ prune=true;
+ delimiter=" ";
+ }
+
string[] S;
int last=0;
int i;
int N=length(delimiter);
+ int n=length(s);
while((i=find(s,delimiter,last)) >= 0) {
if(i >= last)
S.push(substr(s,last,i-last));
last=i+N;
+ if(prune)
+ while(substr(s,last,1) == " ") ++last;
}
- int n=length(s);
- if(n >= last)
+ if(n > last || (n == last && !prune))
S.push(substr(s,last,n-last));
return S;
}
@@ -150,7 +159,12 @@ string graphic(string name, string options="")
{
if(latex()) {
if(options != "") options="["+options+"]";
- return "\includegraphics"+options+"{"+name+"}";
+ if(find(name," ") < 0)
+ return "\includegraphics"+options+"{"+name+"}";
+ else {
+ return "\includegraphics"+options+
+ (pdf() ? "{\""+stripextension(name)+"\".pdf}" : "{\""+name+"\"}");
+ }
}
if(settings.tex != "context")
notimplemented("graphic");
@@ -160,10 +174,10 @@ string graphic(string name, string options="")
string minipage(string s, real width=100bp)
{
if(latex())
- return "\begin{minipage}{"+(string) (width*pt)+"pt}"+s+"\end{minipage}";
+ return "\begin{minipage}{"+(string) (width/pt)+"pt}"+s+"\end{minipage}";
if(settings.tex != "context")
notimplemented("minipage");
- return "\startframedtext[none][frame=off,width="+(string) (width*pt)+
+ return "\startframedtext[none][frame=off,width="+(string) (width/pt)+
"pt]"+s+"\stopframedtext";
}
@@ -197,20 +211,3 @@ string phantom(string s)
}
restricted int ocgindex=0;
-
-void begin(picture pic=currentpicture, string name, string id="",
- bool visible=true)
-{
- if(!latex() || !pdf()) return;
- settings.twice=true;
- if(id == "") id=string(++ocgindex);
- tex(pic,"\begin{ocg}{"+name+"}{"+id+"}{"+(visible ? "1" : "0")+"}");
- layer(pic);
-}
-
-void end(picture pic=currentpicture)
-{
- if(!latex() || !pdf()) return;
- tex(pic,"\end{ocg}");
- layer(pic);
-}
diff --git a/Build/source/utils/asymptote/base/slide.asy b/Build/source/utils/asymptote/base/slide.asy
index 579253b182f..dc1608203e1 100644
--- a/Build/source/utils/asymptote/base/slide.asy
+++ b/Build/source/utils/asymptote/base/slide.asy
@@ -526,8 +526,8 @@ void item(string s, pen p=itempen, bool step=itemstep)
frame b;
label(b,bullet,(0,0),p);
real bulletwidth=max(b).x-min(b).x;
- remark(bullet+"\hangindent"+(string) bulletwidth+"pt$\,$"+s,p,
- -bulletwidth*pt,step=step);
+ remark(bullet+"\hangindent"+(string) (bulletwidth/pt)+"pt$\,$"+s,p,
+ -bulletwidth,step=step);
}
void subitem(string s, pen p=itempen)
diff --git a/Build/source/utils/asymptote/base/solids.asy b/Build/source/utils/asymptote/base/solids.asy
index db98af3088c..7e5883bc884 100644
--- a/Build/source/utils/asymptote/base/solids.asy
+++ b/Build/source/utils/asymptote/base/solids.asy
@@ -1,7 +1,7 @@
import three;
import graph3;
-pen defaultbackpen=linetype("4 4",4,scale=false);
+pen defaultbackpen=linetype(new real[] {4,4},4,scale=false);
// A solid geometry package.
diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy
index 5c3a404514f..cc7534f6f79 100644
--- a/Build/source/utils/asymptote/base/three.asy
+++ b/Build/source/utils/asymptote/base/three.asy
@@ -17,6 +17,7 @@ real linegranularity=0.005;
int linesectors=8; // Number of angular sectors.
real dotgranularity=0.0001;
real angleprecision=1e-5; // Precision for centering perspective projections.
+int maxangleiterations=25;
real rendermargin=0.02;
string defaultembed3Doptions;
@@ -190,7 +191,9 @@ projection operator * (transform3 t, projection P)
projection P=P.copy();
if(!P.absolute) {
P.camera=t*P.camera;
+ P.normal=t*(P.target+P.normal);
P.target=t*P.target;
+ P.normal -= P.target;
P.calculate();
}
return P;
@@ -268,9 +271,10 @@ projection oblique(real angle=45)
t[0][2]=-c2;
t[1][2]=-s2;
t[2][2]=1;
- return projection((c2,s2,1),up=Y,
+ t[2][3]=-1;
+ return projection((c2,s2,1),up=Y,normal=(0,0,1),
new transformation(triple,triple,triple) {
- return transformation(t,oblique=true);});
+ return transformation(t);});
}
projection obliqueZ(real angle=45) {return oblique(angle);}
@@ -287,9 +291,10 @@ projection obliqueX(real angle=45)
t[1][2]=1;
t[2][2]=0;
t[2][0]=1;
- return projection((1,c2,s2),
+ t[2][3]=-1;
+ return projection((1,c2,s2),normal=(1,0,0),
new transformation(triple,triple,triple) {
- return transformation(t,oblique=true);});
+ return transformation(t);});
}
projection obliqueY(real angle=45)
@@ -302,9 +307,10 @@ projection obliqueY(real angle=45)
t[1][2]=1;
t[2][1]=-1;
t[2][2]=0;
- return projection((c2,-1,s2),
+ t[2][3]=-1;
+ return projection((c2,-1,s2),normal=(0,-1,0),
new transformation(triple,triple,triple) {
- return transformation(t,oblique=true);});
+ return transformation(t);});
}
projection oblique=oblique();
@@ -349,13 +355,13 @@ triple invert(pair z, triple normal, triple point,
// Map pair to a triple on the projection plane.
triple invert(pair z, projection P=currentprojection)
{
- return invert(z,P.vector(),P.target,P);
+ return invert(z,P.normal,P.target,P);
}
// Map pair dir to a triple direction at point v on the projection plane.
triple invert(pair dir, triple v, projection P=currentprojection)
{
- return invert(project(v,P)+dir,P.vector(),v,P)-v;
+ return invert(project(v,P)+dir,P.normal,v,P)-v;
}
pair xypart(triple v)
@@ -675,24 +681,22 @@ guide3 operator .. (... guide3[] g)
};
}
-guide3 operator ::(... guide3[] a)
-{
- if(a.length == 0) return nullpath3;
- guide3 g=a[0];
- for(int i=1; i < a.length; ++i)
- g=g.. tension atleast 1 ..a[i];
- return g;
-}
+typedef guide3 interpolate3(... guide3[]);
-guide3 operator ---(... guide3[] a)
+interpolate3 join3(tensionSpecifier t)
{
- if(a.length == 0) return nullpath3;
- guide3 g=a[0];
- for(int i=1; i < a.length; ++i)
- g=g.. tension atleast infinity ..a[i];
- return g;
+ return new guide3(... guide3[] a) {
+ if(a.length == 0) return nullpath3;
+ guide3 g=a[0];
+ for(int i=1; i < a.length; ++i)
+ g=g..t..a[i];
+ return g;
+ };
}
+interpolate3 operator ::=join3(operator tension(1,true));
+interpolate3 operator ---=join3(operator tension(infinity,true));
+
flatguide3 operator cast(guide3 g)
{
flatguide3 f;
@@ -1056,12 +1060,12 @@ path3 invert(path p, triple normal, triple point,
path3 invert(path p, triple point, projection P=currentprojection)
{
- return path3(p,new triple(pair z) {return invert(z,P.vector(),point,P);});
+ return path3(p,new triple(pair z) {return invert(z,P.normal,point,P);});
}
path3 invert(path p, projection P=currentprojection)
{
- return path3(p,new triple(pair z) {return invert(z,P.vector(),P.target,P);});
+ return path3(p,new triple(pair z) {return invert(z,P.normal,P.target,P);});
}
// Construct a path from a path3 by applying P to each control point.
@@ -1250,7 +1254,7 @@ path3 solve(flatguide3 g)
path nurb(path3 p, projection P, int ninterpolate=P.ninterpolate)
{
triple f=P.camera;
- triple u=unit(P.vector());
+ triple u=unit(P.normal);
transform3 t=P.t;
path nurb(triple v0, triple v1, triple v2, triple v3) {
@@ -1360,8 +1364,9 @@ triple normal(path3 p)
if(abs(a) >= fuzz && abs(b) >= fuzz) {
triple n=cross(unit(a),unit(b));
real absn=abs(n);
+ if(absn < sqrtEpsilon) return false;
n=unit(n);
- if(absnormal > 0 && absn > sqrtEpsilon &&
+ if(absnormal > 0 &&
abs(normal-n) > sqrtEpsilon && abs(normal+n) > sqrtEpsilon)
return true;
else {
@@ -1455,7 +1460,7 @@ private transform3 flip(transform3 t, triple X, triple Y, triple Z,
return scale(s(v.x),s(v.y),s(v.z));
}
- triple u=unit(P.vector());
+ triple u=unit(P.normal);
triple up=unit(perp(P.up,u));
bool upright=dot(Z,u) >= 0;
if(dot(Y,up) < 0) {
@@ -1782,7 +1787,7 @@ transform3 align(triple u)
// return a rotation that maps X,Y to the projection plane.
transform3 transform3(projection P=currentprojection)
{
- triple w=unit(P.oblique ? P.camera : P.vector());
+ triple w=unit(P.normal);
triple v=unit(perp(P.up,w));
if(v == O) v=cross(perp(w),w);
triple u=cross(v,w);
@@ -1864,8 +1869,6 @@ triple max(explicit path3[] p)
return maxp;
}
-typedef guide3 interpolate3(... guide3[]);
-
path3 randompath3(int n, bool cumulate=true, interpolate3 join=operator ..)
{
guide3 g;
@@ -1948,8 +1951,9 @@ path3 arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW)
real t1=t1[0];
real t2=t2[0];
int n=length(unitcircle3);
- if(t1 >= t2 && direction) t1 -= n;
- if(t2 >= t1 && !direction) t2 -= n;
+ if(direction) {
+ if (t1 >= t2) t1 -= n;
+ } else if(t2 >= t1) t2 -= n;
path3 p=subpath(unitcircle3,t1,t2);
if(align) p=T*p;
@@ -1986,11 +1990,6 @@ path3 plane(triple u, triple v, triple O=O)
return O--O+u--O+u+v--O+v--cycle;
}
-triple size3(frame f)
-{
- return max3(f)-min3(f);
-}
-
// PRC/OpenGL support
include three_light;
@@ -2125,10 +2124,9 @@ draw=new void(frame f, path3 g, material p=currentpen,
} else _draw(f,g,q,name);
} else _draw(f,g,q,name);
}
- string type=linetype(adjust(q,arclength(g),cyclic(g)));
- if(length(type) == 0) drawthick(g);
+ real[] dash=linetype(adjust(q,arclength(g),cyclic(g)));
+ if(dash.length == 0) drawthick(g);
else {
- real[] dash=(real[]) split(type," ");
if(sum(dash) > 0) {
dash.cyclic=true;
real offset=offset(q);
@@ -2418,7 +2416,7 @@ for(i=0; i < count; i++) {
var center=new Array(
";
for(int i=0; i < center.length; ++i)
- s += "Vector3("+format(center[i],",")+"),
+ s += "Vector3("+format(center[i]/cm,",")+"),
";
s += ");
@@ -2457,7 +2455,7 @@ runtime.addEventHandler(billboardHandler);
runtime.refresh();
";
-return s;
+ return s;
}
string lightscript(light light) {
@@ -2674,7 +2672,6 @@ struct scene
// 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;
@@ -2695,7 +2692,7 @@ struct scene
}
diff=abs(s-lasts);
++i;
- } while (diff > angleprecision && i < maxiterations);
+ } while (diff > angleprecision && i < maxangleiterations);
real aspect=width > 0 ? height/width : 1;
real rx=-r.x*aspect;
real Rx=R.x*aspect;
@@ -2767,7 +2764,9 @@ object embed(string label="", string text=label, string prefix=defaultfilename,
bool preview=settings.render > 0;
if(prc) {
// The movie15.sty package cannot handle spaces or dots in filenames.
- prefix=replace(prefix,new string[][]{{" ","_"},{".","_"}});
+ string dir=stripfile(prefix);
+ prefix=dir+replace(stripdirectory(prefix),
+ new string[][]{{" ","_"},{".","_"}});
if(settings.embed || nativeformat() == "pdf")
prefix += "+"+(string) file3.length;
} else
@@ -3063,7 +3062,8 @@ frame[] fit3(string prefix="", picture[] pictures, picture all,
for(int i=settings.reverse ? pictures.length-1 : 0;
i >= 0 && i < pictures.length && !settings.interrupt;
settings.reverse ? --i : ++i) {
- embedder(prefix,out[i],format,view,options,script,light,S.P);
+ frame f=embedder(prefix,out[i],format,view,options,script,light,S.P);
+ if(!settings.loop) out[i]=f;
}
if(!settings.loop) break;
}
diff --git a/Build/source/utils/asymptote/base/three_arrows.asy b/Build/source/utils/asymptote/base/three_arrows.asy
index 161a716bca6..b4319219284 100644
--- a/Build/source/utils/asymptote/base/three_arrows.asy
+++ b/Build/source/utils/asymptote/base/three_arrows.asy
@@ -70,6 +70,7 @@ struct arrowhead3
real size(pen p)=arrowsize;
real arcsize(pen p)=arcarrowsize;
real gap=1;
+ real size;
bool splitpath=true;
surface head(path3 g, position position=EndPoint,
@@ -88,7 +89,7 @@ struct arrowhead3
if(filltype == null) filltype=FillDraw(p);
bool draw=filltype.type != filltype.Fill;
triple v=point(s,length(s));
- triple N=normal == O ? P.vector() : normal;
+ triple N=normal == O ? P.normal : normal;
triple w=unit(v-point(s,0));
transform3 t=transform3(w,unit(cross(w,N)));
path3[] H=t*path3(h);
@@ -180,8 +181,7 @@ arrowhead3 HookHead3(real dir=arrowdir, real barb=arrowbarb)
filltype filltype=null, bool forwards=true,
projection P=currentprojection) {
if(size == 0) size=a.size(p);
-
- angle=min(angle*arrowhookfactor,45);
+
bool relative=position.relative;
real position=position.position.x;
if(relative) position=reltime(g,position);
@@ -265,11 +265,12 @@ arrowhead3 DefaultHead2(triple normal=O) {
projection P=currentprojection) {
if(size == 0) size=a.size(p);
path h=a.project(g,forwards,P);
+ a.size=min(size,arclength(h));
path[] H=a.align(DefaultHead.head(h,p,size,angle),h);
H=forwards ? yscale(-1)*H : H;
return a.surface(g,position,size,H,p,filltype,normal,P);
};
- a.gap=1.01;
+ a.gap=1.005;
return a;
}
arrowhead3 DefaultHead2=DefaultHead2();
@@ -282,14 +283,14 @@ arrowhead3 HookHead2(real dir=arrowdir, real barb=arrowbarb, triple normal=O)
filltype filltype=null, bool forwards=true,
projection P=currentprojection) {
if(size == 0) size=a.size(p);
- angle=min(angle*arrowhookfactor,45);
path h=a.project(g,forwards,P);
+ a.size=min(size,arclength(h));
path[] H=a.align(HookHead.head(h,p,size,angle),h);
H=forwards ? yscale(-1)*H : H;
return a.surface(g,position,size,H,p,filltype,normal,P);
};
a.arrowhead2=HookHead;
- a.gap=1.01;
+ a.gap=1.005;
return a;
}
arrowhead3 HookHead2=HookHead2();
@@ -302,6 +303,7 @@ arrowhead3 TeXHead2(triple normal=O) {
bool forwards=true, projection P=currentprojection) {
if(size == 0) size=a.size(p);
path h=a.project(g,forwards,P);
+ a.size=min(size,arclength(h));
h=rotate(-degrees(dir(h,length(h)),warn=false))*h;
path[] H=TeXHead.head(h,p,size,angle);
H=forwards ? yscale(-1)*H : H;
@@ -311,8 +313,8 @@ arrowhead3 TeXHead2(triple normal=O) {
};
a.arrowhead2=TeXHead;
a.size=TeXHead.size;
- a.gap=0.83;
a.splitpath=false;
+ a.gap=1.005;
return a;
}
arrowhead3 TeXHead2=TeXHead2();
@@ -359,6 +361,7 @@ void drawarrow(picture pic, arrowhead3 arrowhead=DefaultHead3,
path3 r=subpath(g,position,0);
size=min(arrowsizelimit*arclength(r),size);
surface head=arrowhead.head(g,position,q,size,angle,filltype,forwards,P);
+ if(arrowhead.size > 0) size=arrowhead.size;
bool endpoint=position > L-sqrtEpsilon;
if(arrowhead.splitpath || endpoint) {
if(position > 0) {
@@ -526,11 +529,11 @@ void bar(picture pic, triple a, triple d, triple perp=O,
pic.add(new void(frame f, transform3 t, picture pic2, projection P) {
picture opic;
triple A=t*a;
- triple v=d == O ? abs(perp)*unit(cross(P.vector(),perp)) : d;
+ triple v=d == O ? abs(perp)*unit(cross(P.normal,perp)) : d;
draw(opic,A-v--A+v,p,light);
add(f,opic.fit3(identity4,pic2,P));
});
- triple v=d == O ? cross(currentprojection.vector(),perp) : d;
+ triple v=d == O ? cross(currentprojection.normal,perp) : d;
pen q=(pen) p;
triple m=min3(q);
triple M=max3(q);
diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy
index e520229371f..1be7eb0259c 100644
--- a/Build/source/utils/asymptote/base/three_surface.asy
+++ b/Build/source/utils/asymptote/base/three_surface.asy
@@ -1343,7 +1343,7 @@ triple rectify(triple dir)
path3[] align(path3[] g, transform3 t=identity4, triple position,
triple align, pen p=currentpen)
{
- if(determinant(t) == 0) return g;
+ if(determinant(t) == 0 || g.length == 0) return g;
triple m=min(g);
triple dir=rectify(inverse(t)*-align);
triple a=m+realmult(dir,max(g)-m);
@@ -1353,7 +1353,7 @@ path3[] align(path3[] g, transform3 t=identity4, triple position,
surface align(surface s, transform3 t=identity4, triple position,
triple align, pen p=currentpen)
{
- if(determinant(t) == 0) return s;
+ if(determinant(t) == 0 || s.s.length == 0) return s;
triple m=min(s);
triple dir=rectify(inverse(t)*-align);
triple a=m+realmult(dir,max(s)-m);
@@ -1367,17 +1367,11 @@ surface surface(Label L, triple position=O)
shift(position)*L.T3*s;
}
-path[] path(Label L, pair z=0, projection P)
+private path[] path(Label L, pair z=0, projection P)
{
path[] g=texpath(L);
- if(L.defaulttransform3) {
- return L.align.is3D ? align(g,z,project(L.align.dir3,P)-project(O,P),L.p) :
- shift(z)*g;
- } else {
- path3[] G=path3(g);
- return L.align.is3D ? shift(z)*project(align(G,L.T3,O,L.align.dir3,L.p),P) :
- shift(z)*project(L.T3*G,P);
- }
+ return L.align.is3D ? align(g,z,project(L.align.dir3,P)-project(O,P),L.p) :
+ shift(z)*g;
}
void label(frame f, Label L, triple position, align align=NoAlign,
@@ -1397,15 +1391,18 @@ void label(frame f, Label L, triple position, align align=NoAlign,
for(patch S : surface(L,position).s)
draw3D(f,S,position,L.p,light,partname(name,++i),interaction);
} else {
- if(L.filltype == NoFill)
- fill(f,path(L,project(position,P.t),P),
- color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
- else {
- frame d;
- fill(d,path(L,project(position,P.t),P),
- color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
- add(f,d,L.filltype);
- }
+ pen p=color(L.T3*Z,L.p,light,shiftless(P.T.modelview));
+ if(L.defaulttransform3) {
+ if(L.filltype == NoFill)
+ fill(f,path(L,project(position,P.t),P),p);
+ else {
+ frame d;
+ fill(d,path(L,project(position,P.t),P),p);
+ add(f,d,L.filltype);
+ }
+ } else
+ for(patch S : surface(L,position).s)
+ fill(f,project(S.external(),P,1),p);
}
}
@@ -1430,22 +1427,30 @@ void label(picture pic=currentpicture, Label L, triple position,
L.T=L.T*scale(abs(P.camera-v)/abs(P.vector()));
if(L.defaulttransform3)
L.T3=transform3(P);
+
if(is3D()) {
int i=-1;
for(patch S : surface(L,v).s)
draw3D(f,S,v,L.p,light,partname(name,++i),interaction);
}
+
if(pic != null) {
- if(L.filltype == NoFill)
- fill(project(v,P.t),pic,path(L,P),
- color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
- else {
- picture d;
- fill(project(v,P.t),d,path(L,P),
- color(L.T3*Z,L.p,light,shiftless(P.T.modelview)));
- add(pic,d,L.filltype);
- }
+ pen p=color(L.T3*Z,L.p,light,shiftless(P.T.modelview));
+ if(L.defaulttransform3) {
+ if(L.filltype == NoFill)
+ fill(project(v,P.t),pic,path(L,P),p);
+ else {
+ picture d;
+ fill(project(v,P.t),d,path(L,P),p);
+ add(pic,d,L.filltype);
+ }
+ } else
+ pic.add(new void(frame f, transform T) {
+ for(patch S : surface(L,v).s)
+ fill(f,T*project(S.external(),P,1),p);
+ });
}
+
},!L.defaulttransform3);
Label L=L.copy();
@@ -1678,12 +1683,28 @@ void dot(picture pic=currentpicture, Label L, triple v, align align=NoAlign,
label(pic,L,v);
}
+pair minbound(triple[] A, projection P)
+{
+ pair b=project(A[0],P);
+ for(triple v : A)
+ b=minbound(b,project(v,P.t));
+ return b;
+}
+
+pair maxbound(triple[] A, projection P)
+{
+ pair b=project(A[0],P);
+ for(triple v : A)
+ b=maxbound(b,project(v,P.t));
+ return b;
+}
+
pair minbound(triple[][] A, projection P)
{
pair b=project(A[0][0],P);
for(triple[] a : A) {
for(triple v : a) {
- b=minbound(b,project(v,P));
+ b=minbound(b,project(v,P.t));
}
}
return b;
@@ -1694,7 +1715,7 @@ pair maxbound(triple[][] A, projection P)
pair b=project(A[0][0],P);
for(triple[] a : A) {
for(triple v : a) {
- b=maxbound(b,project(v,P));
+ b=maxbound(b,project(v,P.t));
}
}
return b;
@@ -1711,6 +1732,24 @@ triple[][] operator / (triple[][] a, real[][] b)
return A;
}
+// Draw a NURBS curve.
+void draw(picture pic=currentpicture, triple[] P, real[] knot,
+ real[] weights=new real[], pen p=currentpen, string name="")
+{
+ P=copy(P);
+ knot=copy(knot);
+ weights=copy(weights);
+ pic.add(new void(frame f, transform3 t, picture pic, projection Q) {
+ if(is3D()) {
+ triple[] P=t*P;
+ draw(f,P,knot,weights,p,name);
+ if(pic != null)
+ pic.addBox(minbound(P,Q),maxbound(P,Q));
+ }
+ },true);
+ pic.addBox(minbound(P),maxbound(P));
+}
+
// Draw a NURBS surface.
void draw(picture pic=currentpicture, triple[][] P, real[] uknot, real[] vknot,
real[][] weights=new real[][], material m=currentpen,
@@ -1734,14 +1773,11 @@ void draw(picture pic=currentpicture, triple[][] P, real[] uknot, real[] vknot,
PRCshininess=PRCshininess(m.shininess);
draw(f,P,uknot,vknot,weights,m.p,m.opacity,m.shininess,PRCshininess,
granularity,colors,lighton,name);
- if(pic != null) {
- triple[][] R=weights.length > 0 ? P/weights : P;
- pic.addBox(minbound(R,Q),maxbound(R,Q));
- }
+ if(pic != null)
+ pic.addBox(minbound(P,Q),maxbound(P,Q));
}
},true);
- triple[][] R=weights.length > 0 ? P/weights : P;
- pic.addBox(minbound(R),maxbound(R));
+ pic.addBox(minbound(P),maxbound(P));
}
// A structure to subdivide two intersecting patches about their intersection.