Asymptote: Parametric surface isn't working with Spline interpolation

Question

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));

Answer 1

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.