Revision history [back]

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2))$, $-1\leq x,y,z\leq1$;

$-1=x_0<x_1&lt;...<x_l=1$, $-1="y_0&lt;y_1&lt;...&lt;y_m=1$," $-1="z_0&lt;x_1&lt;...&lt;x_n=1$;&lt;/p">

$0=c_0<x_1&lt;...<c_k=1$, where="" $[0,1]$="" is="" the="" range="" of="" values="" taken="" by="" $f$.<="" p="">

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y') sage: def f0(x,y,n): return x, y, exp(-(x2 + y2)) sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2) sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1) sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3) sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4) sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5) sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5) sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5) sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5) sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5) sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5) sage: P=P0+P1+P2+P3+P4 sage: P.show(viewer='jmol') Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$P=P1+\dots+P99$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100 x,y,n=var('x,y,n') def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1))

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5) ; P=sum(Pplot(n), n in range(N)) P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and jmol discussed in the link given above would be solved completely.

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2))$, $$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2))$, $-1\leq x,y,z\leq1$;

$-1=x_0<x_1&lt;...<x_l=1$, $-1="y_0&lt;y_1&lt;...&lt;y_m=1$," $-1="z_0&lt;x_1&lt;...&lt;x_n=1$;&lt;/p">

$0=c_0<x_1&lt;...<c_k=1$, where="" $[0,1]$="" is="" the="" range="" of="" values="" taken="" by="" $f$.<="" x,y,z\leq1;$$

$$-1=x_0<x_1&lt;...<x_l=1, \quad="" -1="y_0&lt;y_1&lt;...&lt;y_m=1," \quad="" -1="z_0&lt;x_1&lt;...&lt;x_n=1;$$&lt;/p">

$$0=c_0<x_1&lt;...<c_k=1,$$ <="" p="">

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y') x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2)) 2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2) 2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1) 2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3) 2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4) 2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5) 2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5) opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5) opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5) opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5) opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5) sage: P=P0+P1+P2+P3+P4 sage: P.show(viewer='jmol') opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$P=P1+\dots+P99$$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100 x,y,n=var('x,y,n') N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1))y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5) ; opacity=0.5);

P=sum(Pplot(n), n in range(N)) range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and jmol discussed in the link given above would be solved completely.

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2))$, $-1\leq x,y,z\leq1;$$

$$-1=x_0<x_1&lt;...<x_l=1, \quad="" -1="y_0&lt;y_1&lt;...&lt;y_m=1," \quad="" -1="z_0&lt;x_1&lt;...&lt;x_n=1;$$&lt;/p">

$$0=c_0<x_1&lt;...<c_k=1,$$ <="" p="">

where $[0,1]$ is the range of values taken by $F$.$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and jmol the use of jmol, as discussed in the link given above above, would be solved in a very satisfactory way.

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and the use of jmol, as discussed in the link given above, would be solved in a very satisfactory way.

I believe that the following (tested) code provides a simple, yet complete answer to my question:

M=5

funcsurfplot=[0 for i in range(M)]

x, y = var ('x, y', domain='real')

n=var('n', domain='integer')

funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), (x, -2-n, 2+n), (y, -2-n, 2+n), color=rainbow(M)[n], opacity=0.3, mesh=True) for n in range(M)] range(M)]

P=sum(funcsurfplot[i] for i in range(M)) P.show(viewer='jmol') x, y = var ('x, y', domain='real') n=var('n', domain='integer') funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), (x, -2-n, 2+n), (y, -2-n, 2+n), color=rainbow(M)[n], opacity=0.3, mesh=True) for n in range(M)] P=sum(funcsurfplot[i] for i in range(M)) P.show(viewer='jmol') range(M))

P.show(viewer='jmol')

