Search code examples
pythonmathmatplotlibcoordinate-systems

Ellipsoid creation in Python


I have ran into a problem relating to the drawing of the Ellipsoid.

The ellipsoid that I am drawing to draw is the following:

x**2/16 + y**2/16 + z**2/16 = 1.

So I saw a lot of references relating to calculating and plotting of an Ellipse void and in multiple questions a cartesian to spherical or vice versa calculation was mentioned.

Ran into a website that had a calculator for it, but I had no idea on how to successfully perform this calculation. Also I am not sure as to what the linspaces should be set to. Have seen the ones that I have there as defaults, but as I got no previous experience with these libraries, I really don't know what to expect from it.

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=plt.figaspect(1))  # Square figure
ax = fig.add_subplot(111, projection='3d')

multip = (1, 1, 1) 
# Radii corresponding to the coefficients:
rx, ry, rz = 1/np.sqrt(multip)

# Spherical Angles
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# Cartesian coordinates

#Lots of uncertainty.
#x = 
#y = 
#z = 

# Plot:
ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')

# Axis modifications
max_radius = max(rx, ry, rz)
for axis in 'xyz':
    getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))

plt.show()

Solution

  • Your ellipsoid is not just an ellipsoid, it's a sphere.

    Notice that if you use the substitution formulas written below for x, y and z, you'll get an identity. It is in general easier to plot such a surface of revolution in a different coordinate system (spherical in this case), rather than attempting to solve an implicit equation (which in most plotting programs ends up jagged, unless you take some countermeasures).

    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    import numpy as np
    
    phi = np.linspace(0,2*np.pi, 256).reshape(256, 1) # the angle of the projection in the xy-plane
    theta = np.linspace(0, np.pi, 256).reshape(-1, 256) # the angle from the polar axis, ie the polar angle
    radius = 4
    
    # Transformation formulae for a spherical coordinate system.
    x = radius*np.sin(theta)*np.cos(phi)
    y = radius*np.sin(theta)*np.sin(phi)
    z = radius*np.cos(theta)
    
    fig = plt.figure(figsize=plt.figaspect(1))  # Square figure
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, color='b')
    

    enter image description here