"""Original Gaussian displacement ellipsoids; Python3.14, NumPy2.3.5, Matplotlib3.10.7.
Run from the mirrored media directory; output is written to the working directory.
"""
import math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

def chi3_cdf(c):
 return math.erf(math.sqrt(c/2))-math.sqrt(2*c/math.pi)*math.exp(-c/2)
lo,hi=0.0,20.0
for _ in range(80):
 mid=(lo+hi)/2
 if chi3_cdf(mid)<0.5:lo=mid
 else:hi=mid
c=(lo+hi)/2
assert abs(chi3_cdf(c)-0.5)<1e-12
u=np.linspace(0,2*np.pi,55);v=np.linspace(0,np.pi,27)
unit=np.stack([np.outer(np.cos(u),np.sin(v)),np.outer(np.sin(u),np.sin(v)),np.outer(np.ones_like(u),np.cos(v))])
fig=plt.figure(figsize=(8,3.8),dpi=100,facecolor='white')
for k,(ev,title,color)in enumerate([(np.array([.01,.04,.09]),'Anisotropic: U = diag(0.01, 0.04, 0.09)','#2166ac'),(np.full(3,.14/3),'Isotropic: U = (0.14/3) I','#b35806')],1):
 ax=fig.add_subplot(1,2,k,projection='3d');r=np.sqrt(c*ev);xyz=unit*r[:,None,None]
 ax.plot_wireframe(*xyz,rstride=3,cstride=3,color=color,linewidth=.65)
 for j in range(3):
  p=np.zeros((3,2));p[j]=[-r[j],r[j]];ax.plot(*p,color=color,linewidth=1.6)
 ax.set(xlim=(-.5,.5),ylim=(-.5,.5),zlim=(-.5,.5),xlabel='x (Å)',ylabel='y (Å)',zlabel='z (Å)',title=title)
 ax.set_xticks([-.4,0,.4]);ax.set_yticks([-.4,0,.4]);ax.set_zticks([-.4,0,.4]);ax.set_box_aspect((1,1,1));ax.view_init(elev=22,azim=38);ax.tick_params(labelsize=8);ax.title.set_size(9)
fig.suptitle('Same Ueq = 0.04667 Ų; different directional spread',fontsize=12,y=.98)
fig.text(.5,.025,'Both surfaces enclose 50% probability in a three-dimensional Gaussian model.',ha='center',fontsize=9)
fig.subplots_adjust(left=.015,right=.94,top=.83,bottom=.12,wspace=.09)
fig.savefig('atomic-displacement-ellipsoid-comparison.png',dpi=100,facecolor='white')
print(f'c={c:.12f}; semiaxes_anisotropic={np.sqrt(c*np.array([.01,.04,.09]))}; dimensions=800x380')