Asymptote: Parametric surface isn't working with Spline interpolation

I have a parametric surface that is defined in terms of the Lambert W function and which I want to display with Asymptote. For that the Lambert W function was implemented using Newton's method and the original (closed) surface had to be split up into two open surfaces to avoid divide by zero issues. Here is the MWE:

``````settings.render=8;
settings.prc=false;
settings.outformat="pdf";
import graph3;

size(200);

currentprojection=orthographic(40,10,10);

// contour value
real c = 0.006;

// parameter ranges
real umin = asin(1.5*exp(1)*sqrt(c*(sqrt(2*pi))));
real umax = pi-asin(1.5*exp(1)*sqrt(c*(sqrt(2*pi))));
real vmin = 0;
real vmax = 1;

// Lambert W function
real w1(real w, real z, int i){return z<-1/exp(1) - 0.00001 ? 1/0 : z<-1/exp(1) ? -1 : i>0 && abs((w*exp(w)-z)/(exp(w)+w*exp(w))) > 1e-7 ? w1(w-(w*exp(w)-z)/(exp(w)+w*exp(w)),z,i-1) : w-(w*exp(w)-z)/(exp(w)+w*exp(w));};

// auxiliary functions
real y5(real h, real p){return (1/4.)*sqrt(15./pi) * sin(2*p) * sin(h)**2;};

real r1(real y){return -6*w1(-2,-sqrt(c*9*sqrt(30)/abs(y))/4,200);};
real r2(real y){return -6*w1(1,-sqrt(c*9*sqrt(30)/abs(y))/4,200);};

// x, y, and z coordinates of the surfaces
real x11(real u, real v){return r1(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2))))*sin(u)*cos(v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)));};
real y11(real u, real v){return r1(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2))))*sin(u)*sin(v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)));};
real z11(real u, real v){return r1(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)) ))*cos(u);};

real x12(real u, real v){return r2(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2))))*sin(u)*cos(v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)));};
real y12(real u, real v){return r2(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2))))*sin(u)*sin(v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)));};
real z12(real u, real v){return r2(y5(u,v*pi/2+(0.5-v)*asin(exp(1)**2*9.*c*sqrt(2.*pi)/(4.*sin(u)**2)) ))*cos(u);};

triple f11(pair p){return (x11(p.x,p.y),y11(p.x,p.y),z11(p.x,p.y));};
triple f12(pair p){return (x12(p.x,p.y),y12(p.x,p.y),z12(p.x,p.y));};

surface s11 = surface(f=f11,a=(umin, vmin),b=(umax,vmax));  // this works
surface s12 = surface(f=f12,a=(umin, vmin),b=(umax,vmax));  // this works

// surface s11 = surface(f=f11,a=(umin, vmin),b=(umax,vmax),Spline);  // this doesn't work
// surface s12 = surface(f=f12,a=(umin, vmin),b=(umax,vmax),Spline);  // this doesn't work

draw(s11, red+opacity(0.5));
draw(s12, red+opacity(0.5));
``````

If I add the `Spline` directive to make the surface look smooth Asymptote crashes with the error `some path/graph_splinetype.asy: 89.10: function values are not periodic`. I tried to understand what went wrong in `graph_splinetype.asy` and `graph3.asy` but unfortunately I'm not proficient enough to succeed. So my question is: Is there a chance to get `Spline` working with this parametric surface or maybe another way to make it look smooth?

What makes it even more puzzling is that `Spline` works just fine for a similar parametric surface (although this one has rotational symmetry around the z axis which might be important), namely this one:

``````settings.render=8;
settings.prc=false;
settings.outformat="pdf";
import graph3;

size(200);

currentprojection=orthographic(40,10,10);

// contour value
real c = 0.006;

// parameter ranges
real umin = 0;
real umax = acos(2*exp(1)*c*sqrt(2*pi))-0.0000001;
real vmin = 0;
real vmax = 2*pi;

// Lambert W function
real w1(real w, real z, int i){return z<-1/exp(1) - 0.00001 ? 1/0 : z<-1/exp(1) ? -1 : i>0 && abs((w*exp(w)-z)/(exp(w)+w*exp(w))) > 1e-7 ? w1(w-(w*exp(w)-z)/(exp(w)+w*exp(w)),z,i-1) : w-(w*exp(w)-z)/(exp(w)+w*exp(w));};

// auxiliary functions
real y1(real h){return sqrt(3./pi)/2.*cos(h);};

real r1(real y){return -2*w1(-2,-c*sqrt(6)/abs(y),200);};
real r2(real y){return -2*w1(1,-c*sqrt(6)/abs(y),200);};

// x, y, and z coordinates of the surfaces
real x11(real u, real v){return r1(y1(u))*sin(u)*cos(v);};
real y11(real u, real v){return r1(y1(u))*sin(u)*sin(v);};
real z11(real u, real v){return r1(y1(u))*cos(u);};

real x12(real u, real v){return r2(y1(u))*sin(u)*cos(v);};
real y12(real u, real v){return r2(y1(u))*sin(u)*sin(v);};
real z12(real u, real v){return r2(y1(u))*cos(u);};

triple f11(pair p){return (x11(p.x,p.y),y11(p.x,p.y),z11(p.x,p.y));};
triple f12(pair p){return (x12(p.x,p.y),y12(p.x,p.y),z12(p.x,p.y));};

surface s11 = surface(f11,(umin,vmin),(umax,vmax),50,Spline);
surface s12 = surface(f12,(umin,vmin),(umax,vmax),50,Spline);

draw(s11, red+opacity(0.5));
draw(s12, red+opacity(0.5));
``````

Could you precise your system and Asymptote version ? On my computer (Linux 64 bits, Sid, Asymptote svn) I have not problem and your code produces a smooth surface.

Since the choice between periodic or not_a_knot Spline is implemented in the surface routine it seems that it is less restrictive than the final parametric spline routine (I wrote a long time ago a part of this smooth surface routine). It is possible to force a "not a knot" choice with the following code

``````splinetype[] Notaknot={notaknot,notaknot,notaknot};
surface s11 = surface(f=f11,a=(umin,vmin),b=(umax,vmax),8,16,Notaknot,Notaknot);
surface s12 = surface(f=f12,a=(umin,vmin),b=(umax,vmax),8,16,Notaknot,Notaknot);
``````

I hope it works on your system.

O.G.

• Thank you very much for looking into my problem and your answer. Your code snippet works (after I corrected some typos). Sorry for not providing information about my system and Asymptote version. I work on a Windows 7 (64bit) system with Asymptote 2.29. – Philipp Oct 14 '14 at 20:12
• By the way, do you happen to know if there is a way to combine the two surfaces into one object? Because right now the transparancy doesn't look good with the two surfaces handled seperately. – Philipp Oct 14 '14 at 20:15
• surface s; for(int i=0;i<s11.s.length;++i) { s.push(s11.s[i]); } for(int i=0;i<s12.s.length;++i) { s.push(s12.s[i]); } – O.G. Oct 14 '14 at 20:24
• I just installed Asymptote 2.32 (the most recent Windows version I could find) on my computer. It shows the same result: Doesn't work with my original code but works with your code snippet included. – Philipp Oct 14 '14 at 20:25
• Wow, that was fast. Thanks again, that worked very well too. – Philipp Oct 14 '14 at 20:29