SD15 -- 3D Plotting using Sage
system:sage

<p>&nbsp;</p>
<h1 style="text-align: center;">3D Plotting using Sage</h1>
<h2 style="text-align: center; color: gray;">Bill Cauchois</h2>
<p>&nbsp;</p>

{{{id=1|

///
}}}

<h2>How It All Began</h2>
<p>We couldn't figure out how to embed fast high-quality 3D graphics in the web-based Sage notebook interface until <strong>William Stein discovers Jmol in 2007.</strong></p>
<p><strong>What is Jmol?</strong> Jmol is "a Java viewer for chemical structures in 3D".</p>
<p>Developed as an <strong>open-source alternative to proprietary chemical viewers</strong> and later with the express intent to allow users to <strong>display molecules embedded in web pages</strong>.</p>
<ul>
<li>Renders arbitrary meshes and isosurfaces.</li>
<li>Written in Java, can be embedded as an applet in a web page.</li>
<li>Fast software rendering engine allows Jmol to work even on old computers.</li>
<li>Free, open-source software!</li>
</ul>
<p>In January 2008, <strong>William Stein and Robert Bradshaw release the first version of 3D plotting in Sage</strong>. Includes IndexFaceSet, which can render an arbitrary set of vertices and faces.</p>

{{{id=3|

///
}}}

<h2>Recent Work</h2>
<ul>
<li>Carl Witty contributes code that implements plotting contours of functions (like $f(x, y, z) = 0$).</li>
<li>Me and Jason Grout are adding a feature to graph only parts of 3D plots (analagous to RegionFunction in Mathematica).</li>
<li>Robert Bradshaw writes doctests for plot3d/base.pyx.</li>
<li>Your name here! 3D plotting has lots of room for improvement. Can you help us add that feature you really want?</li>
</ul>

{{{id=5|

///
}}}

<h2>On To the Examples</h2>

<p>Let's begin with an example from calculus. We calculate and plot the plane tangent to
a point.</p>

