summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/base
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2009-07-03 00:02:22 +0000
committerKarl Berry <karl@freefriends.org>2009-07-03 00:02:22 +0000
commitc20b61928150008d080bab09e967fcaebc93da81 (patch)
treebffa7c8765003c0e8f2a40c72479e5caa6f748bd /Build/source/utils/asymptote/base
parent9ec802a431dfe1b27063bac2f8dd9bc58cf997b0 (diff)
asy 1.79
git-svn-id: svn://tug.org/texlive/trunk@14078 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/base')
-rw-r--r--Build/source/utils/asymptote/base/ode.asy255
-rw-r--r--Build/source/utils/asymptote/base/plain_picture.asy22
-rw-r--r--Build/source/utils/asymptote/base/three.asy56
-rw-r--r--Build/source/utils/asymptote/base/three_surface.asy22
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];