5

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

enter image description here

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

enter image description here

4

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
  • 1
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.