{{{id=11|
var('x y')
f(x, y) = x^2 + y^2
fx, fy = f.diff(x), f.diff(y)

# Don't be scared by this syntax! interact is really handy.
@interact
def _(x0=slider(-2, 2, 0.5, 0), y0=slider(-2, 2, 0.5, 0)):
    z0 = f(x0, y0)
    # Now we have the point (x0, y0, z0).
    
    # Plug into the formula for the tangent plane at a point.
    tangent_plane(x, y) = fx(x0, y0)*(x-x0) + fy(x0, y0)*(y-y0) + z0
    
    G = plot3d(f, (x,-2,2), (y,-2,2))
    # Now add the plane and the point to the graph G
    G = G + plot3d(tangent_plane, (x,-2,2), (y,-2,2), opacity=0.5, color='green')
    G = G + point((x0, y0, z0), color='red')
    show(G)
///

<html><!--notruncate--><div padding=6 id='div-interact-11'> <table width=800px height=20px bgcolor='#c5c5c5'
                 cellpadding=15><tr><td bgcolor='#f9f9f9' valign=top align=left><table><tr><td align=right><font color="black">x0&nbsp;</font></td><td><table><tr><td>
        <div id='slider-x0-11' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
        </td><td><font color='black' id='slider-x0-11-lbl'></font></td></tr></table><script>(function(){ var values = ["-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000"]; setTimeout(function() {
    $('#slider-x0-11').slider({
        stepping: 1, min: 0, max: 8, startValue: 4,
        change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-x0-11-lbl').text(values[position]); interact(11, "sage.server.notebook.interact.update(11, \"x0\", 5, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(11)"); },
        slide: function(e,ui) { if(values!=null) $('#slider-x0-11-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-x0-11-lbl').text(values[$('#slider-x0-11').slider('value')]);
    }, 1); })();</script></td></tr>
<tr><td align=right><font color="black">y0&nbsp;</font></td><td><table><tr><td>
        <div id='slider-y0-11' class='ui-slider ui-slider-3' style='margin:0px;'><span class='ui-slider-handle'></span></div>
        </td><td><font color='black' id='slider-y0-11-lbl'></font></td></tr></table><script>(function(){ var values = ["-2.00000000000000","-1.50000000000000","-1.00000000000000","-0.500000000000000","0.000000000000000","0.500000000000000","1.00000000000000","1.50000000000000","2.00000000000000"]; setTimeout(function() {
    $('#slider-y0-11').slider({
        stepping: 1, min: 0, max: 8, startValue: 4,
        change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-y0-11-lbl').text(values[position]); interact(11, "sage.server.notebook.interact.update(11, \"y0\", 6, sage.server.notebook.interact.standard_b64decode(\""+encode64(position)+"\"), globals());sage.server.notebook.interact.recompute(11)"); },
        slide: function(e,ui) { if(values!=null) $('#slider-y0-11-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-y0-11-lbl').text(values[$('#slider-y0-11').slider('value')]);
    }, 1); })();</script></td></tr>
</table><div id='cell-interact-11'><?__SAGE__START>
        <table border=0 bgcolor='#white' width=100% height=100%>
        <tr><td bgcolor=white align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
        <tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
        </table><?__SAGE__END></div></td>
                 </tr></table></div>
                 </html>
}}}

<p>Use <strong>parametric_plot3d</strong> to create parametric plots.</p>

{{{id=9|
var('u v')
cone = parametric_plot3d([v*cos(u), v*sin(u), v], (u, 0, 10), (v, 0, 10))
cone
///
}}}

<p>You can apply transformations to plots.</p>
<ul>
<li>G.rotate(v, theta)</em> -- to rotate about a vector</li>
<li>G.rotateX(theta), G.rotateY(theta), G.rotateZ(theta)</li>
<li>G.translate(x, y, z)</li>
<li>G.scale(x, y, z)</li>
</ul>

{{{id=6|
# Take the sum of a list comprehension whose elements are 3D plots!
sum(cone.rotateY(theta) for theta in srange(0, pi*2, pi/2))
///
}}}

<p>There are a number of 3D objects built into Sage.</p>
<ul>
<li>sphere(center, size)</li>
<li>cube(center, size)</li>
<li><em>Platonic solids:</em> dodecahedron, icosahedron, octahedron, tetrahedron</li>
<li>text3d(string, position)</li>
</ul>

{{{id=7|
shapes = [sphere(size=2, color='green'),
          cube(size=4, color='red'),
          icosahedron(size=3, color='blue')]
          
sum(shape.translate(i * 5, 0, 0) for (i, shape) in enumerate(shapes)).show(aspect_ratio=1)
///
}}}

<p>There is also a Tachyon backend, which renders the plot using ray tracing.</p>
<p>Other backends include X3D and Wavefront OBJ (which might be relevant if you are trying to export your plot to use in another 3D program).</p>

{{{id=15|
var('x y')
g(x, y) = sin(pi*(x^2+y^2)) / 2
plot3d(g, (x,-1,1), (y,-1,1)).show(viewer='tachyon')
///
}}}

<p>Graph contours of functions using implicit_plot3d (new in Sage 4.0).</p>

{{{id=8|
var('x y z')
implicit_plot3d(x^2 + y^2 + z^2, (-2, 2), (-2, 2), (-2, 2), contour=4)
///
}}}

<p>You can optionally specify a <strong>region</strong>, in which case the graph will be restricted to a certain set of points defined by a predicate.</p>

{{{id=10|
def R(x, y, z):
    return x <= 0.2 or y >= 0.2 or z <= 0.2
var('x y z')
implicit_plot3d((x^2 + y^2 + z^2), (-2, 2), (-2, 2), (-2, 2),
                plot_points=60, contour=[1,3,5], region=R).show(viewer='tachyon')
///
}}}

<p>This functionality will be included in the other plotting functions soon.</p>

{{{id=14|
var('u v')
parametric_plot3d([u,v,u^2+v^2], (-2, 2), (-2, 2), region=lambda u,v: u^2+v^2>1)
///
}}}

<h2>And Now For Some Cool Plots...</h2>

{{{id=18|
var('u v')
# Boy's surface is an immersion of the real projective plane in 3-space found by Werner Boy in 1901.
boy_x = 2/3* (cos(u)* cos(2*v) + sqrt(2)* sin(u)* cos(v))* cos(u) / (sqrt(2) - sin(2*u)* sin(3*v))
boy_y = 2/3* (cos(u)* sin(2*v) - sqrt(2)* sin(u)* sin(v))* cos(u) / (sqrt(2) - sin(2*u)* sin(3*v))
boy_z = sqrt(2)* cos(u)* cos(u) / (sqrt(2) - sin(2*u)* sin(3*v))
show(parametric_plot3d([boy_x, boy_y, boy_z], (u, -2*pi, 2*pi), (v, 0, pi), plot_points = [90,90], color='orange'))
///
}}}

{{{id=20|
var('x y z')
implicit_plot3d(81*(x^3+y^3+z^3)-189*(x^2*y+x^2*z+y^2*x+y^2*z+z^2*x+z^2*y) +54*x*y*z+126*(x*y+x*z+y*z)-9*(x^2+y^2+z^2)-9*(x+y+z)+1, (x, -1, 1), (y, -1, 1), (z, -1, 1))
///
}}}

{{{id=21|
var('x y z')
implicit_plot3d(cos(x) * sin(y) + cos(y) * sin(z) + cos(z) * sin(x), (x, -4, 4), (y, -4, 4), (z, -4, 4), plot_points=30)
///
}}}

{{{id=26|
# Yoda! -- over 50,000 triangles.
from scipy import io
X = io.loadmat(DATA + 'yodapose.mat')
from sage.plot.plot3d.index_face_set import IndexFaceSet
V = X['V']; F3=X['F3']-1; F4=X['F4']-1
Y = IndexFaceSet(F3,V,color='green') + IndexFaceSet(F4,V,color='green')
Y = Y.rotateX(-1)
Y.show(aspect_ratio=[1,1,1], frame=False, figsize=4)
html('"Use the source, Luke..."')
///
}}}

{{{id=27|

///
}}}