First of all full disclaimer, I am new to python and my question may be silly, but as far as my research has taught me it seems there is no supported way in matplotlib to do what I want.
I have to plot trajectories of the Three Body Problem and I am trying to do it with python to be less dependent of Matlab (which is what they teach at my University). In these plots appear both Earth and Moon, and it is important for me that they look spherical. The next code gives a basic example of the kind of thing that I have to plot:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D, proj3d # Needed to 3Dplots
plt.close('all') # Close plots of previous scripts
#Data
Er = 0.0166 # Earth radius
Mr = 0.0045 # Moon radius
mu = 0.01215 # Mass parameter
#
def plotCC(ax,rad,pos): # Plot celestial bodies
u = np.linspace(0, 2 * np.pi, 53)
v = np.linspace(0, np.pi, 53)
x = pos + rad * np.outer(np.cos(u), np.sin(v))
y = rad * np.outer(np.sin(u), np.sin(v))
z = rad * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z)
# Dummy Halo orbit
u = np.linspace(0, 2 * np.pi, 30)
xh = (1-mu)*np.ones(len(u))
yh = 20*Mr*np.sin(u)
zh = 0.12 + 30*Mr*np.cos(u)
# Dummy trajectory
xt = np.linspace(1-mu,-mu,20)
yt = (xt-(1-mu))*(xt+mu)*-1.2
zt = (xt-(1-mu))*(xt+mu)*-0.8
Using that I can plot all I need with:
# Plots
fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d',proj_type = 'ortho')
ax.plot(xh,yh,zh, lw=1.2) # Plot Halo orbit
ax.plot(xt,yt,zt, lw=1.2) # Plot trajectory
plotCC(ax,Mr,1-mu) # Plot Moon
plotCC(ax,Er,-mu) # Plot Earth
Which gives as a result this:
It is far of what I need, as none of the celestial bodies are represented as spheres. I thought that the aspect-ratio could be configured to match that of the represented data, but none of the supported methods that I found seem to do just that. This methods are: ax.axis('equal'), ax.set_aspect('equal'), ax.axis('scaled'), ax.set_aspect('equal', adjustable='box'), ax.set_aspect('equal', adjustable='datalim'). None of these do the job.
The first closest thing that I found is this answer, which basically creates a box so that all axis measure the same. This does the trick, but in my case one of the dimensions is much bigger than the others, and this results in a lot of wasted space:
Plot with the solution that creates a box
Finally, the best thing that I could find is the solution in here, which requires to modify the file axes.py so that get_proj() accepts an aspect ratio that you give to it with:
ay = (ax.get_ylim3d()[1]-ax.get_ylim3d()[0])/(ax.get_xlim3d()[1]-ax.get_xlim3d()[0])
az = (ax.get_zlim3d()[1]-ax.get_zlim3d()[0])/(ax.get_xlim3d()[1]-ax.get_xlim3d()[0])
ax.pbaspect = list(1.8*np.array([1,ay,az]))
The result is this. If you look closely the plot isn't centered, and depending in the orientation not all the data is inside the drawing area. For now I will use this last method, but after all the research that I have done I still can't believe that something like this isn't supported as it should by a plotting tool (which is fantastic in everything else as far as I know), and thus my question: Am I missing something obvious?
Thanks in advance :)
