"""Original unilateral gear-contact model. Python 3.14, NumPy 2.3.5, Matplotlib 3.10.7.
Run from output cwd: python ../gear-engagement-tangent-stiffness.py
"""
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

d=np.linspace(0,7,701)
f1=80*d
f2=120*np.maximum(d-3,0)
f=f1+f2
assert np.isclose(f[300],240)
assert np.isclose(80*4.8+120*(4.8-3),600)
plt.rcParams.update({"font.size":10,"axes.spines.top":False,"axes.spines.right":False})
fig,(ax,bx)=plt.subplots(1,2,figsize=(8.5,3.6),dpi=100,layout="constrained")
fig.patch.set_facecolor("white")
fig.suptitle("Gear engagement: force is continuous, tangent stiffness jumps",fontsize=12)
ax.plot(d,f,color="#243747",lw=2,label="Total force")
ax.plot(d,f1,color="#3777a6",lw=1.6,label="Pair 1")
ax.plot(d,f2,color="#27865c",lw=1.6,label="Pair 2")
ax.axvline(3,color="#8b6a34",ls="--",lw=1)
ax.scatter([3],[240],color="#243747",s=25,zorder=4)
ax.annotate("Pair 2 starts here",xy=(3,240),xytext=(.4,740),arrowprops={"arrowstyle":"->","color":"#8b6a34"},color="#765628",fontsize=9)
ax.set(xlabel="Normal approach (micrometres)",ylabel="Normal force (N)",xlim=(0,7),ylim=(0,1100))
ax.legend(loc="lower right",fontsize=8,frameon=False)
bx.plot([0,3],[80,80],color="#245c8b",lw=2)
bx.plot([3,7],[200,200],color="#245c8b",lw=2,label="One-sided tangent stiffness")
bx.scatter([3,3],[80,200],facecolors="white",edgecolors="#245c8b",s=42,zorder=4)
bx.scatter([3],[140],marker="x",color="#b04a32",s=48,zorder=4)
bx.annotate("Centered difference: 140",xy=(3,140),xytext=(.3,160),fontsize=9,color="#963d2b",arrowprops={"arrowstyle":"->","color":"#b04a32"})
bx.axvline(3,color="#8b6a34",ls="--",lw=1)
bx.text(4.4,90,"No unique tangent\nat the threshold",ha="center",fontsize=9)
bx.set(xlabel="Normal approach (micrometres)",ylabel="Tangent stiffness (N/micrometre)",xlim=(0,7),ylim=(40,230))
for a in (ax,bx):a.grid(alpha=.17)
fig.savefig(Path.cwd()/"gear-engagement-tangent-stiffness.png",facecolor="white",dpi=100)
plt.close(fig)