diff options
author | Karl Berry <karl@freefriends.org> | 2010-04-19 23:59:33 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2010-04-19 23:59:33 +0000 |
commit | dddce3148a2bc785603576b18ffcf72d39adbe47 (patch) | |
tree | 14f4c79ef11bcf0820a8d4ff2ab0c5ac009c365c /Build/source/utils/asymptote/base | |
parent | a01e51b01f5819b6091af48cdca581e9f2a9282e (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')
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. |