diff options
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r-- | Build/source/utils/asymptote/base/ode.asy | 255 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/plain_picture.asy | 22 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three.asy | 56 | ||||
-rw-r--r-- | Build/source/utils/asymptote/base/three_surface.asy | 22 |
4 files changed, 252 insertions, 103 deletions
diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy index fe65f2f1c8e..77ab8d6e589 100644 --- a/Build/source/utils/asymptote/base/ode.asy +++ b/Build/source/utils/asymptote/base/ode.asy @@ -1,7 +1,8 @@ -real stepfactor=2.0; // Maximum dynamic step size adjustment factor. +real stepfactor=2; // Maximum dynamic step size adjustment factor. struct coefficients { + real[] steps; real[] factors; real[][] weights; real[] highOrderWeights; @@ -11,7 +12,6 @@ struct coefficients struct RKTableau { int order; - real[] steps; coefficients a; void stepDependence(real h, real c, coefficients a) {} @@ -25,8 +25,8 @@ struct RKTableau return sum(weights[i]);},weights.length), void stepDependence(real, real, coefficients)=null) { this.order=order; - this.steps=steps; - a.factors=array(steps.length+1,1); + a.steps=steps; + a.factors=array(a.steps.length+1,1); a.weights=weights; a.highOrderWeights=highOrderWeights; a.lowOrderWeights=lowOrderWeights; @@ -39,24 +39,113 @@ struct RKTableau } } +real[] Coeff={1,1/2,1/6,1/24,1/120,1/720,1/5040,1/40320,1/362880,1/3628800, + 1/39916800.0,1/479001600.0,1/6227020800.0,1/87178291200.0, + 1/1307674368000.0,1/20922789888000.0,1/355687428096000.0, + 1/6402373705728000.0,1/121645100408832000.0, + 1/2432902008176640000.0,1/51090942171709440000.0, + 1/1124000727777607680000.0}; + +real phi1(real x) {return x != 0 ? expm1(x)/x : 1;} + +real phi2(real x) +{ + real x2=x*x; + if(fabs(x) > 1) return (exp(x)-x-1)/x2; + real x3=x2*x; + real x5=x2*x3; + if(fabs(x) < 0.1) + return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5] + +x5*Coeff[6]+x3*x3*Coeff[7]+x5*x2*Coeff[8]+x5*x3*Coeff[9]; + else { + real x7=x5*x2; + real x8=x7*x; + return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5] + +x5*Coeff[6]+x3*x3*Coeff[7]+x7*Coeff[8]+x8*Coeff[9] + +x8*x*Coeff[10]+x5*x5*Coeff[11]+x8*x3*Coeff[12]+x7*x5*Coeff[13]+ + x8*x5*Coeff[14]+x7*x7*Coeff[15]+x8*x7*Coeff[16]+x8*x8*Coeff[17]; + } +} + +real phi3(real x) +{ + real x2=x*x; + real x3=x2*x; + if(fabs(x) > 1.6) return (exp(x)-0.5*x2-x-1)/x3; + real x5=x2*x3; + if(fabs(x) < 0.1) + return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5] + +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9] + +x5*x3*Coeff[10]; + else { + real x7=x5*x2; + real x8=x7*x; + real x16=x8*x8; + return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5] + +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9] + +x5*x3*Coeff[10]+x8*x*Coeff[11] + +x5*x5*Coeff[12]+x8*x3*Coeff[13]+x7*x5*Coeff[14] + +x8*x5*Coeff[15]+x7*x7*Coeff[16]+x8*x7*Coeff[17]+x16*Coeff[18] + +x16*x*Coeff[19]+x16*x2*Coeff[20]; + } +} + +void expfactors(real x, coefficients a) +{ + for(int i=0; i < a.steps.length; ++i) + a.factors[i]=exp(x*a.steps[i]); + a.factors[a.steps.length]=exp(x); +} + // First-Order Euler RKTableau Euler=RKTableau(1,new real[][], new real[] {1}); +// First-Order Exponential Euler RKTableau E_Euler=RKTableau(1,new real[][], new real[] {1}, - new void (real h, real c, coefficients a) { + new void(real h, real c, coefficients a) { real x=-c*h; - a.factors[0]=exp(x); - a.highOrderWeights[0]=x != 0 ? expm1(x)/x : 1; + expfactors(x,a); + a.highOrderWeights[0]=phi1(x); }); // Second-Order Runge-Kutta RKTableau RK2=RKTableau(2,new real[][] {{1/2}}, - new real[] {0,1}); + new real[] {0,1}, // 2nd order + new real[] {1,0}); // 1st order + +// Second-Order Exponential Runge-Kutta +RKTableau E_RK2=RKTableau(2,new real[][] {{1/2}}, + new real[] {0,1}, // 2nd order + new real[] {1,0}, // 1st order + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + a.weights[0][0]=1/2*phi1(x/2); + real w=phi1(x); + a.highOrderWeights[0]=0; + a.highOrderWeights[1]=w; + a.lowOrderWeights[0]=w; + }); // Second-Order Predictor-Corrector RKTableau PC=RKTableau(2,new real[][] {{1}}, - new real[] {1/2,1/2}); + new real[] {1/2,1/2}, // 2nd order + new real[] {1,0}); // 1st order + +// Second-Order Exponential Predictor-Corrector +RKTableau E_PC=RKTableau(2,new real[][] {{1}}, + new real[] {1/2,1/2}, // 2nd order + new real[] {1,0}, // 1st order + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + real w=phi1(x); + a.weights[0][0]=w; + a.highOrderWeights[0]=w/2; + a.highOrderWeights[1]=w/2; + a.lowOrderWeights[0]=w; + }); // Third-Order Classical Runge-Kutta RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}}, @@ -67,6 +156,30 @@ RKTableau RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}}, new real[] {2/9,1/3,4/9}, // 3rd order new real[] {7/24,1/4,1/3,1/8}); // 2nd order +// Third-Order Exponential Bogacki-Shampine Runge-Kutta +RKTableau E_RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}}, + new real[] {2/9,1/3,4/9}, // 3rd order + new real[] {7/24,1/4,1/3,1/8}, // 2nd order + new void(real h, real c, coefficients a) { + real x=-c*h; + expfactors(x,a); + real w=phi1(x); + real w2=phi2(x); + a.weights[0][0]=1/2*phi1(x/2); + real a11=9/8*phi2(3/4*x)+3/8*phi2(x/2); + a.weights[1][0]=3/4*phi1(3/4*x)-a11; + a.weights[1][1]=a11; + real a21=1/3*w; + real a22=4/3*w2-2/9*w; + a.highOrderWeights[0]=w-a21-a22; + a.highOrderWeights[1]=a21; + a.highOrderWeights[2]=a22; + a.lowOrderWeights[0]=w-17/12*w2; + a.lowOrderWeights[1]=w2/2; + a.lowOrderWeights[2]=2/3*w2; + a.lowOrderWeights[3]=w2/4; + }); + // Fourth-Order Classical Runge-Kutta RKTableau RK4=RKTableau(4,new real[][] {{1/2},{0,1/2},{0,0,1}}, new real[] {1/6,1/3,1/3,1/6}); @@ -108,7 +221,7 @@ RKTableau RK5DP=RKTableau(5,new real[][] {{1/5}, real error(real error, real initial, real lowOrder, real norm, real diff) { - if(initial != 0.0 && lowOrder != initial) { + if(initial != 0 && lowOrder != initial) { static real epsilon=realMin/realEpsilon; real denom=max(abs(norm),abs(initial))+epsilon; return max(error,max(abs(diff)/denom)); @@ -116,46 +229,39 @@ real error(real error, real initial, real lowOrder, real norm, real diff) return error; } -real adjust(real h, real error, real t, real c, real tolmin, real tolmax, - real dtmin, real dtmax, RKTableau tableau, bool verbose=true) +void report(real old, real h, real t) { - real dt=h; - void report(real t) { - if(h != dt) { - tableau.stepDependence(h,c,tableau.a); - if(verbose) - write("Time step changed from "+(string) dt+" to "+(string) h+" at t="+ - (string) t+"."); - } - } - if(error > tolmax) { - h=max(h*max((tolmin/error)^tableau.pshrink,1/stepfactor),dtmin); - report(t); - return h; - } - if(error > 0 && error < tolmin) { - h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax); - report(t+dt); - } + write("Time step changed from "+(string) old+" to "+(string) h+" at t="+ + (string) t+"."); +} + +real adjust(real h, real error, real tolmin, real tolmax, RKTableau tableau) +{ + if(error > tolmax) + h *= max((tolmin/error)^tableau.pshrink,1/stepfactor); + else if(error > 0 && error < tolmin) + h *= min((tolmin/error)^tableau.pgrow,stepfactor); return h; } // Integrate dy/dt+cy=f(t,y) from a to b using initial conditions y, // specifying either the step size h or the number of steps n. -real integrate(real y, real c=0, real g(real t, real y), real a, real b=a, - real h=0, int n=0, bool dynamic=false, real tolmin=0, - real tolmax=0, real dtmin=0, real dtmax=realMax, - RKTableau tableau, bool verbose=false) +real[] integrate(real y, real c=0, real g(real t, real y), real a, real b=a, + real h=0, int n=0, bool dynamic=false, real tolmin=0, + real tolmax=0, real dtmin=0, real dtmax=realMax, + RKTableau tableau, bool verbose=false) { - real f(real t, real y)=(c == 0 || tableau.exponential) ? g : - new real(real t, real y) {return g(t,y)-c*y;}; + real[] Y={y}; if(h == 0) { - if(b == a) return y; + if(b == a) return Y; if(n == 0) abort("Either n or h must be specified"); else h=(b-a)/n; } + real f(real t, real y)=(c == 0 || tableau.exponential) ? g : + new real(real t, real y) {return g(t,y)-c*y;}; + tableau.stepDependence(h,c,tableau.a); real t=a; @@ -165,58 +271,64 @@ real integrate(real y, real c=0, real g(real t, real y), real a, real b=a, (tableau.a.lowOrderWeights.length > tableau.a.highOrderWeights.length); if(fsal) f0=f(t,y); + real dt=h; while(t < b) { + h=min(h,b-t); + if(t+h == t) break; + if(h != dt) { + if(verbose) report(dt,h,t); + tableau.stepDependence(h,c,tableau.a); + dt=h; + } + real[] predictions={fsal ? f0 : f(t,y)}; - for(int i=0; i < tableau.steps.length; ++i) - predictions.push(f(t+h*tableau.steps[i], + for(int i=0; i < tableau.a.steps.length; ++i) + predictions.push(f(t+h*tableau.a.steps[i], tableau.a.factors[i]*y+h*dot(tableau.a.weights[i], predictions))); real highOrder=h*dot(tableau.a.highOrderWeights,predictions); - real Y=tableau.a.factors[tableau.steps.length]*y; + real y0=tableau.a.factors[tableau.a.steps.length]*y; if(dynamic) { real f1; if(fsal) { - f1=f(t+h,Y+highOrder); + f1=f(t+h,y0+highOrder); predictions.push(f1); } real lowOrder=h*dot(tableau.a.lowOrderWeights,predictions); real error; - error=error(error,y,Y+lowOrder,Y+highOrder,highOrder-lowOrder); - real dt=h; - h=adjust(h,error,t,c,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose); + error=error(error,y,y0+lowOrder,y0+highOrder,highOrder-lowOrder); + h=adjust(h,error,tolmin,tolmax,tableau); if(h >= dt) { t += dt; - y=Y+highOrder; + y=y0+highOrder; + Y.push(y); f0=f1; } + h=min(max(h,dtmin),dtmax); } else { t += h; - y=Y+highOrder; - } - real remain=b-t; - if(h > remain) { - h=remain; - tableau.stepDependence(h,c,tableau.a); + y=y0+highOrder; + Y.push(y); } - if(t >= b || t+h == t) break; } - return y; + return Y; } // Integrate a set of equations, dy/dt=f(t,y), from a to b using initial // conditions y, specifying either the step size h or the number of steps n. -real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a, - real h=0, int n=0, bool dynamic=false, - real tolmin=0, real tolmax=0, real dtmin=0, real dtmax=realMax, - RKTableau tableau, bool verbose=false) +real[][] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a, + real h=0, int n=0, bool dynamic=false, + real tolmin=0, real tolmax=0, real dtmin=0, + real dtmax=realMax, RKTableau tableau, bool verbose=false) { + real[][] Y={copy(y)}; + if(h == 0) { - if(b == a) return y; + if(b == a) return Y; if(n == 0) abort("Either n or h must be specified"); else h=(b-a)/n; } - real[] y=copy(y); real t=a; real[] f0; if(tableau.a.lowOrderWeights.length == 0) dynamic=false; @@ -224,10 +336,18 @@ real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a, (tableau.a.lowOrderWeights.length > tableau.a.highOrderWeights.length); if(fsal) f0=f(t,y); + real dt=h; while(t < b) { + h=min(h,b-t); + if(t+h == t) break; + if(h != dt) { + if(verbose) report(dt,h,t); + dt=h; + } + real[][] predictions={fsal ? f0 : f(t,y)}; - for(int i=0; i < tableau.steps.length; ++i) - predictions.push(f(t+h*tableau.steps[i], + for(int i=0; i < tableau.a.steps.length; ++i) + predictions.push(f(t+h*tableau.a.steps[i], y+h*tableau.a.weights[i]*predictions)); real[] highOrder=h*tableau.a.highOrderWeights*predictions; @@ -242,21 +362,21 @@ real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a, for(int i=0; i < y.length; ++i) error=error(error,y[i],y[i]+lowOrder[i],y[i]+highOrder[i], highOrder[i]-lowOrder[i]); - real dt=h; - h=adjust(h,error,t,0,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose); + h=adjust(h,error,tolmin,tolmax,tableau); if(h >= dt) { t += dt; y += highOrder; + Y.push(y); f0=f1; } + h=min(max(h,dtmin),dtmax); } else { t += h; y += highOrder; + Y.push(y); } - h=min(h,b-t); - if(t >= b || t+h == t) break; } - return y; + return Y; } real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t, @@ -292,7 +412,8 @@ real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0, real[] guess, RKTableau tableau, int iterations=100) { real[] g(real[] t) { - return discrepancy(integrate(initial(t),f,a,b,h,n,tableau)); + real[][] y=integrate(initial(t),f,a,b,h,n,tableau); + return discrepancy(y[y.length-1]); } real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);} return initial(newton(iterations,g,jacobian,guess)); diff --git a/Build/source/utils/asymptote/base/plain_picture.asy b/Build/source/utils/asymptote/base/plain_picture.asy index 4d894dc6dea..b3836d7ef51 100644 --- a/Build/source/utils/asymptote/base/plain_picture.asy +++ b/Build/source/utils/asymptote/base/plain_picture.asy @@ -408,16 +408,18 @@ struct projection { bool infinity; bool oblique; bool absolute=false; - triple camera; - triple up; - triple target; + triple camera; // Position of camera. + triple up; // A vector that should be projected to direction (0,1). + triple target; // Point where camera is looking at. + pair viewportshift; // Fractional viewport shift. + real zoom=1; // Zoom factor. + real angle; // Lens angle (for perspective projection). bool showtarget=true; // Expand bounding volume to include target? typedef transformation projector(triple camera, triple up, triple target); projector projector; bool autoadjust=true; // Adjust camera to lie outside bounding volume? - bool center=false; // Center target within bounding volume? - real angle; // Lens angle (currently only used by PRC viewpoint). - int ninterpolate; // Used for projecting nurbs to 2D Bezier curves. + bool center=false; // Center target within bounding volume? + int ninterpolate; // Used for projecting nurbs to 2D Bezier curves. void calculate() { transformation T=projector(camera,up,target); @@ -438,11 +440,15 @@ struct projection { } void operator init(triple camera, triple up=(0,0,1), triple target=(0,0,0), + real zoom=1, real angle=0, pair viewportshift=0, bool showtarget=true, bool autoadjust=true, bool center=false, projector projector) { this.camera=camera; this.up=up; this.target=target; + this.zoom=zoom; + this.angle=angle; + this.viewportshift=viewportshift; this.showtarget=showtarget; this.autoadjust=autoadjust; this.center=center; @@ -459,11 +465,13 @@ struct projection { P.camera=camera; P.up=up; P.target=target; + P.zoom=zoom; + P.angle=angle; + P.viewportshift=viewportshift; P.showtarget=showtarget; P.autoadjust=autoadjust; P.center=center; P.projector=projector; - P.angle=angle; P.ninterpolate=ninterpolate; return P; } diff --git a/Build/source/utils/asymptote/base/three.asy b/Build/source/utils/asymptote/base/three.asy index d42ea928b16..63a497c497c 100644 --- a/Build/source/utils/asymptote/base/three.asy +++ b/Build/source/utils/asymptote/base/three.asy @@ -17,7 +17,6 @@ real linegranularity=0.01; real tubegranularity=0.003; real dotgranularity=0.0001; real viewportfactor=1.002; // Factor used to expand orthographic viewport. -real viewportpadding=1; // Offset used to expand PRC viewport. real angleprecision=1e-3; // Precision for centering perspective projections. real anglefactor=max(1.01,1+angleprecision); // Factor used to expand perspective viewport. @@ -127,10 +126,6 @@ transform3 reflect(triple u, triple v, triple w) }*shift(-u); } -bool operator != (real[][] a, real[][] b) { - return !(a == b); -} - // Project u onto v. triple project(triple u, triple v) { @@ -219,38 +214,44 @@ pair dir(triple v, triple dir, projection P) // 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,showtarget,autoadjust,center, + 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,showtarget,autoadjust,center); + 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,showtarget,center=center, + 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, bool showtarget=true, - bool center=false) + triple target=O, real zoom=1, pair viewportshift=0, + bool showtarget=true, bool center=false) { - return orthographic((x,y,z),up,target,showtarget,center=center); + return orthographic((x,y,z),up,target,zoom,viewportshift,showtarget, + center=center); } projection oblique(real angle=45) @@ -2278,20 +2279,27 @@ string embed3D(string label="", string text=label, string prefix, real viewplanesize; if(P.infinity) { triple lambda=max3(f)-min3(f); - pair margin=viewportpadding*(1,1)+viewportmargin((lambda.x,lambda.y)); - viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/cm; + pair margin=viewportmargin((lambda.x,lambda.y)); + viewplanesize=(max(lambda.x+2*margin.x,lambda.y+2*margin.y))/(cm*P.zoom); } else - if(!P.absolute) angle=2*aTan(Tan(0.5*angle)-viewportpadding/P.target.z); + if(!P.absolute) angle=2*aTan(Tan(0.5*angle)); string name=prefix+".js"; writeJavaScript(name,lightscript+projection(P.infinity,viewplanesize),script); shipout3(prefix,f); - + prefix += ".prc"; if(!settings.inlinetex) file3.push(prefix); + triple target=P.target; + if(P.infinity && 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()/cm; triple u=unit(v); triple w=Z-u.z*u; @@ -2310,7 +2318,7 @@ string embed3D(string label="", string text=label, string prefix, ",toolbar="+(settings.toolbar ? "true" : "false")+ ",3Daac="+Format(P.absolute ? P.angle : angle)+ ",3Dc2c="+Format(u)+ - ",3Dcoo="+Format(P.target/cm)+ + ",3Dcoo="+Format(target/cm)+ ",3Droll="+Format(roll)+ ",3Droo="+Format(abs(v))+ ",3Dbg="+Format(light.background()); @@ -2437,9 +2445,7 @@ object embed(string label="", string text=label, viewportmargin=viewportmargin((lambda.x,lambda.y)); width=lambda.x+2*viewportmargin.x; height=lambda.y+2*viewportmargin.y; - - triple s=(-0.5(m.x+M.x),-0.5*(m.y+M.y),0); - f=shift(s)*f; // Eye will be at (0,0,0). + f=shift((-0.5(m.x+M.x),-0.5*(m.y+M.y),0))*f; // Eye will be at (0,0,0). } else { transform3 T=identity4; // Choose the angle to be just large enough to view the entire image: @@ -2539,9 +2545,10 @@ object embed(string label="", string text=label, M += margin; m -= margin; } else if(M.z >= 0) abort("camera too close"); - + shipout3(prefix,f,preview ? nativeformat() : format, - width,height,P.infinity ? 0 : angle,m,M, + width,height,P.infinity ? 0 : angle,P.zoom,m,M, + prc && !P.infinity ? 0 : P.viewportshift, tinv*inverse(modelview)*shift(0,0,zcenter),light.background(), P.absolute ? (modelview*light).position : light.position, light.diffuse,light.ambient,light.specular, @@ -2557,8 +2564,13 @@ object embed(string label="", string text=label, if(!settings.inlinetex) file3.push(image); image=graphic(image,"hiresbb"); } - if(prc) F.L=embed3D(label,text=image,prefix,f,format, + if(prc) { + if(!P.infinity && P.viewportshift != 0) + write("warning: PRC does not support off-axis projections; use pan instead of shift"); + F.L=embed3D(label,text=image,prefix,f,format, width,height,angle,options,script,light,Q); + } + } if(!is3D) { diff --git a/Build/source/utils/asymptote/base/three_surface.asy b/Build/source/utils/asymptote/base/three_surface.asy index 597104ad202..2f5d80aacbe 100644 --- a/Build/source/utils/asymptote/base/three_surface.asy +++ b/Build/source/utils/asymptote/base/three_surface.asy @@ -731,14 +731,22 @@ struct surface { for(int i=0; i < L; ++i) { path3 h=subpath(g,i,i+1); path3 r=reverse(h); - triple max=max(h); - triple min=min(h); - triple perp=perp(max-c,axis); - real fuzz=epsilon*max(abs(max),abs(min)); - if(abs(perp) < fuzz) - perp=perp(min-c,axis); + path3 H=shift(-c)*h; + real M=0; + triple perp; + void test(real[] t) { + for(int i=0; i < 3; ++i) { + triple v=point(H,t[i]); + triple V=v-dot(v,axis)*axis; + real a=abs(V); + if(a > M) {M=a; perp=V;} + } + } + test(maxtimes(H)); + test(mintimes(H)); + perp=unit(perp); - triple normal=cross(axis,perp); + triple normal=unit(cross(axis,perp)); triple dir(real j) {return Cos(j)*normal-Sin(j)*perp;} real j=angle1; transform3 Tk=T[0]; |