import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Line3DCollection
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.size'] = 15
# Set up 1x2 subplots
# Create a figure and a standard 2D axes array
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
# Remove the second (right-side) axis
axes[1].remove()
# Add a new 3D axis in the same position
ax3d = fig.add_subplot(1, 2, 2, projection='3d')
# --- 3D Ellipsoid Plot Setup ---
R = 0.25 # Half-distance between foci
f1 = np.array([0, 0, R])
f2 = np.array([0, 0, -R])
# Set viewing angle
view_azimuth = 45
view_elevation = 30
ax3d.view_init(elev=view_elevation, azim=view_azimuth)
# Calculate view direction vector
view_az_rad = np.radians(view_azimuth)
view_el_rad = np.radians(view_elevation)
view_dir = np.array([
np.cos(view_el_rad) * np.cos(view_az_rad),
np.cos(view_el_rad) * np.sin(view_az_rad),
np.sin(view_el_rad)
])
def ellipsoidal_to_cartesian(mu, nu, phi, R):
"""Convert prolate spheroidal coordinates to Cartesian"""
x = R * np.sqrt((mu**2 - 1) * (1 - nu**2)) * np.cos(phi)
y = R * np.sqrt((mu**2 - 1) * (1 - nu**2)) * np.sin(phi)
z = R * mu * nu
return x, y, z
# --- Generate and plot ellipsoid surface ---
mu_surface = 2.0
nu_vals = np.linspace(-1, 1, 100)
phi_vals = np.linspace(0, 2 * np.pi, 100)
Nu, Phi = np.meshgrid(nu_vals, phi_vals)
X = R * np.sqrt((mu_surface**2 - 1) * (1 - Nu**2)) * np.cos(Phi)
Y = R * np.sqrt((mu_surface**2 - 1) * (1 - Nu**2)) * np.sin(Phi)
Z = R * mu_surface * Nu
ax3d.plot_surface(X, Y, Z, alpha=0.3, color='lightblue', edgecolor='none')
# --- Generate and plot phi lines ---
n_phi_lines = 12
phi_lines_visible = []
phi_lines_hidden = []
for phi in np.linspace(0, 2 * np.pi, n_phi_lines, endpoint=False):
nu = np.linspace(-1, 1, 200)
x, y, z = ellipsoidal_to_cartesian(mu_surface, nu, phi, R)
points = np.column_stack([x, y, z])
visibility = np.dot(points, view_dir)
for i in range(len(nu) - 1):
if (visibility[i] + visibility[i+1]) / 2 > 0:
phi_lines_visible.append(np.array([x[i:i+2], y[i:i+2], z[i:i+2]]).T)
else:
phi_lines_hidden.append(np.array([x[i:i+2], y[i:i+2], z[i:i+2]]).T)
if phi_lines_visible:
phi_vis_coll = Line3DCollection(phi_lines_visible, colors='g', linewidth=1.5, linestyle='-', alpha=0.8)
ax3d.add_collection(phi_vis_coll)
if phi_lines_hidden:
phi_hid_coll = Line3DCollection(phi_lines_hidden, colors='g', linewidth=0.8, linestyle='--', alpha=0.3)
ax3d.add_collection(phi_hid_coll)
# --- Generate and plot nu lines ---
n_nu_lines = 9
nu_lines_visible = []
nu_lines_hidden = []
for nu in np.linspace(-0.9, 0.9, n_nu_lines):
phi = np.linspace(0, 2 * np.pi, 200)
x, y, z = ellipsoidal_to_cartesian(mu_surface, nu * np.ones_like(phi), phi, R)
points = np.column_stack([x, y, z])
visibility = np.dot(points, view_dir)
for i in range(len(phi) - 1):
if (visibility[i] + visibility[i+1]) / 2 > 0:
nu_lines_visible.append(np.array([x[i:i+2], y[i:i+2], z[i:i+2]]).T)
else:
nu_lines_hidden.append(np.array([x[i:i+2], y[i:i+2], z[i:i+2]]).T)
if nu_lines_visible:
nu_vis_coll = Line3DCollection(nu_lines_visible, colors='r', linewidth=1.5, linestyle='-', alpha=0.8)
ax3d.add_collection(nu_vis_coll)
if nu_lines_hidden:
nu_hid_coll = Line3DCollection(nu_lines_hidden, colors='r', linewidth=0.8, linestyle='--', alpha=0.3)
ax3d.add_collection(nu_hid_coll)
# --- Finalize plot ---
ax3d.scatter(*f1, color='r', s=100, marker='o', label='Focus 1', zorder=10)
ax3d.scatter(*f2, color='r', s=100, marker='o', label='Focus 2', zorder=10)
ax3d.set_box_aspect(aspect=(1,1,1))
max_range = 0.3
ax3d.set_xlim([-max_range, max_range])
ax3d.set_ylim([-max_range, max_range])
ax3d.set_zlim([-max_range, max_range])
ax3d.set_axis_off()
# You can add a plot to the 2D axes here
# Plot constant mu lines (varying nu)
# --- Generate and plot grid on the 2D subplot (axes) ---
# Project onto the xz-plane, for phi = 0 and phi = 180 degrees
phi_xz_values = [0, np.pi] # 0 and 180 degrees in radians
for phi_xz in phi_xz_values:
# Plot constant mu lines (varying nu)
mu_lines_to_plot = [1., 1.25, 1.5, 1.75, 2.0]
nu_vals = np.linspace(-1, 1, 100)
for mu in mu_lines_to_plot:
x_mu, _, z_mu = ellipsoidal_to_cartesian(mu, nu_vals, phi_xz, R)
axes[0].plot(x_mu, z_mu, 'b-', alpha=0.6)
# Add label for the mu line
label_text = rf'$\xi={{{mu:.2f}}}$'
# Position label near the line end. Adjust spacing as needed.
if (phi_xz == 0): axes[0].text(x_mu[50]+0.025, z_mu[50]+0.015, label_text, fontsize=8, rotation=-90, ha='center', color='blue')
# Plot constant nu lines (varying mu)
nu_lines_to_plot = [-0.9, -0.6, -0.3, 0, 0.3, 0.6, 0.9]
mu_vals = np.linspace(1.0, 2.0, 100) # mu > 1
for nu in nu_lines_to_plot:
x_nu, _, z_nu = ellipsoidal_to_cartesian(mu_vals, nu, phi_xz, R)
axes[0].plot(x_nu, z_nu, 'r-', alpha=0.6)
axes[0].text(-0.225, 0.505, rf'$\eta=0.9$', fontsize=8, rotation=-61, ha='center', va='center', color='red')
axes[0].text(-0.41, 0.34, rf'$\eta=0.6$', fontsize=8, rotation=-35, ha='center', va='center', color='red')
axes[0].text(-0.48, 0.168, rf'$\eta=0.3$', fontsize=8, rotation=-15, ha='center', va='center', color='red')
axes[0].text(-0.5, 0, rf'$\eta=0$', fontsize=8, rotation=0, ha='center', va='center', color='red')
axes[0].text(-0.5, -0.18, rf'$\eta=-0.3$', fontsize=8, rotation=15, ha='center', va='center', color='red')
axes[0].text(-0.42, -0.35, rf'$\eta=-0.6$', fontsize=8, rotation=35, ha='center', va='center', color='red')
axes[0].text(-0.22, -0.525, rf'$\eta=-0.9$', fontsize=8, rotation=61, ha='center', va='center', color='red')
axes[0].scatter([0,0], [R,-R], color='r', s=100, marker='o', label='Focus 1', zorder=10)
axes[0].set_xlabel(r'$x$')
axes[0].set_ylabel(r'$z$')
axes[0].set_box_aspect(1) # Equal aspect ratio
max_range = 0.6
axes[0].set_xlim([-max_range, max_range])
axes[0].set_ylim([-max_range, max_range])
plt.tight_layout()
plt.show()
