summaryrefslogtreecommitdiff
path: root/Master/texmf/asymptote/three_tube.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf/asymptote/three_tube.asy')
-rw-r--r--Master/texmf/asymptote/three_tube.asy223
1 files changed, 134 insertions, 89 deletions
diff --git a/Master/texmf/asymptote/three_tube.asy b/Master/texmf/asymptote/three_tube.asy
index c79ae3f619f..2eb5e4cf716 100644
--- a/Master/texmf/asymptote/three_tube.asy
+++ b/Master/texmf/asymptote/three_tube.asy
@@ -1,5 +1,6 @@
-void render(path3 s, real granularity=linegranularity, void f(path3, real))
+void render(path3 s, void f(path3, real), render render=defaultrender)
{
+ real granularity=render.tubegranularity;
void Split(triple z0, triple c0, triple c1, triple z1, real t0=0, real t1=1,
real depth=mantissaBits) {
if(depth > 0) {
@@ -65,8 +66,8 @@ rmf[] rmf(path3 g, real[] t)
return R;
}
-surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
- real[] x, real[] y, bool[][] cond={})
+private real[][][] bispline0(real[][] z, real[][] p, real[][] q, real[][] r,
+ real[] x, real[] y, bool[][] cond={})
{ // z[i][j] is the value at (x[i],y[j])
// p and q are the first derivatives with respect to x and y, respectively
// r is the second derivative ddu/dxdy
@@ -87,62 +88,49 @@ surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
}
}
- surface s=surface(count);
- s.index=new int[n][m];
- int k=-1;
+ real[][][] s=new real[count][][];
+ int k=0;
for(int i=0; i < n; ++i) {
+ int ip=i+1;
bool[] condi=all ? null : cond[i];
real xi=x[i];
real[] zi=z[i];
- real[] zp=z[i+1];
+ real[] zp=z[ip];
real[] ri=r[i];
- real[] rp=r[i+1];
+ real[] rp=r[ip];
real[] pi=p[i];
- real[] pp=p[i+1];
+ real[] pp=p[ip];
real[] qi=q[i];
- real[] qp=q[i+1];
- real xp=x[i+1];
+ real[] qp=q[ip];
+ real xp=x[ip];
real hx=(xp-xi)/3;
- int[] indexi=s.index[i];
for(int j=0; j < m; ++j) {
real yj=y[j];
- real yp=y[j+1];
+ int jp=j+1;
+ real yp=y[jp];
if(all || condi[j]) {
- triple[][] P=array(4,array(4,O));
real hy=(yp-yj)/3;
real hxy=hx*hy;
- // x and y directions
- for(int k=0; k < 4; ++k) {
- P[0][k] += xi*X;
- P[k][0] += yj*Y;
- P[1][k] += (xp+2*xi)/3*X;
- P[k][1] += (yp+2*yj)/3*Y;
- P[2][k] += (2*xp+xi)/3*X;
- P[k][2] += (2*yp+yj)/3*Y;
- P[3][k] += xp*X;
- P[k][3] += yp*Y;
- }
- // z: value
- P[0][0] += zi[j]*Z;
- P[3][0] += zp[j]*Z;
- P[0][3] += zi[j+1]*Z;
- P[3][3] += zp[j+1]*Z;
- // z: first derivative
- P[1][0] += (P[0][0].z+hx*pi[j])*Z;
- P[1][3] += (P[0][3].z+hx*pi[j+1])*Z;
- P[2][0] += (P[3][0].z-hx*pp[j])*Z;
- P[2][3] += (P[3][3].z-hx*pp[j+1])*Z;
- P[0][1] += (P[0][0].z+hy*qi[j])*Z;
- P[3][1] += (P[3][0].z+hy*qp[j])*Z;
- P[0][2] += (P[0][3].z-hy*qi[j+1])*Z;
- P[3][2] += (P[3][3].z-hy*qp[j+1])*Z;
- // z: second derivative
- P[1][1] += (P[0][1].z+P[1][0].z-P[0][0].z+hxy*ri[j])*Z;
- P[1][2] += (P[0][2].z+P[1][3].z-P[0][3].z-hxy*ri[j+1])*Z;
- P[2][1] += (P[2][0].z+P[3][1].z-P[3][0].z-hxy*rp[j])*Z;
- P[2][2] += (P[2][3].z+P[3][2].z-P[3][3].z+hxy*rp[j+1])*Z;
- s.s[++k]=patch(P);
- indexi[j]=k;
+ real zij=zi[j];
+ real zip=zi[jp];
+ real zpj=zp[j];
+ real zpp=zp[jp];
+ real pij=hx*pi[j];
+ real ppj=hx*pp[j];
+ real qip=hy*qi[jp];
+ real qpp=hy*qp[jp];
+ real zippip=zip+hx*pi[jp];
+ real zppmppp=zpp-hx*pp[jp];
+ real zijqij=zij+hy*qi[j];
+ real zpjqpj=zpj+hy*qp[j];
+
+ s[k]=new real[][] {{zij,zijqij,zip-qip,zip},
+ {zij+pij,zijqij+pij+hxy*ri[j],
+ zippip-qip-hxy*ri[jp],zippip},
+ {zpj-ppj,zpjqpj-ppj-hxy*rp[j],
+ zppmppp-qpp+hxy*rp[jp],zppmppp},
+ {zpj,zpjqpj,zpp-qpp,zpp}};
+ ++k;
}
}
}
@@ -150,11 +138,11 @@ surface bispline(real[][] z, real[][] p, real[][] q, real[][] r,
return s;
}
-// return the surface described by a real matrix f, interpolated with
+// return the surface values described by a real matrix f, interpolated with
// xsplinetype and ysplinetype.
-surface surface(real[][] f, real[] x, real[] y,
- splinetype xsplinetype=null, splinetype ysplinetype=xsplinetype,
- bool[][] cond={})
+real[][][] bispline(real[][] f, real[] x, real[] y,
+ splinetype xsplinetype=null,
+ splinetype ysplinetype=xsplinetype, bool[][] cond={})
{
real epsilon=sqrtEpsilon*norm(y);
if(xsplinetype == null)
@@ -176,9 +164,7 @@ surface surface(real[][] f, real[] x, real[] y,
real[][] p=transpose(tp);
for(int i=0; i < n; ++i)
r[i]=clamped(d1[i],d2[i])(y,p[i]);
- surface s=bispline(f,p,q,r,x,y,cond);
- if(xsplinetype == periodic) s.ucyclic(true);
- if(ysplinetype == periodic) s.vcyclic(true);
+ real[][][] s=bispline0(f,p,q,r,x,y,cond);
return s;
}
@@ -248,11 +234,11 @@ surface surface(triple f(pair z), real[] u, real[] v,
vperiodic(fz) ? periodic : notaknot};
} else if(vsplinetype.length != 3) abort("vsplinetype must have length 3");
- surface sx=surface(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
- surface sy=surface(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
- surface sz=surface(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
+ real[][][] sx=bispline(fx,ipt,jpt,usplinetype[0],vsplinetype[0],active);
+ real[][][] sy=bispline(fy,ipt,jpt,usplinetype[1],vsplinetype[1],active);
+ real[][][] sz=bispline(fz,ipt,jpt,usplinetype[2],vsplinetype[2],active);
- surface s=surface(sx.s.length);
+ surface s=surface(sx.length);
s.index=new int[nu][nv];
int k=-1;
for(int i=0; i < nu; ++i) {
@@ -261,12 +247,20 @@ surface surface(triple f(pair z), real[] u, real[] v,
indexi[j]=++k;
}
- for(int k=0; k < sx.s.length; ++k) {
+ for(int k=0; k < sx.length; ++k) {
triple[][] Q=new triple[4][];
- for(int i=0; i < 4 ; ++i)
- Q[i]=sequence(new triple(int j) {
- return (sx.s[k].P[i][j].z,sy.s[k].P[i][j].z,sz.s[k].P[i][j].z);
- },4);
+ real[][] Px=sx[k];
+ real[][] Py=sy[k];
+ real[][] Pz=sz[k];
+ for(int i=0; i < 4 ; ++i) {
+ real[] Pxi=Px[i];
+ real[] Pyi=Py[i];
+ real[] Pzi=Pz[i];
+ Q[i]=new triple[] {(Pxi[0],Pyi[0],Pzi[0]),
+ (Pxi[1],Pyi[1],Pzi[1]),
+ (Pxi[2],Pyi[2],Pzi[2]),
+ (Pxi[3],Pyi[3],Pzi[3])};
+ }
s.s[k]=patch(Q);
}
@@ -282,8 +276,8 @@ surface surface(triple f(pair z), real[] u, real[] v,
path3 interp(path3 a, path3 b, real t)
{
int n=size(a);
- return path3(sequence(new triple(int i) {return interp(precontrol(a,i),
- precontrol(b,i),t);},n),
+ return path3(sequence(new triple(int i) {
+ return interp(precontrol(a,i),precontrol(b,i),t);},n),
sequence(new triple(int i) {return interp(point(a,i),point(b,i),t);},n),
sequence(new triple(int i) {return interp(postcontrol(a,i),
postcontrol(b,i),t);},n),
@@ -294,12 +288,15 @@ path3 interp(path3 a, path3 b, real t)
struct tube
{
surface s;
- path3 center;
+ path3 center; // tube axis
- void operator init(path3 p, real width, int sectors=4,
- real granularity=linegranularity) {
- sectors += sectors % 2; // Must be even.
- int h=quotient(sectors,2);
+ void Null(transform3) {}
+ void Null(transform3, bool) {}
+
+ void operator init(path3 p, real width, render render=defaultrender,
+ void cylinder(transform3)=Null,
+ void sphere(transform3, bool half)=Null,
+ void tube(path3, path3)=null) {
real r=0.5*width;
void generate(path3 p) {
@@ -308,31 +305,41 @@ struct tube
for(int i=0; i < n; ++i) {
triple v=point(p,i);
triple u=point(p,i+1)-v;
- s.append(shift(v)*align(unit(u))*scale(r,r,abs(u))*unitcylinder);
+ transform3 t=shift(v)*align(unit(u))*scale(r,r,abs(u));
+ s.append(t*unitcylinder);
+ cylinder(t);
}
center=center&p;
} else {
real[] T;
path3 G;
for(int i=0; i < n; ++i)
- render(subpath(p,i,i+1),granularity,
+ render(subpath(p,i,i+1),
new void(path3 g, real s) {
G=G&g;
T.push(i+s);
- });
+ },render);
T.push(n);
T.cyclic=cyclic(p);
rmf[] rmf=rmf(p,T);
triple f(pair t) {
rmf R=rmf[round(t.x)];
- return point(G,t.x)+r*(R.r*cos(t.y)-R.s*sin(t.y));
+ int n=round(t.y);
+ static real[] x={1,0,-1,0};
+ static real[] y={0,1,0,-1};
+ return point(G,t.x)+r*(R.r*x[n]-R.s*y[n]);
}
- real[] v=uniform(0,2pi,sectors);
+ static real[] v={0,1,2,3,0};
+ static real[] circular(real[] x, real[] y) {
+ static real a=8/3*(sqrt(2)-1);
+ return a*periodic(x,y);
+ }
+
static splinetype[] Monotonic={monotonic,monotonic,monotonic};
- static splinetype[] Periodic={periodic,periodic,periodic};
+ static splinetype[] Circular={circular,circular,circular};
if(T.length > 0) {
- surface S=surface(f,sequence(T.length),v,Monotonic,Periodic);
+ surface S=surface(f,sequence(T.length),v,Monotonic,Circular);
s.append(S);
// Compute center of tube:
@@ -341,19 +348,52 @@ struct tube
triple[] pre=new triple[n+1];
triple[] point=new triple[n+1];
triple[] post=new triple[n+1];
+
int[] index=S.index[0];
- pre[0]=point[0]=0.5*(S.s[index[0]].P[0][0]+S.s[index[h]].P[0][0]);
+ triple Point;
+ for(int m=0; m < 4; ++m)
+ Point += S.s[index[m]].P[0][0];
+ pre[0]=point[0]=0.25*Point;
+
for(int i=0; i < n; ++i) {
index=S.index[i];
- triple [][] P=S.s[index[0]].P;
- triple [][] Q=S.s[index[h]].P;
- post[i]=0.5*(P[1][0]+Q[1][0]);
- pre[i+1]=0.5*(P[2][0]+Q[2][0]);
- point[i+1]=0.5*(P[3][0]+Q[3][0]);
+ triple Pre,Point,Post;
+ for(int m=0; m < 4; ++m) {
+ triple [][] P=S.s[index[m]].P;
+ Post += P[1][0];
+ Pre += P[2][0];
+ Point += P[3][0];
+ }
+ post[i]=0.25*Post;
+ pre[i+1]=0.25*Pre;
+ point[i+1]=0.25*Point;
+
}
+
index=S.index[n-1];
- post[n]=0.5*(S.s[index[0]].P[3][0]+S.s[index[h]].P[3][0]);
- center=center&path3(pre,point,post,array(n+1,false),T.cyclic);
+ triple Post;
+ for(int m=0; m < 4; ++m)
+ Post += S.s[index[m]].P[3][0];
+ post[n]=0.25*Post;
+
+ bool[] b=array(n+1,false);
+ path3 Center=path3(pre,point,post,b,T.cyclic);
+ center=center&Center;
+
+ if(tube != null) { // Compute path along tube
+ triple[] pre=new triple[n+1];
+ triple[] point=new triple[n+1];
+ triple[] post=new triple[n+1];
+ pre[0]=point[0]=S.s[S.index[0][0]].P[0][0];
+ for(int i=0; i < n; ++i) {
+ triple [][] P=S.s[S.index[i][0]].P;
+ post[i]=P[1][0];
+ pre[i+1]=P[2][0];
+ point[i+1]=P[3][0];
+ }
+ post[n]=S.s[S.index[n-1][0]].P[3][0];
+ tube(Center,path3(pre,point,post,b,T.cyclic));
+ }
}
}
}
@@ -363,11 +403,16 @@ struct tube
int begin=0;
int n=length(p);
for(int i=cyclic ? 0 : 1; i < n; ++i)
- if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) {
- generate(subpath(p,begin,i));
- s.append(shift(point(p,i))*t*align(dir(p,i,-1))*unithemisphere);
- begin=i;
- }
+ if(abs(dir(p,i,1)-dir(p,i,-1)) > sqrtEpsilon) {
+ generate(subpath(p,begin,i));
+ triple dir=dir(p,i,-1);
+ s.append(shift(point(p,i))*t*align(dir)*
+ (dir != O ? unithemisphere : unitsphere));
+ int L=length(center);
+ sphere(shift(point(center,L))*t*align(dir(center,L,-1)),
+ half=straight(p,i-1) && straight(p,i));
+ begin=i;
+ }
generate(subpath(p,begin,n));
}
}