From 752012c605d34cd943795527a9738475a6958fcc Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Sun, 7 Apr 2013 18:19:31 +0000 Subject: texmf -> texmf-dist: start with unique dirs from texmf git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751 --- Master/texmf-dist/asymptote/three.asy | 3319 +++++++++++++++++++++++++++++++++ 1 file changed, 3319 insertions(+) create mode 100644 Master/texmf-dist/asymptote/three.asy (limited to 'Master/texmf-dist/asymptote/three.asy') diff --git a/Master/texmf-dist/asymptote/three.asy b/Master/texmf-dist/asymptote/three.asy new file mode 100644 index 00000000000..5d7ae1a8ba5 --- /dev/null +++ b/Master/texmf-dist/asymptote/three.asy @@ -0,0 +1,3319 @@ +private import math; + +if(inXasyMode) settings.render=0; + +if(prc0()) { + if(settings.tex == "context") settings.prc=false; + else { + access embed; + Embed=embed.embed; + Link=embed.link; + } +} + +// Useful lossy compression values. +restricted real Zero=0; +restricted real Low=0.0001; +restricted real Medium=0.001; +restricted real High=0.01; + +restricted int PRCsphere=0; // Renders slowly but produces smaller PRC files. +restricted int NURBSsphere=1; // Renders fast but produces larger PRC files. + +struct render +{ + // PRC parameters: + real compression; // lossy compression parameter (0=no compression) + real granularity; // PRC rendering granularity + + bool closed; // use one-sided rendering? + bool tessellate; // use tessellated mesh to store straight patches? + + bool3 merge; // merge nodes before rendering, for faster but + // lower quality PRC rendering (the value default means + // merge opaque patches only). + + int sphere; // PRC sphere type (PRCsphere or NURBSsphere). + + // General parameters: + real margin; // shrink amount for rendered openGL viewport, in bp. + real tubegranularity; // granularity for rendering tubes + bool labelfill; // fill subdivision cracks in unlighted labels + + static render defaultrender; + + void operator init(real compression=defaultrender.compression, + real granularity=defaultrender.granularity, + bool closed=defaultrender.closed, + bool tessellate=defaultrender.tessellate, + bool3 merge=defaultrender.merge, + int sphere=defaultrender.sphere, + real margin=defaultrender.margin, + real tubegranularity=defaultrender.tubegranularity, + bool labelfill=defaultrender.labelfill) + { + this.compression=compression; + this.granularity=granularity; + this.closed=closed; + this.tessellate=tessellate; + this.merge=merge; + this.sphere=sphere; + this.margin=margin; + this.tubegranularity=tubegranularity; + this.labelfill=labelfill; + } +} + +render operator init() {return render();} + +render defaultrender=render.defaultrender=new render; +defaultrender.compression=High; +defaultrender.granularity=Medium; +defaultrender.closed=false; +defaultrender.tessellate=false; +defaultrender.merge=false; +defaultrender.margin=0.02; +defaultrender.tubegranularity=0.005; +defaultrender.sphere=NURBSsphere; +defaultrender.labelfill=true; + +real defaultshininess=0.25; + +real angleprecision=1e-5; // Precision for centering perspective projections. +int maxangleiterations=25; + +string defaultembed3Doptions; +string defaultembed3Dscript; +real defaulteyetoview=63mm/1000mm; + +string partname(int i=0) +{ + return string(i+1); +} + +triple O=(0,0,0); +triple X=(1,0,0), Y=(0,1,0), Z=(0,0,1); + +// A translation in 3D space. +transform3 shift(triple v) +{ + transform3 t=identity(4); + t[0][3]=v.x; + t[1][3]=v.y; + t[2][3]=v.z; + return t; +} + +// Avoid two parentheses. +transform3 shift(real x, real y, real z) +{ + return shift((x,y,z)); +} + +transform3 shift(transform3 t) +{ + transform3 T=identity(4); + T[0][3]=t[0][3]; + T[1][3]=t[1][3]; + T[2][3]=t[2][3]; + return T; +} + +// A 3D scaling in the x direction. +transform3 xscale3(real x) +{ + transform3 t=identity(4); + t[0][0]=x; + return t; +} + +// A 3D scaling in the y direction. +transform3 yscale3(real y) +{ + transform3 t=identity(4); + t[1][1]=y; + return t; +} + +// A 3D scaling in the z direction. +transform3 zscale3(real z) +{ + transform3 t=identity(4); + t[2][2]=z; + return t; +} + +// A 3D scaling by s in the v direction. +transform3 scale(triple v, real s) +{ + v=unit(v); + s -= 1; + return new real[][] { + {1+s*v.x^2, s*v.x*v.y, s*v.x*v.z, 0}, + {s*v.x*v.y, 1+s*v.y^2, s*v.y*v.z, 0}, + {s*v.x*v.z, s*v.y*v.z, 1+s*v.z^2, 0}, + {0, 0, 0, 1}}; +} + +// A transformation representing rotation by an angle in degrees about +// an axis v through the origin (in the right-handed direction). +transform3 rotate(real angle, triple v) +{ + if(v == O) abort("cannot rotate about the zero vector"); + v=unit(v); + real x=v.x, y=v.y, z=v.z; + real s=Sin(angle), c=Cos(angle), t=1-c; + + return new real[][] { + {t*x^2+c, t*x*y-s*z, t*x*z+s*y, 0}, + {t*x*y+s*z, t*y^2+c, t*y*z-s*x, 0}, + {t*x*z-s*y, t*y*z+s*x, t*z^2+c, 0}, + {0, 0, 0, 1}}; +} + +// A transformation representing rotation by an angle in degrees about +// the line u--v (in the right-handed direction). +transform3 rotate(real angle, triple u, triple v) +{ + return shift(u)*rotate(angle,v-u)*shift(-u); +} + +// Reflects about the plane through u, v, and w. +transform3 reflect(triple u, triple v, triple w) +{ + triple n=unit(cross(v-u,w-u)); + if(n == O) + abort("points determining reflection plane cannot be colinear"); + + return new real[][] { + {1-2*n.x^2, -2*n.x*n.y, -2*n.x*n.z, u.x}, + {-2*n.x*n.y, 1-2*n.y^2, -2*n.y*n.z, u.y}, + {-2*n.x*n.z, -2*n.y*n.z, 1-2*n.z^2, u.z}, + {0, 0, 0, 1} + }*shift(-u); +} + +// Project u onto v. +triple project(triple u, triple v) +{ + v=unit(v); + return dot(u,v)*v; +} + +// Return a unit vector perpendicular to a given unit vector v. +triple perp(triple v) +{ + triple u=cross(v,Y); + if(abs(u) > sqrtEpsilon) return unit(u); + u=cross(v,Z); + return (abs(u) > sqrtEpsilon) ? unit(u) : X; +} + +// Return the transformation corresponding to moving the camera from the target +// (looking in the negative z direction) to the point 'eye' (looking at target), +// orienting the camera so that direction 'up' points upwards. +// Since, in actuality, we are transforming the points instead of the camera, +// we calculate the inverse matrix. +// Based on the gluLookAt implementation in the OpenGL manual. +transform3 look(triple eye, triple up=Z, triple target=O) +{ + triple f=unit(target-eye); + if(f == O) + f=-Z; // The eye is already at the origin: look down. + + triple s=cross(f,up); + + // If the eye is pointing either directly up or down, there is no + // preferred "up" direction. Pick one arbitrarily. + s=s != O ? unit(s) : perp(f); + + triple u=cross(s,f); + + transform3 M={{ s.x, s.y, s.z, 0}, + { u.x, u.y, u.z, 0}, + {-f.x, -f.y, -f.z, 0}, + { 0, 0, 0, 1}}; + + return M*shift(-eye); +} + +// Return a matrix to do perspective distortion based on a triple v. +transform3 distort(triple v) +{ + transform3 t=identity(4); + real d=length(v); + if(d == 0) return t; + t[3][2]=-1/d; + t[3][3]=0; + return t; +} + +projection operator * (transform3 t, projection P) +{ + projection P=P.copy(); + if(!P.absolute) { + P.camera=t*P.camera; + triple target=P.target; + P.target=t*P.target; + if(P.infinity) + P.normal=t*(target+P.normal)-P.target; + else + P.normal=P.vector(); + P.calculate(); + } + return P; +} + +// With this, save() and restore() in plain also save and restore the +// currentprojection. +addSaveFunction(new restoreThunk() { + projection P=currentprojection.copy(); + return new void() { + currentprojection=P; + }; + }); + +pair project(triple v, projection P=currentprojection) +{ + return project(v,P.t); +} + +pair dir(triple v, triple dir, projection P) +{ + return unit(project(v+0.5dir,P)-project(v-0.5*dir,P)); +} + +// Uses the homogenous coordinate to perform perspective distortion. +// When combined with a projection to the XY plane, this effectively maps +// points in three space to a plane through target and +// perpendicular to the vector camera-target. +projection perspective(triple camera, triple up=Z, triple target=O, + real zoom=1, real angle=0, pair viewportshift=0, + bool showtarget=true, bool autoadjust=true, + bool center=autoadjust) +{ + if(camera == target) + abort("camera cannot be at target"); + return projection(camera,up,target,zoom,angle,viewportshift, + showtarget,autoadjust,center, + new transformation(triple camera, triple up, triple target) + {return transformation(look(camera,up,target), + distort(camera-target));}); +} + +projection perspective(real x, real y, real z, triple up=Z, triple target=O, + real zoom=1, real angle=0, pair viewportshift=0, + bool showtarget=true, bool autoadjust=true, + bool center=autoadjust) +{ + return perspective((x,y,z),up,target,zoom,angle,viewportshift,showtarget, + autoadjust,center); +} + +projection orthographic(triple camera, triple up=Z, triple target=O, + real zoom=1, pair viewportshift=0, + bool showtarget=true, bool center=false) +{ + return projection(camera,up,target,zoom,viewportshift,showtarget, + center=center,new transformation(triple camera, triple up, + triple target) { + return transformation(look(camera,up,target));}); +} + +projection orthographic(real x, real y, real z, triple up=Z, + triple target=O, real zoom=1, pair viewportshift=0, + bool showtarget=true, bool center=false) +{ + return orthographic((x,y,z),up,target,zoom,viewportshift,showtarget, + center=center); +} + +// Compute camera position with x axis below the horizontal at angle alpha, +// y axis below the horizontal at angle beta, and z axis up. +triple camera(real alpha, real beta) +{ + real denom=Tan(alpha+beta); + real Tanalpha=Tan(alpha); + real Tanbeta=Tan(beta); + return (sqrt(Tanalpha/denom),sqrt(Tanbeta/denom),sqrt(Tanalpha*Tanbeta)); +} + +projection oblique(real angle=45) +{ + transform3 t=identity(4); + real c2=Cos(angle)^2; + real s2=1-c2; + t[0][2]=-c2; + t[1][2]=-s2; + t[2][2]=1; + t[2][3]=-1; + return projection((c2,s2,1),up=Y,normal=(0,0,1), + new transformation(triple,triple,triple) { + return transformation(t);}); +} + +projection obliqueZ(real angle=45) {return oblique(angle);} + +projection obliqueX(real angle=45) +{ + transform3 t=identity(4); + real c2=Cos(angle)^2; + real s2=1-c2; + t[0][0]=-c2; + t[1][0]=-s2; + t[1][1]=0; + t[0][1]=1; + t[1][2]=1; + t[2][2]=0; + t[2][0]=1; + t[2][3]=-1; + return projection((1,c2,s2),normal=(1,0,0), + new transformation(triple,triple,triple) { + return transformation(t);}); +} + +projection obliqueY(real angle=45) +{ + transform3 t=identity(4); + real c2=Cos(angle)^2; + real s2=1-c2; + t[0][1]=c2; + t[1][1]=s2; + t[1][2]=1; + t[2][1]=-1; + t[2][2]=0; + t[2][3]=-1; + return projection((c2,-1,s2),normal=(0,-1,0), + new transformation(triple,triple,triple) { + return transformation(t);}); +} + +projection oblique=oblique(); +projection obliqueX=obliqueX(), obliqueY=obliqueY(), obliqueZ=obliqueZ(); + +projection LeftView=orthographic(-X,showtarget=true); +projection RightView=orthographic(X,showtarget=true); +projection FrontView=orthographic(-Y,showtarget=true); +projection BackView=orthographic(Y,showtarget=true); +projection BottomView=orthographic(-Z,up=-Y,showtarget=true); +projection TopView=orthographic(Z,up=Y,showtarget=true); + +currentprojection=perspective(5,4,2); + +projection projection() +{ + projection P; + real[] a=_projection(); + if(a.length == 0 || a[10] == 0.0) return currentprojection; + int k=0; + return a[0] == 1 ? + orthographic((a[++k],a[++k],a[++k]),(a[++k],a[++k],a[++k]), + (a[++k],a[++k],a[++k]),a[++k],(a[k += 2],a[++k])) : + perspective((a[++k],a[++k],a[++k]),(a[++k],a[++k],a[++k]), + (a[++k],a[++k],a[++k]),a[++k],a[++k],(a[++k],a[++k])); +} + +// Map pair z to a triple by inverting the projection P onto the +// plane perpendicular to normal and passing through point. +triple invert(pair z, triple normal, triple point, + projection P=currentprojection) +{ + transform3 t=P.t; + real[][] A={{t[0][0]-z.x*t[3][0],t[0][1]-z.x*t[3][1],t[0][2]-z.x*t[3][2]}, + {t[1][0]-z.y*t[3][0],t[1][1]-z.y*t[3][1],t[1][2]-z.y*t[3][2]}, + {normal.x,normal.y,normal.z}}; + real[] b={z.x*t[3][3]-t[0][3],z.y*t[3][3]-t[1][3],dot(normal,point)}; + real[] x=solve(A,b,warn=false); + return x.length > 0 ? (x[0],x[1],x[2]) : P.camera; +} + +// Map pair to a triple on the projection plane. +triple invert(pair z, projection P=currentprojection) +{ + 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.normal,v,P)-v; +} + +pair xypart(triple v) +{ + return (v.x,v.y); +} + +struct control { + triple post,pre; + bool active=false; + bool straight=true; + void operator init(triple post, triple pre, bool straight=false) { + this.post=post; + this.pre=pre; + active=true; + this.straight=straight; + } +} + +control nocontrol; + +control operator * (transform3 t, control c) +{ + control C; + C.post=t*c.post; + C.pre=t*c.pre; + C.active=c.active; + C.straight=c.straight; + return C; +} + +void write(file file, control c) +{ + write(file,".. controls "); + write(file,c.post); + write(file," and "); + write(file,c.pre); +} + +struct Tension { + real out,in; + bool atLeast; + bool active; + void operator init(real out=1, real in=1, bool atLeast=false, + bool active=true) { + real check(real val) { + if(val < 0.75) abort("tension cannot be less than 3/4"); + return val; + } + this.out=check(out); + this.in=check(in); + this.atLeast=atLeast; + this.active=active; + } +} + +Tension operator init() +{ + return Tension(); +} + +Tension noTension; +noTension.active=false; + +void write(file file, Tension t) +{ + write(file,"..tension "); + if(t.atLeast) write(file,"atleast "); + write(file,t.out); + write(file," and "); + write(file,t.in); +} + +struct dir { + triple dir; + real gamma=1; // endpoint curl + bool Curl; // curl specified + bool active() { + return dir != O || Curl; + } + void init(triple v) { + this.dir=v; + } + void init(real gamma) { + if(gamma < 0) abort("curl cannot be less than 0"); + this.gamma=gamma; + this.Curl=true; + } + void init(dir d) { + dir=d.dir; + gamma=d.gamma; + Curl=d.Curl; + } + void default(triple v) { + if(!active()) init(v); + } + void default(dir d) { + if(!active()) init(d); + } + dir copy() { + dir d=new dir; + d.init(this); + return d; + } +} + +void write(file file, dir d) +{ + if(d.dir != O) { + write(file,"{"); write(file,unit(d.dir)); write(file,"}"); + } else if(d.Curl) { + write(file,"{curl "); write(file,d.gamma); write(file,"}"); + } +} + +dir operator * (transform3 t, dir d) +{ + dir D=d.copy(); + D.init(unit(shiftless(t)*d.dir)); + return D; +} + +void checkEmpty(int n) { + if(n == 0) + abort("nullpath3 has no points"); +} + +int adjustedIndex(int i, int n, bool cycles) +{ + checkEmpty(n); + if(cycles) + return i % n; + else if(i < 0) + return 0; + else if(i >= n) + return n-1; + else + return i; +} + +struct flatguide3 { + triple[] nodes; + bool[] cyclic; // true if node is really a cycle + control[] control; // control points for segment starting at node + Tension[] Tension; // Tension parameters for segment starting at node + dir[] in,out; // in and out directions for segment starting at node + + bool cyclic() {int n=cyclic.length; return n > 0 ? cyclic[n-1] : false;} + bool precyclic() {int i=find(cyclic); return i >= 0 && i < cyclic.length-1;} + + int size() { + return cyclic() ? nodes.length-1 : nodes.length; + } + + void node(triple v, bool b=false) { + nodes.push(v); + control.push(nocontrol); + Tension.push(noTension); + in.push(new dir); + out.push(new dir); + cyclic.push(b); + } + + void control(triple post, triple pre) { + if(control.length > 0) { + control c=control(post,pre,false); + control[control.length-1]=c; + } + } + + void Tension(real out, real in, bool atLeast) { + if(Tension.length > 0) + Tension[Tension.length-1]=Tension(out,in,atLeast,true); + } + + void in(triple v) { + if(in.length > 0) { + in[in.length-1].init(v); + } + } + + void out(triple v) { + if(out.length > 0) { + out[out.length-1].init(v); + } + } + + void in(real gamma) { + if(in.length > 0) { + in[in.length-1].init(gamma); + } + } + + void out(real gamma) { + if(out.length > 0) { + out[out.length-1].init(gamma); + } + } + + void cycleToken() { + if(nodes.length > 0) + node(nodes[0],true); + } + + // Return true if outgoing direction at node i is known. + bool solved(int i) { + return out[i].active() || control[i].active; + } +} + +void write(file file, string s="", explicit flatguide3 x, suffix suffix=none) +{ + write(file,s); + if(x.size() == 0) write(file,""); + else for(int i=0; i < x.nodes.length; ++i) { + if(i > 0) write(file,endl); + if(x.cyclic[i]) write(file,"cycle"); + else write(file,x.nodes[i]); + if(i < x.nodes.length-1) { + // Explicit control points trump other specifiers + if(x.control[i].active) + write(file,x.control[i]); + else { + write(file,x.out[i]); + if(x.Tension[i].active) write(file,x.Tension[i]); + } + write(file,".."); + if(!x.control[i].active) write(file,x.in[i]); + } + } + write(file,suffix); +} + +void write(string s="", flatguide3 x, suffix suffix=endl) +{ + write(stdout,s,x,suffix); +} + +// A guide3 is most easily represented as something that modifies a flatguide3. +typedef void guide3(flatguide3); + +restricted void nullpath3(flatguide3) {}; + +guide3 operator init() {return nullpath3;} + +guide3 operator cast(triple v) +{ + return new void(flatguide3 f) { + f.node(v); + }; +} + +guide3 operator cast(cycleToken) { + return new void(flatguide3 f) { + f.cycleToken(); + }; +} + +guide3 operator controls(triple post, triple pre) +{ + return new void(flatguide3 f) { + f.control(post,pre); + }; +}; + +guide3 operator controls(triple v) +{ + return operator controls(v,v); +} + +guide3 operator cast(tensionSpecifier t) +{ + return new void(flatguide3 f) { + f.Tension(t.out, t.in, t.atLeast); + }; +} + +guide3 operator cast(curlSpecifier spec) +{ + return new void(flatguide3 f) { + if(spec.side == JOIN_OUT) f.out(spec.value); + else if(spec.side == JOIN_IN) f.in(spec.value); + else + abort("invalid curl specifier"); + }; +} + +guide3 operator spec(triple v, int side) +{ + return new void(flatguide3 f) { + if(side == JOIN_OUT) f.out(v); + else if(side == JOIN_IN) f.in(v); + else + abort("invalid direction specifier"); + }; +} + +guide3 operator -- (... guide3[] g) +{ + return new void(flatguide3 f) { + if(g.length > 0) { + for(int i=0; i < g.length-1; ++i) { + g[i](f); + f.out(1); + f.in(1); + } + g[g.length-1](f); + } + }; +} + +guide3 operator .. (... guide3[] g) +{ + return new void(flatguide3 f) { + for(int i=0; i < g.length; ++i) + g[i](f); + }; +} + +typedef guide3 interpolate3(... guide3[]); + +interpolate3 join3(tensionSpecifier t) +{ + 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; + g(f); + return f; +} + +flatguide3[] operator cast(guide3[] g) +{ + flatguide3[] p=new flatguide3[g.length]; + for(int i=0; i < g.length; ++i) { + flatguide3 f; + g[i](f); + p[i]=f; + } + return p; +} + +// A version of asin that tolerates numerical imprecision +real asin1(real x) +{ + return asin(min(max(x,-1),1)); +} + +// A version of acos that tolerates numerical imprecision +real acos1(real x) +{ + return acos(min(max(x,-1),1)); +} + +struct Controls { + triple c0,c1; + + // 3D extension of John Hobby's control point formula + // (cf. The MetaFont Book, page 131), + // as described in John C. Bowman and A. Hammerlindl, + // TUGBOAT: The Communications of the TeX Users Group 29:2 (2008). + + void operator init(triple v0, triple v1, triple d0, triple d1, real tout, + real tin, bool atLeast) { + triple v=v1-v0; + triple u=unit(v); + real L=length(v); + d0=unit(d0); + d1=unit(d1); + real theta=acos1(dot(d0,u)); + real phi=acos1(dot(d1,u)); + if(dot(cross(d0,v),cross(v,d1)) < 0) phi=-phi; + c0=v0+d0*L*relativedistance(theta,phi,tout,atLeast); + c1=v1-d1*L*relativedistance(phi,theta,tin,atLeast); + } +} + +private triple cross(triple d0, triple d1, triple reference) +{ + triple normal=cross(d0,d1); + return normal == O ? reference : normal; +} + +private triple dir(real theta, triple d0, triple d1, triple reference) +{ + triple normal=cross(d0,d1,reference); + if(normal == O) return d1; + return rotate(degrees(theta),dot(normal,reference) >= 0 ? normal : -normal)* + d1; +} + +private real angle(triple d0, triple d1, triple reference) +{ + real theta=acos1(dot(unit(d0),unit(d1))); + return dot(cross(d0,d1,reference),reference) >= 0 ? theta : -theta; +} + +// 3D extension of John Hobby's angle formula (The MetaFont Book, page 131). +// Notational differences: here psi[i] is the turning angle at z[i+1], +// beta[i] is the tension for segment i, and in[i] is the incoming +// direction for segment i (where segment i begins at node i). + +real[] theta(triple[] v, real[] alpha, real[] beta, + triple dir0, triple dirn, real g0, real gn, triple reference) +{ + real[] a,b,c,f,l,psi; + int n=alpha.length; + bool cyclic=v.cyclic; + for(int i=0; i < n; ++i) + l[i]=1/length(v[i+1]-v[i]); + int i0,in; + if(cyclic) {i0=0; in=n;} + else {i0=1; in=n-1;} + for(int i=0; i < in; ++i) + psi[i]=angle(v[i+1]-v[i],v[i+2]-v[i+1],reference); + if(cyclic) { + l.cyclic=true; + psi.cyclic=true; + } else { + psi[n-1]=0; + if(dir0 == O) { + real a0=alpha[0]; + real b0=beta[0]; + real chi=g0*(b0/a0)^2; + a[0]=0; + b[0]=3a0-a0/b0+chi; + real C=chi*(3a0-1)+a0/b0; + c[0]=C; + f[0]=-C*psi[0]; + } else { + a[0]=c[0]=0; + b[0]=1; + f[0]=angle(v[1]-v[0],dir0,reference); + } + if(dirn == O) { + real an=alpha[n-1]; + real bn=beta[n-1]; + real chi=gn*(an/bn)^2; + a[n]=chi*(3bn-1)+bn/an; + b[n]=3bn-bn/an+chi; + c[n]=f[n]=0; + } else { + a[n]=c[n]=0; + b[n]=1; + f[n]=angle(v[n]-v[n-1],dirn,reference); + } + } + + for(int i=i0; i < n; ++i) { + real in=beta[i-1]^2*l[i-1]; + real A=in/alpha[i-1]; + a[i]=A; + real B=3*in-A; + real out=alpha[i]^2*l[i]; + real C=out/beta[i]; + b[i]=B+3*out-C; + c[i]=C; + f[i]=-B*psi[i-1]-C*psi[i]; + } + + return tridiagonal(a,b,c,f); +} + +triple reference(triple[] v, int n, triple d0, triple d1) +{ + triple[] V=sequence(new triple(int i) { + return cross(v[i+1]-v[i],v[i+2]-v[i+1]); + },n-1); + if(n > 0) { + V.push(cross(d0,v[1]-v[0])); + V.push(cross(v[n]-v[n-1],d1)); + } + + triple max=V[0]; + real M=abs(max); + for(int i=1; i < V.length; ++i) { + triple vi=V[i]; + real a=abs(vi); + if(a > M) { + M=a; + max=vi; + } + } + + triple reference; + for(int i=0; i < V.length; ++i) { + triple u=unit(V[i]); + reference += dot(u,max) < 0 ? -u : u; + } + + return reference; +} + +// Fill in missing directions for n cyclic nodes. +void aim(flatguide3 g, int N) +{ + bool cyclic=true; + int start=0, end=0; + + // If the cycle contains one or more direction specifiers, break the loop. + for(int k=0; k < N; ++k) + if(g.solved(k)) {cyclic=false; end=k; break;} + for(int k=N-1; k >= 0; --k) + if(g.solved(k)) {cyclic=false; start=k; break;} + while(start < N && g.control[start].active) ++start; + + int n=N-(start-end); + if(n <= 1 || (cyclic && n <= 2)) return; + + triple[] v=new triple[cyclic ? n : n+1]; + real[] alpha=new real[n]; + real[] beta=new real[n]; + for(int k=0; k < n; ++k) { + int K=(start+k) % N; + v[k]=g.nodes[K]; + alpha[k]=g.Tension[K].out; + beta[k]=g.Tension[K].in; + } + if(cyclic) { + v.cyclic=true; + alpha.cyclic=true; + beta.cyclic=true; + } else v[n]=g.nodes[(start+n) % N]; + int final=(end-1) % N; + + triple d0=g.out[start].dir; + triple d1=g.in[final].dir; + + triple reference=reference(v,n,d0,d1); + + real[] theta=theta(v,alpha,beta,d0,d1,g.out[start].gamma,g.in[final].gamma, + reference); + + v.cyclic=true; + theta.cyclic=true; + + for(int k=1; k < (cyclic ? n+1 : n); ++k) { + triple w=dir(theta[k],v[k]-v[k-1],v[k+1]-v[k],reference); + g.in[(start+k-1) % N].init(w); + g.out[(start+k) % N].init(w); + } + + if(g.out[start].dir == O) + g.out[start].init(dir(theta[0],v[0]-g.nodes[(start-1) % N],v[1]-v[0], + reference)); + if(g.in[final].dir == O) + g.in[final].init(dir(theta[n],v[n-1]-v[n-2],v[n]-v[n-1],reference)); +} + +// Fill in missing directions for the sequence of nodes i...n. +void aim(flatguide3 g, int i, int n) +{ + int j=n-i; + if(j > 1 || g.out[i].dir != O || g.in[i].dir != O) { + triple[] v=new triple[j+1]; + real[] alpha=new real[j]; + real[] beta=new real[j]; + for(int k=0; k < j; ++k) { + v[k]=g.nodes[i+k]; + alpha[k]=g.Tension[i+k].out; + beta[k]=g.Tension[i+k].in; + } + v[j]=g.nodes[n]; + + triple d0=g.out[i].dir; + triple d1=g.in[n-1].dir; + + triple reference=reference(v,j,d0,d1); + + real[] theta=theta(v,alpha,beta,d0,d1,g.out[i].gamma,g.in[n-1].gamma, + reference); + + for(int k=1; k < j; ++k) { + triple w=dir(theta[k],v[k]-v[k-1],v[k+1]-v[k],reference); + g.in[i+k-1].init(w); + g.out[i+k].init(w); + } + if(g.out[i].dir == O) { + triple w=dir(theta[0],g.in[i].dir,v[1]-v[0],reference); + if(i > 0) g.in[i-1].init(w); + g.out[i].init(w); + } + if(g.in[n-1].dir == O) { + triple w=dir(theta[j],g.out[n-1].dir,v[j]-v[j-1],reference); + g.in[n-1].init(w); + g.out[n].init(w); + } + } +} + +private real Fuzz=10*realEpsilon; + +triple XYplane(pair z) {return (z.x,z.y,0);} +triple YZplane(pair z) {return (0,z.x,z.y);} +triple ZXplane(pair z) {return (z.y,0,z.x);} + +bool cyclic(guide3 g) {flatguide3 f; g(f); return f.cyclic();} +int size(guide3 g) {flatguide3 f; g(f); return f.size();} +int length(guide3 g) {flatguide3 f; g(f); return f.nodes.length-1;} + +triple dir(path3 p) +{ + return dir(p,length(p)); +} + +triple dir(path3 p, path3 h) +{ + return 0.5*(dir(p)+dir(h)); +} + +// return the point on path3 p at arclength L +triple arcpoint(path3 p, real L) +{ + return point(p,arctime(p,L)); +} + +// return the direction on path3 p at arclength L +triple arcdir(path3 p, real L) +{ + return dir(p,arctime(p,L)); +} + +// return the time on path3 p at the relative fraction l of its arclength +real reltime(path3 p, real l) +{ + return arctime(p,l*arclength(p)); +} + +// return the point on path3 p at the relative fraction l of its arclength +triple relpoint(path3 p, real l) +{ + return point(p,reltime(p,l)); +} + +// return the direction of path3 p at the relative fraction l of its arclength +triple reldir(path3 p, real l) +{ + return dir(p,reltime(p,l)); +} + +// return the initial point of path3 p +triple beginpoint(path3 p) +{ + return point(p,0); +} + +// return the point on path3 p at half of its arclength +triple midpoint(path3 p) +{ + return relpoint(p,0.5); +} + +// return the final point of path3 p +triple endpoint(path3 p) +{ + return point(p,length(p)); +} + +path3 path3(triple v) +{ + triple[] point={v}; + return path3(point,point,point,new bool[] {false},false); +} + +path3 path3(path p, triple plane(pair)=XYplane) +{ + int n=size(p); + return path3(sequence(new triple(int i) {return plane(precontrol(p,i));},n), + sequence(new triple(int i) {return plane(point(p,i));},n), + sequence(new triple(int i) {return plane(postcontrol(p,i));},n), + sequence(new bool(int i) {return straight(p,i);},n), + cyclic(p)); +} + +path3[] path3(explicit path[] g, triple plane(pair)=XYplane) +{ + return sequence(new path3(int i) {return path3(g[i],plane);},g.length); +} + +path3 invert(path p, triple normal, triple point, + projection P=currentprojection) +{ + return path3(p,new triple(pair z) {return invert(z,normal,point,P);}); +} + +path3 invert(path p, triple point, projection P=currentprojection) +{ + 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.normal,P.target,P);}); +} + +// Construct a path from a path3 by applying P to each control point. +path path(path3 p, pair P(triple)=xypart) +{ + int n=length(p); + if(n < 0) return nullpath; + guide g=P(point(p,0)); + if(n == 0) return g; + for(int i=1; i < n; ++i) + g=straight(p,i-1) ? g--P(point(p,i)) : + g..controls P(postcontrol(p,i-1)) and P(precontrol(p,i))..P(point(p,i)); + + if(straight(p,n-1)) + return cyclic(p) ? g--cycle : g--P(point(p,n)); + + pair post=P(postcontrol(p,n-1)); + pair pre=P(precontrol(p,n)); + return cyclic(p) ? g..controls post and pre..cycle : + g..controls post and pre..P(point(p,n)); +} + +void write(file file, string s="", explicit path3 x, suffix suffix=none) +{ + write(file,s); + int n=length(x); + if(n < 0) write(""); + else { + for(int i=0; i < n; ++i) { + write(file,point(x,i)); + if(i < length(x)) { + if(straight(x,i)) write(file,"--"); + else { + write(file,".. controls "); + write(file,postcontrol(x,i)); + write(file," and "); + write(file,precontrol(x,i+1),newl); + write(file," .."); + } + } + } + if(cyclic(x)) + write(file,"cycle",suffix); + else + write(file,point(x,n),suffix); + } +} + +void write(string s="", explicit path3 x, suffix suffix=endl) +{ + write(stdout,s,x,suffix); +} + +void write(file file, string s="", explicit path3[] x, suffix suffix=none) +{ + 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 path3[] x, suffix suffix=endl) +{ + write(stdout,s,x,suffix); +} + +path3 solve(flatguide3 g) +{ + int n=g.nodes.length-1; + + // If duplicate points occur consecutively, add dummy controls (if absent). + for(int i=0; i < n; ++i) { + if(g.nodes[i] == g.nodes[i+1] && !g.control[i].active) + g.control[i]=control(g.nodes[i],g.nodes[i],straight=true); + } + + // Fill in empty direction specifiers inherited from explicit control points. + for(int i=0; i < n; ++i) { + if(g.control[i].active) { + g.out[i].init(g.control[i].post-g.nodes[i]); + g.in[i].init(g.nodes[i+1]-g.control[i].pre); + } + } + + // Propagate directions across nodes. + for(int i=0; i < n; ++i) { + int next=g.cyclic[i+1] ? 0 : i+1; + if(g.out[next].active()) + g.in[i].default(g.out[next]); + if(g.in[i].active()) { + g.out[next].default(g.in[i]); + g.out[i+1].default(g.in[i]); + } + } + + // Compute missing 3D directions. + // First, resolve cycles + int i=find(g.cyclic); + if(i > 0) { + aim(g,i); + // All other cycles can now be reduced to sequences. + triple v=g.out[0].dir; + for(int j=i; j <= n; ++j) { + if(g.cyclic[j]) { + g.in[j-1].default(v); + g.out[j].default(v); + if(g.nodes[j-1] == g.nodes[j] && !g.control[j-1].active) + g.control[j-1]=control(g.nodes[j-1],g.nodes[j-1]); + } + } + } + + // Next, resolve sequences. + int i=0; + int start=0; + while(i < n) { + // Look for a missing outgoing direction. + while(i <= n && g.solved(i)) {start=i; ++i;} + if(i > n) break; + // Look for the end of the sequence. + while(i < n && !g.solved(i)) ++i; + + while(start < i && g.control[start].active) ++start; + + if(start < i) + aim(g,start,i); + } + + // Compute missing 3D control points. + for(int i=0; i < n; ++i) { + int next=g.cyclic[i+1] ? 0 : i+1; + if(!g.control[i].active) { + control c; + if((g.out[i].Curl && g.in[i].Curl) || + (g.out[i].dir == O && g.in[i].dir == O)) { + // fill in straight control points for path3 functions + triple delta=(g.nodes[i+1]-g.nodes[i])/3; + c=control(g.nodes[i]+delta,g.nodes[i+1]-delta,straight=true); + } else { + Controls C=Controls(g.nodes[i],g.nodes[next],g.out[i].dir,g.in[i].dir, + g.Tension[i].out,g.Tension[i].in, + g.Tension[i].atLeast); + c=control(C.c0,C.c1); + } + g.control[i]=c; + } + } + + // Convert to Knuth's format (control points stored with nodes) + int n=g.nodes.length; + bool cyclic; + if(n > 0) { + cyclic=g.cyclic[n-1]; + if(cyclic) --n; + } + triple[] pre=new triple[n]; + triple[] point=new triple[n]; + triple[] post=new triple[n]; + bool[] straight=new bool[n]; + if(n > 0) { + for(int i=0; i < n-1; ++i) { + point[i]=g.nodes[i]; + post[i]=g.control[i].post; + pre[i+1]=g.control[i].pre; + straight[i]=g.control[i].straight; + } + point[n-1]=g.nodes[n-1]; + if(cyclic) { + pre[0]=g.control[n-1].pre; + post[n-1]=g.control[n-1].post; + straight[n-1]=g.control[n-1].straight; + } else { + pre[0]=point[0]; + post[n-1]=point[n-1]; + straight[n-1]=false; + } + } + + return path3(pre,point,post,straight,cyclic); +} + +path nurb(path3 p, projection P, int ninterpolate=P.ninterpolate) +{ + triple f=P.camera; + triple u=unit(P.normal); + transform3 t=P.t; + + path nurb(triple v0, triple v1, triple v2, triple v3) { + return nurb(project(v0,t),project(v1,t),project(v2,t),project(v3,t), + dot(u,f-v0),dot(u,f-v1),dot(u,f-v2),dot(u,f-v3),ninterpolate); + } + + path g; + + if(straight(p,0)) + g=project(point(p,0),t); + + int last=length(p); + for(int i=0; i < last; ++i) { + if(straight(p,i)) + g=g--project(point(p,i+1),t); + else + g=g&nurb(point(p,i),postcontrol(p,i),precontrol(p,i+1),point(p,i+1)); + } + + int n=length(g); + if(cyclic(p)) g=g&cycle; + + return g; +} + +path project(path3 p, projection P=currentprojection, + int ninterpolate=P.ninterpolate) +{ + guide g; + + int last=length(p); + if(last < 0) return g; + + transform3 t=P.t; + + if(ninterpolate == 1 || piecewisestraight(p)) { + g=project(point(p,0),t); + // Construct the path. + int stop=cyclic(p) ? last-1 : last; + for(int i=0; i < stop; ++i) { + if(straight(p,i)) + g=g--project(point(p,i+1),t); + else { + g=g..controls project(postcontrol(p,i),t) and + project(precontrol(p,i+1),t)..project(point(p,i+1),t); + } + } + } else return nurb(p,P); + + if(cyclic(p)) + g=straight(p,last-1) ? g--cycle : + g..controls project(postcontrol(p,last-1),t) and + project(precontrol(p,last),t)..cycle; + return g; +} + +pair[] project(triple[] v, projection P=currentprojection) +{ + return sequence(new pair(int i) {return project(v[i],P.t);},v.length); +} + +path[] project(explicit path3[] g, projection P=currentprojection) +{ + return sequence(new path(int i) {return project(g[i],P);},g.length); +} + +guide3 operator cast(path3 p) +{ + int last=length(p); + + bool cyclic=cyclic(p); + int stop=cyclic ? last-1 : last; + return new void(flatguide3 f) { + if(last >= 0) { + f.node(point(p,0)); + for(int i=0; i < stop; ++i) { + if(straight(p,i)) { + f.out(1); + f.in(1); + } else + f.control(postcontrol(p,i),precontrol(p,i+1)); + f.node(point(p,i+1)); + } + if(cyclic) { + if(straight(p,stop)) { + f.out(1); + f.in(1); + } else + f.control(postcontrol(p,stop),precontrol(p,last)); + f.cycleToken(); + } + } + }; +} + +// Return a unit normal vector to a planar path p (or O if the path is +// nonplanar). +triple normal(path3 p) +{ + triple normal; + real fuzz=sqrtEpsilon*abs(max(p)-min(p)); + real absnormal; + real theta; + + bool Cross(triple a, triple b) { + 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 && + abs(normal-n) > sqrtEpsilon && abs(normal+n) > sqrtEpsilon) + return true; + else { + int sign=dot(n,normal) >= 0 ? 1 : -1; + theta += sign*asin1(absn); + if(absn > absnormal) { + absnormal=absn; + normal=n; + theta=sign*theta; + } + } + } + return false; + } + + int L=length(p); + if(L <= 0) return O; + + triple zi=point(p,0); + triple v0=zi-precontrol(p,0); + for(int i=0; i < L; ++i) { + triple c0=postcontrol(p,i); + triple c1=precontrol(p,i+1); + triple zp=point(p,i+1); + triple v1=c0-zi; + triple v2=c1-c0; + triple v3=zp-c1; + if(Cross(v0,v1) || Cross(v1,v2) || Cross(v2,v3)) return O; + v0=v3; + zi=zp; + } + return theta >= 0 ? normal : -normal; +} + +// Return a unit normal vector to a polygon with vertices in p. +triple normal(triple[] p) +{ + triple normal; + real fuzz=sqrtEpsilon*abs(maxbound(p)-minbound(p)); + real absnormal; + real theta; + + bool Cross(triple a, triple b) { + if(abs(a) >= fuzz && abs(b) >= fuzz) { + triple n=cross(unit(a),unit(b)); + real absn=abs(n); + n=unit(n); + if(absnormal > 0 && absn > sqrtEpsilon && + abs(normal-n) > sqrtEpsilon && abs(normal+n) > sqrtEpsilon) + return true; + else { + int sign=dot(n,normal) >= 0 ? 1 : -1; + theta += sign*asin1(absn); + if(absn > absnormal) { + absnormal=absn; + normal=n; + theta=sign*theta; + } + } + } + return false; + } + + if(p.length <= 0) return O; + + triple zi=p[0]; + triple v0=zi-p[p.length-1]; + for(int i=0; i < p.length-1; ++i) { + triple zp=p[i+1]; + triple v1=zp-zi; + if(Cross(v0,v1)) return O; + v0=v1; + zi=zp; + } + return theta >= 0 ? normal : -normal; +} + +// Transforms that map XY plane to YX, YZ, ZY, ZX, and XZ planes. +restricted transform3 XY=identity4; +restricted transform3 YX=rotate(-90,O,Z); +restricted transform3 YZ=rotate(90,O,Z)*rotate(90,O,X); +restricted transform3 ZY=rotate(-90,O,X)*YZ; +restricted transform3 ZX=rotate(-90,O,Z)*rotate(-90,O,Y); +restricted transform3 XZ=rotate(-90,O,Y)*ZX; + +private transform3 flip(transform3 t, triple X, triple Y, triple Z, + projection P) +{ + static transform3 flip(triple v) { + static real s(real x) {return x > 0 ? -1 : 1;} + return scale(s(v.x),s(v.y),s(v.z)); + } + + triple u=unit(P.normal); + triple up=unit(perp(P.up,u)); + bool upright=dot(Z,u) >= 0; + if(dot(Y,up) < 0) { + t=flip(Y)*t; + upright=!upright; + } + return upright ? t : flip(X)*t; +} + +restricted transform3 XY(projection P=currentprojection) +{ + return flip(XY,X,Y,Z,P); +} + +restricted transform3 YX(projection P=currentprojection) +{ + return flip(YX,Y,X,Z,P); +} + +restricted transform3 YZ(projection P=currentprojection) +{ + return flip(YZ,Y,Z,X,P); +} + +restricted transform3 ZY(projection P=currentprojection) +{ + return flip(ZY,Z,Y,X,P); +} + +restricted transform3 ZX(projection P=currentprojection) +{ + return flip(ZX,Z,X,Y,P); +} + +restricted transform3 XZ(projection P=currentprojection) +{ + return flip(XZ,X,Z,Y,P); +} + +// Transform3 that projects in direction dir onto plane with normal n +// through point O. +transform3 planeproject(triple n, triple O=O, triple dir=n) +{ + real a=n.x, b=n.y, c=n.z; + real u=dir.x, v=dir.y, w=dir.z; + real delta=1.0/(a*u+b*v+c*w); + real d=-(a*O.x+b*O.y+c*O.z)*delta; + return new real[][] { + {(b*v+c*w)*delta,-b*u*delta,-c*u*delta,-d*u}, + {-a*v*delta,(a*u+c*w)*delta,-c*v*delta,-d*v}, + {-a*w*delta,-b*w*delta,(a*u+b*v)*delta,-d*w}, + {0,0,0,1} + }; +} + +// Transform3 that projects in direction dir onto plane defined by p. +transform3 planeproject(path3 p, triple dir=O) +{ + triple n=normal(p); + return planeproject(n,point(p,0),dir == O ? n : dir); +} + +// Transform for projecting onto plane through point O with normal cross(u,v). +transform transform(triple u, triple v, triple O=O, + projection P=currentprojection) +{ + transform3 t=P.t; + static real[] O={0,0,0,1}; + real[] tO=t*O; + real tO3=tO[3]; + real factor=1/tO3^2; + real[] x=(tO3*t[0]-tO[0]*t[3])*factor; + real[] y=(tO3*t[1]-tO[1]*t[3])*factor; + triple x=(x[0],x[1],x[2]); + triple y=(y[0],y[1],y[2]); + u=unit(u); + v=unit(v); + return (0,0,dot(u,x),dot(v,x),dot(u,y),dot(v,y)); +} + +// Project Label onto plane through point O with normal cross(u,v). +Label project(Label L, triple u, triple v, triple O=O, + projection P=currentprojection) { + Label L=L.copy(); + L.position=project(O,P.t); + L.transform(transform(u,v,O,P)); + return L; +} + +path3 operator cast(guide3 g) {return solve(g);} +path3 operator cast(triple v) {return path3(v);} + +guide3[] operator cast(triple[] v) +{ + return sequence(new guide3(int i) {return v[i];},v.length); +} + +path3[] operator cast(triple[] v) +{ + return sequence(new path3(int i) {return v[i];},v.length); +} + +path3[] operator cast(guide3[] g) +{ + return sequence(new path3(int i) {return solve(g[i]);},g.length); +} + +guide3[] operator cast(path3[] g) +{ + return sequence(new guide3(int i) {return g[i];},g.length); +} + +void write(file file, string s="", explicit guide3[] x, suffix suffix=none) +{ + write(file,s,(path3[]) x,suffix); +} + +void write(string s="", explicit guide3[] x, suffix suffix=endl) +{ + write(stdout,s,(path3[]) x,suffix); +} + +triple point(explicit guide3 g, int t) { + flatguide3 f; + g(f); + int n=f.size(); + return f.nodes[adjustedIndex(t,n,f.cyclic())]; +} + +triple[] dirSpecifier(guide3 g, int t) +{ + flatguide3 f; + g(f); + int n=f.size(); + checkEmpty(n); + if(f.cyclic()) t=t % n; + else if(t < 0 || t >= n-1) return new triple[]; + return new triple[] {f.out[t].dir,f.in[t].dir}; +} + +triple[] controlSpecifier(guide3 g, int t) { + flatguide3 f; + g(f); + int n=f.size(); + checkEmpty(n); + if(f.cyclic()) t=t % n; + else if(t < 0 || t >= n-1) return new triple[]; + control c=f.control[t]; + if(c.active) return new triple[] {c.post,c.pre}; + else return new triple[]; +} + +tensionSpecifier tensionSpecifier(guide3 g, int t) +{ + flatguide3 f; + g(f); + int n=f.size(); + checkEmpty(n); + if(f.cyclic()) t=t % n; + else if(t < 0 || t >= n-1) return operator tension(1,1,false); + Tension T=f.Tension[t]; + return operator tension(T.out,T.in,T.atLeast); +} + +real[] curlSpecifier(guide3 g, int t) +{ + flatguide3 f; + g(f); + int n=f.size(); + checkEmpty(n); + if(f.cyclic()) t=t % n; + else if(t < 0 || t >= n-1) return new real[]; + return new real[] {f.out[t].gamma,f.in[t].gamma}; +} + +guide3 reverse(guide3 g) +{ + flatguide3 f; + bool cyclic=cyclic(g); + g(f); + + if(f.precyclic()) + return reverse(solve(g)); + + int n=f.size(); + checkEmpty(n); + guide3 G; + if(n >= 0) { + int start=cyclic ? n : n-1; + for(int i=start; i > 0; --i) { + G=G..f.nodes[i]; + control c=f.control[i-1]; + if(c.active) + G=G..operator controls(c.pre,c.post); + else { + dir in=f.in[i-1]; + triple d=in.dir; + if(d != O) G=G..operator spec(-d,JOIN_OUT); + else if(in.Curl) G=G..operator curl(in.gamma,JOIN_OUT); + dir out=f.out[i-1]; + triple d=out.dir; + if(d != O) G=G..operator spec(-d,JOIN_IN); + else if(out.Curl) G=G..operator curl(out.gamma,JOIN_IN); + } + } + if(cyclic) G=G..cycle; + else G=G..f.nodes[0]; + } + return G; +} + +triple intersectionpoint(path3 p, path3 q, real fuzz=-1) +{ + real[] t=intersect(p,q,fuzz); + if(t.length == 0) abort("paths do not intersect"); + return point(p,t[0]); +} + +// return an array containing all intersection points of p and q +triple[] intersectionpoints(path3 p, path3 q, real fuzz=-1) +{ + real[][] t=intersections(p,q,fuzz); + return sequence(new triple(int i) {return point(p,t[i][0]);},t.length); +} + +triple[] intersectionpoints(explicit path3[] p, explicit path3[] q, + real fuzz=-1) +{ + triple[] v; + for(int i=0; i < p.length; ++i) + for(int j=0; j < q.length; ++j) + v.append(intersectionpoints(p[i],q[j],fuzz)); + return v; +} + +path3 operator &(path3 p, cycleToken tok) +{ + int n=length(p); + if(n < 0) return nullpath3; + triple a=point(p,0); + triple b=point(p,n); + return subpath(p,0,n-1)..controls postcontrol(p,n-1) and precontrol(p,n).. + cycle; +} + +// return the point on path3 p at arclength L +triple arcpoint(path3 p, real L) +{ + return point(p,arctime(p,L)); +} + +// return the point on path3 p at arclength L +triple arcpoint(path3 p, real L) +{ + return point(p,arctime(p,L)); +} + +// return the direction on path3 p at arclength L +triple arcdir(path3 p, real L) +{ + return dir(p,arctime(p,L)); +} + +// return the time on path3 p at the relative fraction l of its arclength +real reltime(path3 p, real l) +{ + return arctime(p,l*arclength(p)); +} + +// return the point on path3 p at the relative fraction l of its arclength +triple relpoint(path3 p, real l) +{ + return point(p,reltime(p,l)); +} + +// return the direction of path3 p at the relative fraction l of its arclength +triple reldir(path3 p, real l) +{ + return dir(p,reltime(p,l)); +} + +// return the point on path3 p at half of its arclength +triple midpoint(path3 p) +{ + return relpoint(p,0.5); +} + +real relative(Label L, path3 g) +{ + return L.position.relative ? reltime(g,L.relative()) : L.relative(); +} + +// return the linear transformation that maps X,Y,Z to u,v,w. +transform3 transform3(triple u, triple v, triple w=cross(u,v)) +{ + return new real[][] { + {u.x,v.x,w.x,0}, + {u.y,v.y,w.y,0}, + {u.z,v.z,w.z,0}, + {0,0,0,1} + }; +} + +// return the rotation that maps Z to a unit vector u about cross(u,Z), +transform3 align(triple u) +{ + real a=u.x; + real b=u.y; + real c=u.z; + real d=a^2+b^2; + + if(d != 0) { + d=sqrt(d); + real e=1/d; + return new real[][] { + {-b*e,-a*c*e,a,0}, + {a*e,-b*c*e,b,0}, + {0,d,c,0}, + {0,0,0,1}}; + } + return c >= 0 ? identity(4) : diagonal(1,-1,-1,1); +} + +// return a rotation that maps X,Y to the projection plane. +transform3 transform3(projection P=currentprojection) +{ + 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); + return u != O ? transform3(u,v,w) : identity(4); +} + +triple[] triples(real[] x, real[] y, real[] z) +{ + if(x.length != y.length || x.length != z.length) + abort("arrays have different lengths"); + return sequence(new triple(int i) {return (x[i],y[i],z[i]);},x.length); +} + +path3[] operator cast(path3 p) +{ + return new path3[] {p}; +} + +path3[] operator cast(guide3 g) +{ + return new path3[] {(path3) g}; +} + +path3[] operator ^^ (path3 p, path3 q) +{ + return new path3[] {p,q}; +} + +path3[] operator ^^ (path3 p, explicit path3[] q) +{ + return concat(new path3[] {p},q); +} + +path3[] operator ^^ (explicit path3[] p, path3 q) +{ + return concat(p,new path3[] {q}); +} + +path3[] operator ^^ (explicit path3[] p, explicit path3[] q) +{ + return concat(p,q); +} + +path3[] operator * (transform3 t, explicit path3[] p) +{ + return sequence(new path3(int i) {return t*p[i];},p.length); +} + +triple[] operator * (transform3 t, triple[] v) +{ + return sequence(new triple(int i) {return t*v[i];},v.length); +} + +triple[][] operator * (transform3 t, triple[][] v) +{ + triple[][] V=new triple[v.length][]; + for(int i=0; i < v.length; ++i) { + triple[] vi=v[i]; + V[i]=sequence(new triple(int j) {return t*vi[j];},vi.length); + } + return V; +} + +triple min(explicit path3[] p) +{ + checkEmpty(p.length); + triple minp=min(p[0]); + for(int i=1; i < p.length; ++i) + minp=minbound(minp,min(p[i])); + return minp; +} + +triple max(explicit path3[] p) +{ + checkEmpty(p.length); + triple maxp=max(p[0]); + for(int i=1; i < p.length; ++i) + maxp=maxbound(maxp,max(p[i])); + return maxp; +} + +path3 randompath3(int n, bool cumulate=true, interpolate3 join=operator ..) +{ + guide3 g; + triple w; + for(int i=0; i <= n; ++i) { + triple z=(unitrand()-0.5,unitrand()-0.5,unitrand()-0.5); + if(cumulate) w += z; + else w=z; + g=join(g,w); + } + return g; +} + +path3[] box(triple v1, triple v2) +{ + return + (v1.x,v1.y,v1.z)-- + (v1.x,v1.y,v2.z)-- + (v1.x,v2.y,v2.z)-- + (v1.x,v2.y,v1.z)-- + (v1.x,v1.y,v1.z)-- + (v2.x,v1.y,v1.z)-- + (v2.x,v1.y,v2.z)-- + (v2.x,v2.y,v2.z)-- + (v2.x,v2.y,v1.z)-- + (v2.x,v1.y,v1.z)^^ + (v2.x,v2.y,v1.z)-- + (v1.x,v2.y,v1.z)^^ + (v1.x,v2.y,v2.z)-- + (v2.x,v2.y,v2.z)^^ + (v2.x,v1.y,v2.z)-- + (v1.x,v1.y,v2.z); +} + +restricted path3[] unitbox=box(O,(1,1,1)); +restricted path3 unitcircle3=X..Y..-X..-Y..cycle; +restricted path3 unitsquare3=O--X--X+Y--Y--cycle; + +path3 circle(triple c, real r, triple normal=Z) +{ + path3 p=normal == Z ? unitcircle3 : align(unit(normal))*unitcircle3; + return shift(c)*scale3(r)*p; +} + +// return an arc centered at c from triple v1 to v2 (assuming |v2-c|=|v1-c|), +// drawing in the given direction. +// The normal must be explicitly specified if c and the endpoints are colinear. +path3 arc(triple c, triple v1, triple v2, triple normal=O, bool direction=CCW) +{ + v1 -= c; + real r=abs(v1); + v1=unit(v1); + v2=unit(v2-c); + + if(normal == O) { + normal=cross(v1,v2); + if(normal == O) abort("explicit normal required for these endpoints"); + } + + transform3 T; + bool align=normal != Z; + if(align) { + T=align(unit(normal)); + transform3 Tinv=transpose(T); + v1=Tinv*v1; + v2=Tinv*v2; + } + + string invalidnormal="invalid normal vector"; + real fuzz=sqrtEpsilon*max(abs(v1),abs(v2)); + if(abs(v1.z) > fuzz || abs(v2.z) > fuzz) + abort(invalidnormal); + + real[] t1=intersect(unitcircle3,O--2*(v1.x,v1.y,0)); + real[] t2=intersect(unitcircle3,O--2*(v2.x,v2.y,0)); + + if(t1.length == 0 || t2.length == 0) + abort(invalidnormal); + + real t1=t1[0]; + real t2=t2[0]; + int n=length(unitcircle3); + if(direction) { + if (t1 >= t2) t1 -= n; + } else if(t2 >= t1) t2 -= n; + + path3 p=subpath(unitcircle3,t1,t2); + if(align) p=T*p; + return shift(c)*scale3(r)*p; +} + +// return an arc centered at c with radius r from c+r*dir(theta1,phi1) to +// c+r*dir(theta2,phi2) in degrees, drawing in the given direction +// relative to the normal vector cross(dir(theta1,phi1),dir(theta2,phi2)). +// The normal must be explicitly specified if c and the endpoints are colinear. +path3 arc(triple c, real r, real theta1, real phi1, real theta2, real phi2, + triple normal=O, bool direction) +{ + return arc(c,c+r*dir(theta1,phi1),c+r*dir(theta2,phi2),normal,direction); +} + +// return an arc centered at c with radius r from c+r*dir(theta1,phi1) to +// c+r*dir(theta2,phi2) in degrees, drawing drawing counterclockwise +// relative to the normal vector cross(dir(theta1,phi1),dir(theta2,phi2)) +// iff theta2 > theta1 or (theta2 == theta1 and phi2 >= phi1). +// The normal must be explicitly specified if c and the endpoints are colinear. +path3 arc(triple c, real r, real theta1, real phi1, real theta2, real phi2, + triple normal=O) +{ + return arc(c,r,theta1,phi1,theta2,phi2,normal, + theta2 > theta1 || (theta2 == theta1 && phi2 >= phi1) ? CCW : CW); +} + +private real epsilon=1000*realEpsilon; + +// Return a representation of the plane through point O with normal cross(u,v). +path3 plane(triple u, triple v, triple O=O) +{ + return O--O+u--O+u+v--O+v--cycle; +} + +// PRC/OpenGL support + +include three_light; + +void draw(frame f, path3 g, material p=currentpen, light light=nolight, + string name="", render render=defaultrender, + projection P=currentprojection); + +void begingroup3(frame f, string name="", render render=defaultrender, + triple center=O, int interaction=0) +{ + _begingroup3(f,name,render.compression,render.granularity,render.closed, + render.tessellate,render.merge == false, + render.merge == true,center,interaction); +} + +void begingroup3(picture pic=currentpicture, string name="", + render render=defaultrender, + triple center=O, int interaction=0) +{ + pic.add(new void(frame f, transform3, picture pic, projection) { + if(is3D()) + begingroup3(f,name,render,center,interaction); + if(pic != null) + begingroup(pic); + },true); +} + +void endgroup3(picture pic=currentpicture) +{ + pic.add(new void(frame f, transform3, picture pic, projection) { + if(is3D()) + endgroup3(f); + if(pic != null) + endgroup(pic); + },true); +} + +void addPath(picture pic, path3 g, pen p) +{ + if(size(g) > 0) + pic.addBox(min(g),max(g),min3(p),max3(p)); +} + +include three_surface; +include three_margins; + +pair min(path3 p, projection P) +{ + path3 q=P.T.modelview*p; + if(P.infinity) + return xypart(min(q)); + return maxratio(q)/P.T.projection[3][2]; +} + +pair max(path3 p, projection P) +{ + path3 q=P.T.modelview*p; + if(P.infinity) + return xypart(max(q)); + return minratio(q)/P.T.projection[3][2]; +} + +pair min(frame f, projection P) +{ + frame g=P.T.modelview*f; + if(P.infinity) + return xypart(min3(g)); + return maxratio(g)/P.T.projection[3][2]; +} + +pair max(frame f, projection P) +{ + frame g=P.T.modelview*f; + if(P.infinity) + return xypart(max3(g)); + return minratio(g)/P.T.projection[3][2]; +} + +void draw(picture pic=currentpicture, Label L="", path3 g, + align align=NoAlign, material p=currentpen, margin3 margin=NoMargin3, + light light=nolight, string name="", render render=defaultrender) +{ + pen q=(pen) p; + pic.add(new void(frame f, transform3 t, picture pic, projection P) { + path3 G=margin(t*g,q).g; + if(is3D()) { + draw(f,G,p,light,name,render,null); + if(pic != null && size(G) > 0) + pic.addBox(min(G,P),max(G,P),min(q),max(q)); + } + if(pic != null) + draw(pic,project(G,P),q); + },true); + Label L=L.copy(); + L.align(align); + if(L.s != "") { + L.p(q); + label(pic,L,g); + } + addPath(pic,g,q); +} + +include three_tube; + +draw=new void(frame f, path3 g, material p=currentpen, + light light=nolight, string name="", + render render=defaultrender, + projection P=currentprojection) { + pen q=(pen) p; + if(is3D()) { + p=material(p); + real width=linewidth(q); + void drawthick(path3 g) { + if(settings.thick) { + if(width > 0) { + bool prc=prc(); + void cylinder(transform3) {}; + void sphere(transform3, bool half) {}; + void disk(transform3) {}; + void pipe(path3, path3); + if(prc) { + cylinder=new void(transform3 t) {drawPRCcylinder(f,t,p,light);}; + sphere=new void(transform3 t, bool half) + {drawPRCsphere(f,t,half,p,light,render);}; + disk=new void(transform3 t) {draw(f,t*unitdisk,p,light,render);}; + pipe=new void(path3 center, path3 g) + {drawPRCtube(f,center,g,p,light);}; + } + real linecap=linecap(q); + real r=0.5*width; + bool open=!cyclic(g); + int L=length(g); + triple g0=point(g,0); + triple gL=point(g,L); + if(open && L > 0) { + if(linecap == 2) { + g0 -= r*dir(g,0); + gL += r*dir(g,L); + g=g0..g..gL; + L += 2; + } + } + tube T=tube(g,width,render,cylinder,sphere,pipe); + path3 c=T.center; + if(L >= 0) { + if(open) { + int Lc=length(c); + triple c0=point(c,0); + triple cL=point(c,Lc); + triple dir0=dir(g,0); + triple dirL=dir(g,L); + triple dirc0=dir(c,0); + triple dircL=dir(c,Lc); + transform3 t0=shift(g0)*align(-dir0); + transform3 tL=shift(gL)*align(dirL); + transform3 tc0=shift(c0)*align(-dirc0); + transform3 tcL=shift(cL)*align(dircL); + if(linecap == 0 || linecap == 2){ + transform3 scale2r=scale(r,r,1); + T.s.append(t0*scale2r*unitdisk); + disk(tc0*scale2r); + if(L > 0) { + T.s.append(tL*scale2r*unitdisk); + disk(tcL*scale2r); + } + } else if(linecap == 1) { + transform3 scale3r=scale3(r); + T.s.append(t0*scale3r* + (dir0 != O ? unithemisphere : unitsphere)); + sphere(tc0*scale3r,half=straight(c,0)); + if(L > 0) { + T.s.append(tL*scale3r* + (dirL != O ? unithemisphere : unitsphere)); + sphere(tcL*scale3r,half=straight(c,Lc-1)); + } + } + } + if(opacity(q) == 1) + _draw(f,c,q); + } + for(patch s : T.s.s) + draw3D(f,s,p,light,prc=false); + } else _draw(f,g,q); + } else _draw(f,g,q); + } + if(q != nullpen) + begingroup3(f,name == "" ? "curve" : name,render); + if(linetype(q).length == 0) drawthick(g); + else { + real[] dash=linetype(adjust(q,arclength(g),cyclic(g))); + if(sum(dash) > 0) { + dash.cyclic=true; + real offset=offset(q); + real L=arclength(g); + int i=0; + real l=offset; + while(l <= L) { + real t1=arctime(g,l); + l += dash[i]; + real t2=arctime(g,min(l,L)); + drawthick(subpath(g,t1,t2)); + ++i; + l += dash[i]; + ++i; + } + } + } + if(q != nullpen) + endgroup3(f); + } else draw(f,project(g,P),q); +}; + +void draw(frame f, explicit path3[] g, material p=currentpen, + light light=nolight, string name="", + render render=defaultrender, projection P=currentprojection) +{ + if(g.length > 1) + begingroup3(f,name == "" ? "curve" : name,render); + for(int i=0; i < g.length; ++i) + draw(f,g[i],p,light,partname(i),render,P); + if(g.length > 1) + endgroup3(f); +} + +void draw(picture pic=currentpicture, explicit path3[] g, + material p=currentpen, margin3 margin=NoMargin3, light light=nolight, + string name="", render render=defaultrender) +{ + if(g.length > 1) + begingroup3(pic,name == "" ? "curves" : name,render); + for(int i=0; i < g.length; ++i) + draw(pic,g[i],p,margin,light,partname(i),render); + if(g.length > 1) + endgroup3(pic); +} + +include three_arrows; + +void draw(picture pic=currentpicture, Label L="", path3 g, + align align=NoAlign, material p=currentpen, arrowbar3 arrow, + arrowbar3 bar=None, margin3 margin=NoMargin3, light light=nolight, + light arrowheadlight=currentlight, string name="", + render render=defaultrender) +{ + bool group=arrow != None || bar != None; + if(group) + begingroup3(pic,name,render); + bool drawpath=arrow(pic,g,p,margin,light,arrowheadlight); + if(bar(pic,g,p,margin,light,arrowheadlight) && drawpath) + draw(pic,L,g,align,p,margin,light,render); + if(group) + endgroup3(pic); + if(L.s != "") + label(pic,L,g,align,(pen) p); +} + +void draw(frame f, path3 g, material p=currentpen, arrowbar3 arrow, + light light=nolight, light arrowheadlight=currentlight, + string name="", render render=defaultrender, + projection P=currentprojection) +{ + picture pic; + bool group=arrow != None; + if(group) + begingroup3(f,name,render); + if(arrow(pic,g,p,NoMargin3,light,arrowheadlight)) + draw(f,g,p,light,render,P); + add(f,pic.fit()); + if(group) + endgroup3(f); +} + +void add(picture pic=currentpicture, void d(picture,transform3), + bool exact=false) +{ + pic.add(d,exact); +} + +// Fit the picture src using the identity transformation (so user +// coordinates and truesize coordinates agree) and add it about the point +// position to picture dest. +void add(picture dest, picture src, triple position, bool group=true, + bool above=true) +{ + dest.add(new void(picture f, transform3 t) { + f.add(shift(t*position)*src,group,above); + }); +} + +void add(picture src, triple position, bool group=true, bool above=true) +{ + add(currentpicture,src,position,group,above); +} + +// Align an arrow pointing to b from the direction dir. The arrow is +// 'length' PostScript units long. +void arrow(picture pic=currentpicture, Label L="", triple b, triple dir, + real length=arrowlength, align align=NoAlign, + pen p=currentpen, arrowbar3 arrow=Arrow3, margin3 margin=EndMargin3, + light light=nolight, light arrowheadlight=currentlight, + string name="", render render=defaultrender) +{ + Label L=L.copy(); + if(L.defaultposition) L.position(0); + L.align(L.align,dir); + L.p(p); + picture opic; + marginT3 margin=margin(b--b,p); // Extract margin.begin and margin.end + triple a=(margin.begin+length+margin.end)*unit(dir); + draw(opic,L,a--O,align,p,arrow,margin,light,arrowheadlight,name,render); + add(pic,opic,b); +} + +void arrow(picture pic=currentpicture, Label L="", triple b, pair dir, + real length=arrowlength, align align=NoAlign, + pen p=currentpen, arrowbar3 arrow=Arrow3, margin3 margin=EndMargin3, + light light=nolight, light arrowheadlight=currentlight, + string name="", render render=defaultrender, + projection P=currentprojection) +{ + arrow(pic,L,b,invert(dir,b,P),length,align,p,arrow,margin,light, + arrowheadlight,name,render); +} + +triple min3(picture pic, projection P=currentprojection) +{ + return pic.min3(P); +} + +triple max3(picture pic, projection P=currentprojection) +{ + return pic.max3(P); +} + +triple size3(picture pic, bool user=false, projection P=currentprojection) +{ + transform3 t=pic.calculateTransform3(P); + triple M=pic.max(t); + triple m=pic.min(t); + if(!user) return M-m; + t=inverse(t); + return t*M-t*m; +} + +triple point(frame f, triple dir) +{ + triple m=min3(f); + triple M=max3(f); + return m+realmult(rectify(dir),M-m); +} + +triple point(picture pic=currentpicture, triple dir, bool user=true, + projection P=currentprojection) +{ + triple min = pic.userMin(), max = pic.userMax(); + triple v=min+realmult(rectify(dir),max-min); + return user ? v : pic.calculateTransform3(P)*v; +} + +triple truepoint(picture pic=currentpicture, triple dir, bool user=true, + projection P=currentprojection) +{ + transform3 t=pic.calculateTransform3(P); + triple m=pic.min(t); + triple M=pic.max(t); + triple v=m+realmult(rectify(dir),M-m); + return user ? inverse(t)*v : v; +} + +void add(picture dest=currentpicture, object src, pair position=0, pair align=0, + bool group=true, filltype filltype=NoFill, bool above=true) +{ + if(prc()) + label(dest,src,position,align); + else if(settings.render == 0) + plain.add(dest,src,position,align,group,filltype,above); +} + +string cameralink(string label, string text="View Parameters") +{ + if(!prc() || Link == null) return ""; + return Link(label,text,"3Dgetview"); +} + +private struct viewpoint { + triple target,camera,up; + real angle; + void operator init(string s) { + s=replace(s,'\n'," "); + string[] S=split(s); + int pos(string s, string key) { + int pos=find(s,key); + if(pos < 0) return -1; + pos += length(key); + while(substr(s,pos,1) == " ") ++pos; + if(substr(s,pos,1) == "=") + return pos+1; + return -1; + } + triple C2C=X; + real ROO=1; + real ROLL=0; + angle=30; + int pos; + for(int k=0; k < S.length; ++k) { + if((pos=pos(S[k],"COO")) >= 0) + target=((real) substr(S[k],pos),(real) S[++k],(real) S[++k]); + else if((pos=pos(S[k],"C2C")) >= 0) + C2C=((real) substr(S[k],pos),(real) S[++k],(real) S[++k]); + else if((pos=pos(S[k],"ROO")) >= 0) + ROO=(real) substr(S[k],pos); + else if((pos=pos(S[k],"ROLL")) >= 0) + ROLL=(real) substr(S[k],pos); + else if((pos=pos(S[k],"AAC")) >= 0) + angle=(real) substr(S[k],pos); + } + camera=target+ROO*C2C; + triple u=unit(target-camera); + triple w=unit(Z-u.z*u); + up=rotate(ROLL,O,u)*w; + } +} + +projection perspective(string s) +{ + viewpoint v=viewpoint(s); + projection P=perspective(v.camera,v.up,v.target); + P.angle=v.angle; + P.absolute=true; + return P; +} + +private string Format(real x) +{ + // Work around movie15.sty division by zero bug; + // e.g. u=unit((1e-10,1e-10,0.9)); + if(abs(x) < 1e-9) x=0; + assert(abs(x) < 1e17,"Number too large: "+string(x)); + return format("%.18f",x,"C"); +} + +private string Format(triple v, string sep=" ") +{ + return Format(v.x)+sep+Format(v.y)+sep+Format(v.z); +} + +private string Format(real[] c) +{ + return Format((c[0],c[1],c[2])); +} + +private string format(triple v, string sep=" ") +{ + return string(v.x)+sep+string(v.y)+sep+string(v.z); +} + +private string projection(bool infinity, real viewplanesize) +{ + return "activeCamera=scene.cameras.getByIndex(0); +function asyProjection() {"+ + (infinity ? "activeCamera.projectionType=activeCamera.TYPE_ORTHOGRAPHIC;" : + "activeCamera.projectionType=activeCamera.TYPE_PERSPECTIVE;")+" +activeCamera.viewPlaneSize="+string(viewplanesize)+"; +activeCamera.binding=activeCamera.BINDING_"+(infinity ? "MAX" : "VERTICAL")+"; +} + +asyProjection(); + +handler=new CameraEventHandler(); +runtime.addEventHandler(handler); +handler.onEvent=function(event) +{ + asyProjection(); + scene.update(); +}"; +} + +private string projection(bool infinity, real viewplanesize) +{ + return "activeCamera=scene.cameras.getByIndex(0); +function asyProjection() {"+ + (infinity ? "activeCamera.projectionType=activeCamera.TYPE_ORTHOGRAPHIC;" : + "activeCamera.projectionType=activeCamera.TYPE_PERSPECTIVE;")+" +activeCamera.viewPlaneSize="+string(viewplanesize)+"; +activeCamera.binding=activeCamera.BINDING_"+(infinity ? "MAX" : "VERTICAL")+"; +} + +asyProjection(); + +handler=new CameraEventHandler(); +runtime.addEventHandler(handler); +handler.onEvent=function(event) +{ + asyProjection(); + scene.update(); +}"; +} + +private string billboard(int[] index, triple[] center) +{ + if(index.length == 0) return ""; + string s=" +var zero=new Vector3(0,0,0); +var nodes=scene.nodes; +var count=nodes.count; + +var index=new Array(); +for(i=0; i < count; i++) { + var node=nodes.getByIndex(i); + var name=node.name; + end=name.lastIndexOf(\".\")-1; + if(end > 0) { + if(name.substr(end,1) == \"\001\") { + start=name.lastIndexOf(\"-\")+1; + n=end-start; + if(n > 0) { + index[name.substr(start,n)]=i; + node.name=name.substr(0,start-1); + } + } + } +} + +var center=new Array( +"; + for(int i=0; i < center.length; ++i) + s += "Vector3("+format(center[i],",")+"), +"; + s += "); + +billboardHandler=new RenderEventHandler(); +billboardHandler.onEvent=function(event) +{ + var camera=scene.cameras.getByIndex(0); + var position=camera.position; + var direction=position.subtract(camera.targetPosition); + var up=camera.up.subtract(position); + + function f(i,k) { + j=index[i]; + if(j >= 0) { + var node=nodes.getByIndex(j); + var name=node.name; + var R=Matrix4x4(); + R.setView(zero,direction,up); + var c=center[k]; + var T=node.transform; + T.setIdentity(); + T.translateInPlace(c.scale(-1)); + T.multiplyInPlace(R); + T.translateInPlace(c); + } + } +"; + for(int i=0; i < index.length; ++i) + s += "f("+string(i)+","+string(index[i])+"); +"; + s += " + runtime.refresh(); +} + +runtime.addEventHandler(billboardHandler); + +runtime.refresh(); +"; + return s; +} + +string lightscript(light light) { + string script="for(var i=scene.lights.count-1; i >= 0; i--) + scene.lights.removeByIndex(i);"+'\n\n'; + for(int i=0; i < light.position.length; ++i) { + string Li="L"+string(i); + real[] diffuse=light.diffuse[i]; + script += Li+"=scene.createLight();"+'\n'+ + Li+".direction.set("+format(-light.position[i],",")+");"+'\n'+ + Li+".color.set("+format((diffuse[0],diffuse[1],diffuse[2]),",")+");"+'\n'; + } + // Work around initialization bug in Adobe Reader 8.0: + return script +" +scene.lightScheme=scene.LIGHT_MODE_HEADLAMP; +scene.lightScheme=scene.LIGHT_MODE_FILE; +"; +} + +void writeJavaScript(string name, string preamble, string script) +{ + file out=output(name); + write(out,preamble); + if(script != "") { + write(out,endl); + file in=input(script); + while(true) { + string line=in; + if(eof(in)) break; + write(out,line,endl); + } + } + close(out); + if(settings.verbose > 1) write("Wrote "+name); + if(!settings.inlinetex) + file3.push(name); +} + +pair viewportmargin(pair lambda) +{ + return maxbound(0.5*(viewportsize-lambda),viewportmargin); +} + +string embed3D(string label="", string text=label, string prefix, + frame f, string format="", + real width=0, real height=0, + string options="", string script="", + light light=currentlight, projection P=currentprojection) +{ + if(!prc(format) || Embed == null) return ""; + + if(width == 0) width=settings.paperwidth; + if(height == 0) height=settings.paperheight; + + if(script == "") script=defaultembed3Dscript; + + // Adobe Reader doesn't appear to support user-specified viewport lights. + string lightscript=light.on() && light.viewport ? "" : lightscript(light); + + real viewplanesize; + if(P.infinity) { + triple lambda=max3(f)-min3(f); + pair margin=viewportmargin((lambda.x,lambda.y)); + viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/P.zoom; + } else + if(!P.absolute) P.angle=2*aTan(Tan(0.5*P.angle)); + + int[] index; + triple[] center; + shipout3(prefix,f,index,center); + + string name=prefix+".js"; + writeJavaScript(name,lightscript+projection(P.infinity,viewplanesize)+ + billboard(index,center),script); + + if(!settings.inlinetex) + file3.push(prefix+".prc"); + + triple target=P.target; + if(P.viewportshift != 0) { + triple lambda=max3(f)-min3(f); + target -= (P.viewportshift.x*lambda.x/P.zoom, + P.viewportshift.y*lambda.y/P.zoom,0); + } + + triple v=P.vector(); + triple u=unit(v); + triple w=Z-u.z*u; + real roll; + if(abs(w) > sqrtEpsilon) { + w=unit(w); + triple up=unit(perp(P.up,u)); + roll=degrees(acos1(dot(up,w)))*sgn(dot(cross(up,w),u)); + } else roll=0; + + string options3=light.viewport ? "3Dlights=Headlamp" : "3Dlights=File"; + if(defaultembed3Doptions != "") options3 += ","+defaultembed3Doptions; + if((settings.render < 0 || !settings.embed) && settings.auto3D) + options3 += ",poster"; + options3 += ",text={"+text+"},label="+label+ + ",toolbar="+(settings.toolbar ? "true" : "false")+ + ",3Daac="+Format(P.angle)+ + ",3Dc2c="+Format(u)+ + ",3Dcoo="+Format(target)+ + ",3Droll="+Format(roll)+ + ",3Droo="+Format(abs(v))+ + ",3Dbg="+Format(light.background()); + if(options != "") options3 += ","+options; + if(settings.inlinetex) + prefix=jobname(prefix); + options3 += ",3Djscript="+prefix+".js"; + + return Embed(prefix+".prc",options3,width,height); +} + +struct scene +{ + frame f; + transform3 t; + projection P; + bool adjusted; + real width,height; + pair viewportmargin; + transform3 T=identity4; + picture pic2; + + void operator init(frame f, real width, real height, + projection P=currentprojection) { + this.f=f; + this.t=identity4; + this.P=P; + this.width=width; + this.height=height; + } + + void operator init(picture pic, real xsize=pic.xsize, real ysize=pic.ysize, + bool keepAspect=pic.keepAspect, bool is3D=true, + projection P=currentprojection) { + real xsize3=pic.xsize3, ysize3=pic.ysize3, zsize3=pic.zsize3; + bool warn=true; + + if(xsize3 == 0 && ysize3 == 0 && zsize3 == 0) { + xsize3=ysize3=zsize3=max(xsize,ysize); + warn=false; + } + + if(P.absolute) + this.P=P.copy(); + else if(P.showtarget) + draw(pic,P.target,nullpen); + + t=pic.scaling(xsize3,ysize3,zsize3,keepAspect,warn); + adjusted=false; + triple m=pic.min(t); + triple M=pic.max(t); + + if(!P.absolute) { + this.P=t*P; + if(this.P.center && settings.render != 0) { + triple target=0.5*(m+M); + this.P.target=target; + this.P.calculate(); + } + if(this.P.autoadjust || this.P.infinity) + adjusted=adjusted | this.P.adjust(m,M); + } + + bool scale=xsize != 0 || ysize != 0; + bool scaleAdjust=scale && this.P.autoadjust; + bool noAdjust=(this.P.absolute || !scaleAdjust); + + if(pic.bounds3.exact && noAdjust) + this.P.bboxonly=false; + + f=pic.fit3(t,pic.bounds3.exact ? pic2 : null,this.P); + + if(!pic.bounds3.exact) { + if(noAdjust) + this.P.bboxonly=false; + + transform3 s=pic.scale3(f,xsize3,ysize3,zsize3,keepAspect); + t=s*t; + this.P=s*this.P; + f=pic.fit3(t,pic2,this.P); + } + + if(is3D || scale) { + pic2.bounds.exact=true; + transform s=pic2.scaling(xsize,ysize,keepAspect); + + pair m2=pic2.min(s); + pair M2=pic2.max(s); + pair lambda=M2-m2; + viewportmargin=viewportmargin(lambda); + width=ceil(lambda.x+2*viewportmargin.x); + height=ceil(lambda.y+2*viewportmargin.y); + + if(!this.P.absolute) { + if(scaleAdjust) { + pair v=(s.xx,s.yy); + transform3 T=this.P.t; + pair x=project(X,T); + pair y=project(Y,T); + pair z=project(Z,T); + real f(pair a, pair b) { + return b == 0 ? (0.5*(a.x+a.y)) : + (b.x^2*a.x+b.y^2*a.y)/(b.x^2+b.y^2); + } + transform3 s=keepAspect ? scale3(min(f(v,x),f(v,y),f(v,z))) : + xscale3(f(v,x))*yscale3(f(v,y))*zscale3(f(v,z)); + s=shift(this.P.target)*s*shift(-this.P.target); + t=s*t; + this.P=s*this.P; + this.P.bboxonly=false; + if(!is3D) pic2.erase(); + f=pic.fit3(t,is3D ? null : pic2,this.P); + } + + if(this.P.autoadjust || this.P.infinity) + adjusted=adjusted | this.P.adjust(min3(f),max3(f)); + } + } + } + + // Choose the angle to be just large enough to view the entire image. + real angle(projection P) { + T=identity4; + real h=-0.5*P.target.z; + pair r,R; + real diff=realMax; + pair s; + int i; + do { + r=minratio(f); + R=maxratio(f); + pair lasts=s; + if(P.autoadjust) { + s=r+R; + if(s != 0) { + transform3 t=shift(h*s.x,h*s.y,0); + f=t*f; + T=t*T; + adjusted=true; + } + } + diff=abs(s-lasts); + ++i; + } while (diff > angleprecision && i < maxangleiterations); + real aspect=width > 0 ? height/width : 1; + real rx=-r.x*aspect; + real Rx=R.x*aspect; + real ry=-r.y; + real Ry=R.y; + if(!P.autoadjust) { + if(rx > Rx) Rx=rx; + else rx=Rx; + if(ry > Ry) Ry=ry; + else ry=Ry; + } + return (1+angleprecision)*max(aTan(rx)+aTan(Rx),aTan(ry)+aTan(Ry)); + } +} + +object embed(string label="", string text=label, string prefix=defaultfilename, + scene S, string format="", bool view=true, string options="", + string script="", light light=currentlight) +{ + object F; + transform3 modelview; + projection P=S.P; + transform3 tinv=inverse(S.t); + + projection Q; + modelview=P.T.modelview; + if(P.absolute) { + Q=modelview*P; + } else { + triple target=P.target; + S.f=modelview*S.f; + P=modelview*P; + Q=P.copy(); + light=modelview*light; + + if(P.infinity) { + triple m=min3(S.f); + triple M=max3(S.f); + triple lambda=M-m; + S.viewportmargin=viewportmargin((lambda.x,lambda.y)); + S.width=ceil(lambda.x+2*S.viewportmargin.x); + S.height=ceil(lambda.y+2*S.viewportmargin.y); + S.f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*S.f; // Eye will be at (0,0,0) + } else { + if(P.angle == 0) { + P.angle=S.angle(P); + modelview=S.T*modelview; + if(S.viewportmargin.y != 0) + P.angle=2*aTan(Tan(0.5*P.angle)-S.viewportmargin.y/P.target.z); + } + Q.angle=P.angle; + if(settings.verbose > 0) { + transform3 inv=inverse(modelview); + if(S.adjusted) + write("adjusting camera to ",tinv*inv*P.camera); + target=inv*P.target; + } + P=S.T*P; + } + if(settings.verbose > 0) { + if((P.center && settings.render != 0) || (!P.infinity && P.autoadjust)) + write("adjusting target to ",tinv*target); + } + } + + if(prefix == "") prefix=outprefix(); + bool prc=prc(format); + bool preview=settings.render > 0; + if(prc) { + // The movie15.sty package cannot handle spaces or dots in filenames. + string dir=stripfile(prefix); + prefix=dir+replace(stripdirectory(prefix), + new string[][]{{" ","_"},{".","_"}}); + if(settings.embed || nativeformat() == "pdf") + prefix += "+"+(string) file3.length; + } else + preview=false; + if(preview || (!prc && settings.render != 0)) { + frame f=S.f; + triple m,M; + real zcenter; + real r; + if(P.absolute) { + f=modelview*f; + m=min3(f); + M=max3(f); + r=0.5*abs(M-m); + zcenter=0.5*(M.z+m.z); + } else { + m=min3(f); + M=max3(f); + zcenter=P.target.z; + r=P.distance(m,M); + } + M=(M.x,M.y,zcenter+r); + m=(m.x,m.y,zcenter-r); + + if(P.infinity) { + triple margin=(S.viewportmargin.x,S.viewportmargin.y,0); + M += margin; + m -= margin; + } else if(M.z >= 0) abort("camera too close"); + + shipout3(prefix,f,preview ? nativeformat() : format, + S.width-defaultrender.margin,S.height-defaultrender.margin, + P.infinity ? 0 : 2aTan(Tan(0.5*P.angle)*P.zoom), + P.zoom,m,M,P.viewportshift, + tinv*inverse(modelview)*shift(0,0,zcenter),light.background(), + P.absolute ? (modelview*light).position : light.position, + light.diffuse,light.ambient,light.specular, + light.viewport,view && !preview); + if(!preview) return F; + } + + string image; + if((preview || (prc && settings.render == 0)) && settings.embed) { + image=prefix; + if(settings.inlinetex) image += "_0"; + if(!preview && !shipped && !S.pic2.empty2()) { + transform T=S.pic2.scaling(S.width,S.height); + shipout(image,S.pic2.fit(T),newframe,nativeformat(),false,false,null); + } + + image += "."+nativeformat(); + if(!settings.inlinetex) file3.push(image); + image=graphic(image,"hiresbb"); + } + if(prc) { + if(!P.infinity && P.viewportshift != 0) + warning("offaxis", + "PRC does not support off-axis projections; use pan instead of +shift"); + F.L=embed3D(label,text=image,prefix,S.f,format, + S.width-2,S.height-2,options,script,light,Q); + } + return F; +} + +object embed(string label="", string text=label, + string prefix=defaultfilename, + picture pic, string format="", + real xsize=pic.xsize, real ysize=pic.ysize, + bool keepAspect=pic.keepAspect, bool view=true, string options="", + string script="", light light=currentlight, + projection P=currentprojection) +{ + bool is3D=is3D(format); + scene S=scene(pic,xsize,ysize,keepAspect,is3D,P); + if(is3D) + return embed(label,text,prefix,S,format,view,options,script,light); + else { + object F; + transform T=S.pic2.scaling(xsize,ysize,keepAspect); + F.f=pic.fit(scale(S.t[0][0])*T); + add(F.f,S.pic2.fit()); + return F; + } +} + +object embed(string label="", string text=label, + string prefix=defaultfilename, + frame f, string format="", real width=0, real height=0, + bool view=true, string options="", string script="", + light light=currentlight, projection P=currentprojection) +{ + if(is3D(format)) + return embed(label,text,prefix,scene(f,width,height,P),format,view,options, + script,light); + else { + object F; + F.f=f; + return F; + } +} + +embed3=new object(string prefix, frame f, string format, string options, + string script, light light, projection P) { + return embed(prefix=prefix,f,format,options,script,light,P); +}; + +frame embedder(object embedder(string prefix, string format), + string prefix, string format, bool view, light light) +{ + frame f; + bool prc=prc(format); + if(!prc && settings.render != 0 && !view) { + static int previewcount=0; + bool keep=prefix != ""; + prefix=outprefix(prefix)+"+"+(string) previewcount; + ++previewcount; + format=nativeformat(); + if(!keep) file3.push(prefix+"."+format); + } + object F=embedder(prefix,format); + if(prc) + label(f,F.L); + else { + if(settings.render == 0) { + add(f,F.f); + if(light.background != nullpen) + box(f,light.background,Fill,above=false); + } else if(!view) + label(f,graphic(prefix,"hiresbb")); + } + return f; +} + +currentpicture.fitter=new frame(string prefix, picture pic, string format, + real xsize, real ysize, bool keepAspect, + bool view, string options, string script, + light light, projection P) { + frame f; + bool empty3=pic.empty3(); + if(!empty3) f=embedder(new object(string prefix, string format) { + return embed(prefix=prefix,pic,format,xsize,ysize,keepAspect,view, + options,script,light,P); + },prefix,format,view,light); + if(is3D(format) || empty3) add(f,pic.fit2(xsize,ysize,keepAspect)); + return f; +}; + +frame embedder(string prefix, frame f, string format, real width, real height, + bool view, string options, string script, light light, + projection P) +{ + return embedder(new object(string prefix, string format) { + return embed(prefix=prefix,f,format,width,height,view,options,script, + light,P); + },prefix,format,view,light); +} + +projection[][] ThreeViewsUS={{TopView}, + {FrontView,RightView}}; + +projection[][] SixViewsUS={{null,TopView}, + {LeftView,FrontView,RightView,BackView}, + {null,BottomView}}; + +projection[][] ThreeViewsFR={{RightView,FrontView}, + {null,TopView}}; + +projection[][] SixViewsFR={{null,BottomView}, + {RightView,FrontView,LeftView,BackView}, + {null,TopView}}; + +projection[][] ThreeViews={{FrontView,TopView,RightView}}; + +projection[][] SixViews={{FrontView,TopView,RightView}, + {BackView,BottomView,LeftView}}; + +void addViews(picture dest, picture src, projection[][] views=SixViewsUS, + bool group=true, filltype filltype=NoFill) +{ + frame[][] F=array(views.length,new frame[]); + pair[][] M=array(views.length,new pair[]); + pair[][] m=array(views.length,new pair[]); + + for(int i=0; i < views.length; ++i) { + projection[] viewsi=views[i]; + frame[] Fi=F[i]; + pair[] Mi=M[i]; + pair[] mi=m[i]; + for(projection P : viewsi) { + if(P != null) { + frame f=src.fit(P); + mi.push(min(f)); + Mi.push(max(f)); + Fi.push(f); + } else { + pair Infinity=(infinity,infinity); + mi.push(Infinity); + Mi.push(-Infinity); + Fi.push(newframe); + } + } + } + + real[] my=new real[views.length]; + real[] My=new real[views.length]; + + int Nj=0; + for(int i=0; i < views.length; ++i) { + my[i]=minbound(m[i]).y; + My[i]=maxbound(M[i]).y; + Nj=max(Nj,views[i].length); + } + + real[] mx=array(Nj,infinity); + real[] Mx=array(Nj,-infinity); + for(int i=0; i < views.length; ++i) { + pair[] mi=m[i]; + pair[] Mi=M[i]; + for(int j=0; j < views[i].length; ++j) { + mx[j]=min(mx[j],mi[j].x); + Mx[j]=max(Mx[j],Mi[j].x); + } + } + + if(group) begingroup(dest); + + real y; + for(int i=0; i < views.length; ++i) { + real x; + pair[] mi=m[i]; + for(int j=0; j < views[i].length; ++j) { + if(size(F[i][j]) != 0) + add(dest,shift(x-mx[j],y+my[i])*F[i][j],filltype); + x += (Mx[j]-mx[j]); + } + y -= (My[i]-my[i]); + } + + if(group) endgroup(dest); +} + +void addViews(picture src, projection[][] views=SixViewsUS, bool group=true, + filltype filltype=NoFill) +{ + addViews(currentpicture,src,views,group,filltype); +} + +void addStereoViews(picture dest, picture src, bool group=true, + filltype filltype=NoFill, real eyetoview=defaulteyetoview, + bool leftright=true, projection P=currentprojection) +{ + triple v=P.vector(); + triple h=0.5*abs(v)*eyetoview*unit(cross(P.up,v)); + projection leftEye=P.copy(); + leftEye.camera -= h; + leftEye.calculate(); + projection rightEye=P.copy(); + rightEye.camera += h; + rightEye.calculate(); + addViews(dest,src,leftright ? + new projection[][] {{leftEye,rightEye}} : + new projection[][] {{rightEye,leftEye}},group,filltype); +} + +void addStereoViews(picture src, bool group=true, + filltype filltype=NoFill, + real eyetoview=defaulteyetoview, bool leftright=true, + projection P=currentprojection) +{ + addStereoViews(currentpicture,src,group,filltype,eyetoview,leftright,P); +} + +// Fit an array of 3D pictures simultaneously using the sizing of picture all. +frame[] fit3(string prefix="", picture[] pictures, picture all, + string format="", bool view=true, string options="", + string script="", light light=currentlight, + projection P=currentprojection) +{ + frame[] out; + scene S=scene(all,P); + triple m=all.min(S.t); + triple M=all.max(S.t); + out=new frame[pictures.length]; + int i=0; + bool reverse=settings.reverse; + settings.animating=true; + + for(picture pic : pictures) { + picture pic2; + frame f=pic.fit3(S.t,pic2,S.P); + if(settings.interrupt) break; + add(f,pic2.fit2()); + draw(f,m,nullpen); + draw(f,M,nullpen); + out[i]=f; + ++i; + } + + while(!settings.interrupt) { + for(int i=settings.reverse ? pictures.length-1 : 0; + i >= 0 && i < pictures.length && !settings.interrupt; + settings.reverse ? --i : ++i) { + frame f=embedder(prefix,out[i],format,S.width,S.height,view,options, + script,light,S.P); + if(!settings.loop) out[i]=f; + } + if(!settings.loop) break; + } + + settings.animating=false; + settings.interrupt=false; + settings.reverse=reverse; + + return out; +} + +// Fit an array of pictures simultaneously using the size of the first picture. +fit=new frame[](string prefix="", picture[] pictures, string format="", + bool view=true, string options="", string script="", + projection P=currentprojection) { + if(pictures.length == 0) + return new frame[]; + + picture all; + size(all,pictures[0]); + for(picture pic : pictures) + add(all,pic); + + return all.empty3() ? fit2(pictures,all) : + fit3(prefix,pictures,all,format,view,options,script,P); +}; + +// Add frame src to picture dest about position. +void add(picture dest=currentpicture, frame src, triple position) +{ + if(is3D(src)) { + dest.add(new void(frame f, transform3 t, picture, projection) { + add(f,shift(t*position)*src); + },true); + } else { + dest.add(new void(frame, transform3 t, picture pic, projection P) { + if(pic != null) { + pic.add(new void(frame f, transform T) { + add(f,T*shift(project(t*position,P))*src); + },true); + } + },true); + } + dest.addBox(position,position,min3(src),max3(src)); +} + +exitfcn currentexitfunction=atexit(); + +void exitfunction() +{ + if(currentexitfunction != null) currentexitfunction(); + if(!settings.keep) + for(int i=0; i < file3.length; ++i) + delete(file3[i]); + file3=new string[]; +} + +atexit(exitfunction); -- cgit v1.2.3