Note that in this case we have full control over the domains of each of the functions in the family of surfaces, and that these domains can also depend on the integer parameter of the family. The line """funcsurfplot=[0 for i in range(M)]" might be considered redundant, but I prefer it for creation of array-like list in Python with custom size, because, unlike the use of 'append', this approach allows the easy construction of multidimensional arrays (two dimensional arrays would be needed for visualizing the families of coordinate surfaces in the example I gave in the formulation of the question) -- see, e.g.,

http://www.i-programmer.info/programming/python/3942-arrays-in-python.html

Now all the other cases of the question can be resolved in the same simple uniform way, as follows:

Second, meshes of 3D curves (curvilinear coordinates in 3D and on surfaces):

M=5

curve3dplot=[0 for i in range(M)] # creation of array of size M+1 via initialization

t=var ('t', domain='real')

n=var('n', domain='integer')

curve3dplot=[parametric_plot3d((exp(-(n+t)^2), cos(ntpi), 1/(1+n*t^2)), (t, -1, 1), color=rainbow(M)[n], thickness=4) for n in range(M)]

P=sum(curve3dplot[i] for i in range(M))

P.show(viewer='jmol')

Third, families of parametric surfaces (the surfaces chosen here are Klein bottles of the default type in SAGE - see

http://doc.sagemath.org/html/en/reference/riemannian_geometry/sage/geometry/riemannian_manifolds/surface3d_generators.html

and, in order to introduce dependence on n also in the parametric domains, I have slightly 'unzipped' them):

sage: M=5

sage: parasurf3dplot=[0 for i in range(M)]

sage: u, v = var ('u, v', domain='real')

sage: n=var('n', domain='integer')

sage: def fx(u,v,n): return ((n+1)/M+cos(u/2)cos(v)-sin(u/2)sin(2v))cos(u)

sage: def fy(u,v,n): return ((n+1)/M+cos(u/2)cos(v)-sin(u/2)sin(2v))sin(u)

sage: def fz(u,v,n): return sin(u/2)cos(v)+cos(u/2)sin(2*v)

sage: def rangeu(n): return (u, -pi+(n/M^3)pi,pi-(n/M^3)pi)

sage: def rangev(n): return (v, -pi+(n/M^3)pi,pi-(n/M^3)pi)

sage: parasurf3dplot=[parametric_plot3d((fx(u,v,n), fy(u,v,n), fz(u,v,n)), rangeu(n), rangev(n), color=rainbow(M)[n], opacity=0.1, mesh=True) for n in range(M)]

sage: P=sum(parasurf3dplot[i] for i in range(M))

sage: P.show(viewer='jmol')

Fourth, families of isosurfaces:

M=5

isosurfplot=[0 for i in range(M)]

x, y, z = var ('x, y, z', domain='real')

n=var('n', domain='integer')

isosurfplot=[implicit_plot3d(x^2+y^2+z^2-(n+1)^2, (x, -2-2n/(M-1), 2+2n/(M-1)), (y, -2-2n/(M-1),2+2n/(M-1) ), (z, -2-2n/(M-1),2+2n/(M-1)), color=rainbow(M)[n], opacity=0.3) for n in range(M)]

P=sum(isosurfplot[i] for i in range(M)) # NOT sum([...])

P.show(viewer='jmol')

Enjoy the results and feel very welcome to improve and optimize this code. I suggest that in the next version of the SAGE tutorial and/or reference manual there appears a subtopic dedicated to this, as a 3D-upgrade of the topic about contour plots. If, for this to happen, assistance is required from me, I will provide such assistance.

shao-linux :)

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and the use of jmol, as discussed in the link given above, would be solved in a very satisfactory way.

I believe that the following (tested) code provides a simple, yet complete answer to my question:

First, for functional surfaces (z=f(x,y)), the following code provides an alternative to the solution invoking the class ParametrizedSurface3D: here we only use parametric_plot3d:

M=5
funcsurfplot=[0 for i in range(M)]
 
range(M)]
x, y = var ('x, y', domain='real') 
 
 n=var('n', domain='integer')
 
domain='integer')
funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), (x, -2-n, 2+n), (y, -2-n, 2+n),   \   color=rainbow(M)[n], opacity=0.3,  mesh=True) for n in range(M)]
 
range(M)]
P=sum(funcsurfplot[i] for i in range(M))
 
P.show(viewer='jmol')range(M))
P.show(viewer='jmol')

M=5
curve3dplot=[0 for i in range(M)] # creation of array of size M+1 via initialization
 
initialization
t=var ('t', domain='real')
 
domain='real')
n=var('n', domain='integer')
 
domain='integer')
curve3dplot=[parametric_plot3d((exp(-(n+t)^2), cos(ntpi), cos(n*t*pi), 1/(1+n*t^2)),  (t, -1, 1), \               color=rainbow(M)[n], thickness=4) for n in range(M)]
 
range(M)]
P=sum(curve3dplot[i] for i in range(M))
 
P.show(viewer='jmol')range(M))
P.show(viewer='jmol')

M=5 
 
sage:  parasurf3dplot=[0 for i in range(M)] 
 
sage:  u, v = var ('u, v', domain='real') 
 
sage:  n=var('n', domain='integer') 
 
sage:  def fx(u,v,n): return ((n+1)/M+cos(u/2)cos(v)-sin(u/2)sin(2v))cos(u) 
 
sage: ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*cos(u) 
def fy(u,v,n): return ((n+1)/M+cos(u/2)cos(v)-sin(u/2)sin(2v))sin(u) 
 
sage: ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*sin(u) 
def fz(u,v,n): return sin(u/2)cos(v)+cos(u/2)sin(2*v) 
 
sage: sin(u/2)*cos(v)+cos(u/2)*sin(2*v) 
def rangeu(n): return (u, -pi+(n/M^3)pi,pi-(n/M^3)pi) 
 
sage: -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
def rangev(n): return (v, -pi+(n/M^3)pi,pi-(n/M^3)pi) 
 
sage: -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
parasurf3dplot=[parametric_plot3d((fx(u,v,n), fy(u,v,n), fz(u,v,n)), rangeu(n), rangev(n), \        color=rainbow(M)[n], opacity=0.1,  mesh=True) for n in range(M)] 
 
sage:  P=sum(parasurf3dplot[i] for i in range(M)) 
 
sage: P.show(viewer='jmol') 
P.show(viewer='jmol')

M=5
isosurfplot=[0 for i in range(M)] 
 
 x, y, z = var ('x, y, z', domain='real') 
 
 n=var('n', domain='integer')
 
domain='integer')
isosurfplot=[implicit_plot3d(x^2+y^2+z^2-(n+1)^2, \                                                                                                           (x, -2-2n/(M-1), 2+2n/(M-1)), -2-2*n/(M-1), 2+2*n/(M-1)), (y, -2-2n/(M-1),2+2n/(M-1) -2-2*n/(M-1),2+2*n/(M-1) ), (z, -2-2n/(M-1),2+2n/(M-1)), -2-2*n/(M-1),2+2*n/(M-1)), \                          color=rainbow(M)[n], opacity=0.3) for n in range(M)]
 
range(M)]
P=sum(isosurfplot[i] for i in range(M)) # NOT sum([...])
 
P.show(viewer='jmol') 
P.show(viewer='jmol')

I believe that the following (tested) code provides a simple, yet complete answer to my question:

First, for functional surfaces (z=f(x,y)), the following code provides an alternative to the solution invoking the class ParametrizedSurface3D: here we only use parametric_plot3d:

M=5
funcsurfplot=[0 for i in range(M)]
x, y = var ('x, y', domain='real') 
n=var('n', domain='integer')
funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), \                                                                                         (x, -2-n, 2+n), (y, -2-n, 2+n),  \    color=rainbow(M)[n], opacity=0.3,  mesh=True) for n in range(M)]
P=sum(funcsurfplot[i] for i in range(M))
P.show(viewer='jmol')

Note that in this case we have full control over the domains of each of the functions in the family of surfaces, and that these domains can also depend on the integer parameter of the family. The line "funcsurfplot=[0 for i in range(M)]" might be considered redundant, but I prefer it for creation of array-like list in Python with custom size, because, unlike the use of 'append', this approach allows the easy construction of multidimensional arrays (two dimensional arrays would be needed for visualizing the families of coordinate surfaces in the example I gave in the formulation of the question) -- see, e.g.,

http://www.i-programmer.info/programming/python/3942-arrays-in-python.html

Now all the other cases of the question can be resolved in the same simple uniform way, as follows:

Second, meshes of 3D curves (curvilinear coordinates in 3D and on surfaces):

M=5
curve3dplot=[0 for i in range(M)] # creation of array of size M+1 via initialization
t=var ('t', domain='real')
n=var('n', domain='integer')
curve3dplot=[parametric_plot3d((exp(-(n+t)^2), cos(n*t*pi), 1/(1+n*t^2)),  (t, -1, 1), \                color=rainbow(M)[n], thickness=4) for n in range(M)]
P=sum(curve3dplot[i] for i in range(M))
P.show(viewer='jmol')

Third, families of parametric surfaces (the surfaces chosen here are Klein bottles of the default type in SAGE - see

http://doc.sagemath.org/html/en/reference/riemannian_geometry/sage/geometry/riemannian_manifolds/surface3d_generators.html

and, in order to introduce dependence on n also in the parametric domains, I have slightly 'unzipped' them):

M=5 
parasurf3dplot=[0 for i in range(M)] 
u, v = var ('u, v', domain='real') 
n=var('n', domain='integer') 
def fx(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*cos(u) 
def fy(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*sin(u) 
def fz(u,v,n): return sin(u/2)*cos(v)+cos(u/2)*sin(2*v) 
def rangeu(n): return (u, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
def rangev(n): return (v, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
parasurf3dplot=[parametric_plot3d((fx(u,v,n), fy(u,v,n), fz(u,v,n)), rangeu(n), rangev(n), \        color=rainbow(M)[n], opacity=0.1,  mesh=True) for n in range(M)] 
P=sum(parasurf3dplot[i] for i in range(M)) 
P.show(viewer='jmol')

Fourth, families of isosurfaces:

M=5
isosurfplot=[0 for i in range(M)] 
x, y, z = var ('x, y, z', domain='real') 
n=var('n', domain='integer')
isosurfplot=[implicit_plot3d(x^2+y^2+z^2-(n+1)^2, \                                                                                                           (x, -2-2*n/(M-1), (x,-2-2*n/(M-1), 2+2*n/(M-1)), (y, -2-2*n/(M-1),2+2*n/(M-1) (y,-2-2*n/(M-1),2+2*n/(M-1) ), (z, -2-2*n/(M-1),2+2*n/(M-1)), (z,-2-2*n/(M-1),2+2*n/(M-1)), \                          color=rainbow(M)[n], opacity=0.3) for n in range(M)]
P=sum(isosurfplot[i] for i in range(M)) 
P.show(viewer='jmol')

Enjoy the results and feel very welcome to improve and optimize this code. I suggest that in the next version of the SAGE tutorial and/or reference manual there appears a subtopic dedicated to this, as a 3D-upgrade of the topic about contour plots. If, for this to happen, assistance is required from me, I will provide such assistance.

shao-linux :)

My previous text:

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and the use of jmol, as discussed in the link given above, would be solved in a very satisfactory way.

I believe that the following (tested) code provides a simple, yet complete answer to my question:

First, for functional surfaces (z=f(x,y)), the following code provides an alternative to the solution invoking the class ParametrizedSurface3D: here we only use parametric_plot3d:

M=5
funcsurfplot=[0 for i in range(M)]
x, y = var ('x, y', domain='real') 
n=var('n', domain='integer')
funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), \                                                                                         (x, -2-n, 2+n), (y, -2-n, 2+n), \                                                                                                                                     color=rainbow(M)[n], opacity=0.3,  mesh=True) for n in range(M)]
P=sum(funcsurfplot[i] for i in range(M))
P.show(viewer='jmol')

Note that in this case we have full control over the domains of each of the functions in the family of surfaces, and that these domains can also depend on the integer parameter of the family. The line "funcsurfplot=[0 for i in range(M)]" might be considered redundant, but I prefer it for creation of array-like list in Python with custom size, because, unlike the use of 'append', this approach allows the easy construction of multidimensional arrays (two dimensional arrays would be needed for visualizing the families of coordinate surfaces in the example I gave in the formulation of the question) -- see, e.g.,

http://www.i-programmer.info/programming/python/3942-arrays-in-python.html

Now all the other cases of the question can be resolved in the same simple uniform way, as follows:

Second, meshes of 3D curves (curvilinear coordinates in 3D and on surfaces):

M=5
curve3dplot=[0 for i in range(M)] # creation of array of size M+1 via initialization
t=var ('t', domain='real')
n=var('n', domain='integer')
curve3dplot=[parametric_plot3d((exp(-(n+t)^2), cos(n*t*pi), 1/(1+n*t^2)),  (t, -1, 1), \                                 color=rainbow(M)[n], thickness=4) for n in range(M)]
P=sum(curve3dplot[i] for i in range(M))
P.show(viewer='jmol')

Third, families of parametric surfaces (the surfaces chosen here are Klein bottles of the default type in SAGE - see

http://doc.sagemath.org/html/en/reference/riemannian_geometry/sage/geometry/riemannian_manifolds/surface3d_generators.html

and, in order to introduce dependence on n also in the parametric domains, I have slightly 'unzipped' them):

M=5 
parasurf3dplot=[0 for i in range(M)] 
u, v = var ('u, v', domain='real') 
n=var('n', domain='integer') 
def fx(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*cos(u) 
def fy(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*sin(u) 
def fz(u,v,n): return sin(u/2)*cos(v)+cos(u/2)*sin(2*v) 
def rangeu(n): return (u, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
def rangev(n): return (v, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
parasurf3dplot=[parametric_plot3d((fx(u,v,n), fy(u,v,n), fz(u,v,n)), rangeu(n), rangev(n), \        color=rainbow(M)[n], opacity=0.1,  mesh=True) for n in range(M)] 
P=sum(parasurf3dplot[i] for i in range(M)) 
P.show(viewer='jmol')

Fourth, families of isosurfaces:

M=5
isosurfplot=[0 for i in range(M)] 
x, y, z = var ('x, y, z', domain='real') 
n=var('n', domain='integer')
isosurfplot=[implicit_plot3d(x^2+y^2+z^2-(n+1)^2, \                                                                                                           (x,-2-2*n/(M-1), 2+2*n/(M-1)), (y,-2-2*n/(M-1),2+2*n/(M-1) ), (z,-2-2*n/(M-1),2+2*n/(M-1)), (x,-2-2*n/(M-1),2+2*n/(M-1)),(y,-2-2*n/(M-1),2+2*n/(M-1) ),(z,-2-2*n/(M-1),2+2*n/(M-1)), \                           color=rainbow(M)[n], opacity=0.3) for n in range(M)]
P=sum(isosurfplot[i] for i in range(M)) 
P.show(viewer='jmol')

Enjoy the results and feel very welcome to improve and optimize this code. I suggest that in the next version of the SAGE tutorial and/or reference manual there appears a subtopic dedicated to this, as a 3D-upgrade of the topic about contour plots. If, for this to happen, assistance is required from me, I will provide such assistance.

shao-linux :)

My previous text:

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and the use of jmol, as discussed in the link given above, would be solved in a very satisfactory way.

I believe that the following (tested) code provides a simple, yet complete answer to my question:

First, for functional surfaces (z=f(x,y)), the following code provides an alternative to the solution invoking the class ParametrizedSurface3D: here we only use parametric_plot3d:

M=5
funcsurfplot=[0 for i in range(M)]
x, y = var ('x, y', domain='real') 
n=var('n', domain='integer')
funcsurfplot=[parametric_plot3d((x, y, exp(-(x^2+y^2)/(n+1))), \                                                                                         (x, -2-n, 2+n), (y, -2-n, 2+n), \                                                                                                                                     color=rainbow(M)[n], opacity=0.3,  mesh=True) for n in range(M)]
P=sum(funcsurfplot[i] for i in range(M))
P.show(viewer='jmol')

Note that in this case we have full control over the domains of each of the functions in the family of surfaces, and that these domains can also depend on the integer parameter of the family. The line "funcsurfplot=[0 for i in range(M)]" might be considered redundant, but I prefer it for creation of array-like list in Python with custom size, because, unlike the use of 'append', this approach allows the easy construction of multidimensional arrays (two dimensional arrays would be needed for visualizing the families of coordinate surfaces in the example I gave in the formulation of the question) -- see, e.g.,

http://www.i-programmer.info/programming/python/3942-arrays-in-python.html

Now all the other cases of the question can be resolved in the same simple uniform way, as follows:

Second, meshes of 3D curves (curvilinear coordinates in 3D and on surfaces):

M=5
curve3dplot=[0 for i in range(M)] # creation of array of size M+1 via initialization
t=var ('t', domain='real')
n=var('n', domain='integer')
curve3dplot=[parametric_plot3d((exp(-(n+t)^2), cos(n*t*pi), 1/(1+n*t^2)),  (t, -1, 1), \                                 color=rainbow(M)[n], thickness=4) for n in range(M)]
P=sum(curve3dplot[i] for i in range(M))
P.show(viewer='jmol')

Third, families of parametric surfaces (the surfaces chosen here are Klein bottles of the default type in SAGE - see

http://doc.sagemath.org/html/en/reference/riemannian_geometry/sage/geometry/riemannian_manifolds/surface3d_generators.html

and, in order to introduce dependence on n also in the parametric domains, I have slightly 'unzipped' them):

M=5 
parasurf3dplot=[0 for i in range(M)] 
u, v = var ('u, v', domain='real') 
n=var('n', domain='integer') 
def fx(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*cos(u) 
def fy(u,v,n): return ((n+1)/M+cos(u/2)*cos(v)-sin(u/2)*sin(2*v))*sin(u) 
def fz(u,v,n): return sin(u/2)*cos(v)+cos(u/2)*sin(2*v) 
def rangeu(n): return (u, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
def rangev(n): return (v, -pi+(n/M^3)*pi,pi-(n/M^3)*pi) 
parasurf3dplot=[parametric_plot3d((fx(u,v,n), fy(u,v,n), fz(u,v,n)), rangeu(n), rangev(n), \        color=rainbow(M)[n], opacity=0.1,  mesh=True) for n in range(M)] 
P=sum(parasurf3dplot[i] for i in range(M)) 
P.show(viewer='jmol')

Fourth, families of isosurfaces:

M=5
isosurfplot=[0 for i in range(M)] 
x, y, z = var ('x, y, z', domain='real') 
n=var('n', domain='integer')
isosurfplot=[implicit_plot3d(x^2+y^2+z^2-(n+1)^2, \                                                                                                           (x,-2-2*n/(M-1),2+2*n/(M-1)),(y,-2-2*n/(M-1),2+2*n/(M-1) ),(z,-2-2*n/(M-1),2+2*n/(M-1)), (x,-2-2*n/(M-1),2+2*n/(M-1)),(y,-2-2*n/(M-1),2+2*n/(M-1)),(z,-2-2*n/(M-1),2+2*n/(M-1)), \                            color=rainbow(M)[n], opacity=0.3) for n in range(M)]
P=sum(isosurfplot[i] for i in range(M)) 
P.show(viewer='jmol')

Enjoy the results and feel very welcome to improve and optimize this code. I suggest that in the next version of the SAGE tutorial and/or reference manual there appears a subtopic dedicated to this, as a 3D-upgrade of the topic about contour plots. If, for this to happen, assistance is required from me, I will provide such assistance.

shao-linux :)

My previous text:

Thank you for your comment – I apologize for the long and untidy presentation. Here are the latex formulae which got messed up in my text:

$$F(x,y,z)=exp(-x^2)((cos(y\pi))^2)(1/(1+z^2)), \quad -1\leq x,y,z\leq 1.$$

In the example, x_i, y_j, z_k and c_l are strictly increasing knot-vectors between 0 and 1.

Here is a piece of SAGE code for a simple case which is tested and produces the necessary graphical output via jmol:

sage: x,y=var('x,y');

sage: def f0(x,y,n): return x, y, exp(-(x2 + y2));

sage: def f1(x,y): return x, y, exp(-(x2 + y2)/2);

sage: def f0(x,y): return x, y, exp(-(x2 + y2)/1);

sage: def f2(x,y): return x, y, exp(-(x2 + y2)/3);

sage: def f3(x,y): return x, y, exp(-(x2 + y2)/4);

sage: def f4(x,y): return x, y, exp(-(x2 + y2)/5);

sage: P0 = ParametricSurface(f0, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[0], opacity=0.5);

sage: P1 = ParametricSurface(f1, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[1], opacity=0.5);

sage: P2 = ParametricSurface(f2, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[2], opacity=0.5);

sage: P3 = ParametricSurface(f3, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[3], opacity=0.5);

sage: P4 = ParametricSurface(f4, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(5)[4], opacity=0.5);

sage: P=P0+P1+P2+P3+P4;

sage: P.show(viewer='jmol')

Launched jmol viewer for Graphics3d Object

This works, but is very inconvenient in this form, because the increase of the number of curves from 5 to, say, N where N can be large, e.g., N=100, would require the definition of 95 more functions f5,...,f99 and 95 more plots P5,...P99, after which we would need to write explicitly

$$P=P1+\dots+P99$$

Evidently, one should prefer something of the type (this new piece of 'code' is just to clarify my request: clearly it is not going to work in SAGE 'as is'):

N=100;

x,y,n=var('x,y,n');

def f(x,y,n): return x,y,n, exp(-(x^2 + y^2)/(n+1));

for n in range(N): def Pplot(n): return n, ParametricSurface(f, (srange(-2, 2, 0.1), srange(-2, 2, 0.1)), color=rainbow(N)[n], opacity=0.5);

P=sum(Pplot(n), n in range(N));

P.show(viewer='jmol')

One might argue that $exp(-(x^2 + y^2)/(n+1))$ is also SAGE-predefined function, so let us assume that the function f(x,y,n) is NOT a SAGE-predefined function (e.g., defined by an elliptic integral which is being computed numerically for x and y).

Here ParametricSurface is involved; for the cases described in the example, this idea must work also with parametric_plot3d (for curves and for surfaces in 3D) and with implicit_plot3d, preferably both variants A and B as described in my text.

See also the discussion at http://sourceforge.net/p/jmol/mailman/message/23809652/

If you provide answers for the proposed example, the problem about meshes and the use of jmol, as discussed in the link given above, would be solved in a very satisfactory way.