"""Earth-gravity note: spherical coordinates, Legendre series blow-up, spherical-harmonic patterns, (n,m) grid.""" import numpy as np import matplotlib.pyplot as plt from numpy.polynomial import legendre as L from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine def Pnm(n, m, x): return (1 - x**2)**(m/2) * L.Legendre.basis(n).deriv(m)(x) if m else L.legval(x, [0]*n + [1]) # ------------------------------------------------------------------ spherical coordinates (3-D) fig = plt.figure(figsize=(6.2, 5.6)) ax = fig.add_subplot(111, projection="3d") u, v = np.meshgrid(np.linspace(0, 2*np.pi, 60), np.linspace(-np.pi/2, np.pi/2, 30)) for lat in np.radians([-60, -30, 0, 30, 60]): t = np.linspace(0, 2*np.pi, 200) ax.plot(np.cos(lat)*np.cos(t), np.cos(lat)*np.sin(t), np.sin(lat)*np.ones_like(t), color=INK, lw=.6, alpha=.45) for lon in np.radians(np.arange(0, 360, 30)): t = np.linspace(-np.pi/2, np.pi/2, 100) ax.plot(np.cos(t)*np.cos(lon), np.cos(t)*np.sin(lon), np.sin(t), color=INK, lw=.6, alpha=.45) lam, ph, r = np.radians(50), np.radians(35), 1.7 Pt = r*np.array([np.cos(ph)*np.cos(lam), np.cos(ph)*np.sin(lam), np.sin(ph)]) foot = r*np.array([np.cos(ph)*np.cos(lam), np.cos(ph)*np.sin(lam), 0]) ax.plot([0, Pt[0]], [0, Pt[1]], [0, Pt[2]], color=BLUE, lw=2.4) ax.plot([0, foot[0]], [0, foot[1]], [0, 0], color=INK, lw=1, ls="--") ax.plot([foot[0], Pt[0]], [foot[1], Pt[1]], [0, Pt[2]], color=INK, lw=1, ls="--") t = np.linspace(0, lam, 40); ax.plot(.75*np.cos(t), .75*np.sin(t), 0*t, color=GREEN, lw=2.4) t = np.linspace(0, ph, 40); ax.plot(.9*np.cos(t)*np.cos(lam), .9*np.cos(t)*np.sin(lam), .9*np.sin(t), color=AMBER, lw=2.4) ax.plot([0, 2.0], [0, 0], [0, 0], color=INK, lw=1); ax.text(1.6, -.4, -.9, "prime meridian ($\\lambda=0$)", fontsize=9) ax.plot([0, 0], [0, 0], [0, 2.0], color=INK, lw=1.2); ax.text(0, 0, 2.05, "spin axis $z$", fontsize=9) ax.text(*(Pt*1.04), r"$P=(r,\varphi,\lambda)$", color=BLUE, fontsize=11) ax.text(.95, .45, 0, r"$\lambda$", color=GREEN, fontsize=14) ax.text(1.05*np.cos(lam)*np.cos(ph/2)*.9, 1.05*np.sin(lam)*np.cos(ph/2)*.9, .95*np.sin(ph/2)+.04, r"$\varphi$", color=AMBER, fontsize=14) ax.text(*(Pt*.5+[.0, .0, .12]), r"$r$", color=BLUE, fontsize=14) ax.set_box_aspect((1, 1, 1), zoom=1.5); ax.set_xlim(-1.4, 1.9); ax.set_ylim(-1.4, 1.9); ax.set_zlim(-1.2, 2.1) ax.set_axis_off(); ax.view_init(22, 28) fig.savefig("img/fig-grav-sphcoords.png", dpi=160, bbox_inches="tight"); plt.close(fig) # ------------------------------------------------------------------ Legendre ODE series: only integer n is regular at the poles def series(nu, x, N=4000): a = np.zeros(N); a[0] = 1.0 for j in range(0, N-2, 2): a[j+2] = (j*(j+1) - nu*(nu+1)) / ((j+1)*(j+2)) * a[j] return np.polyval(a[::-1], x) fig, (a1, a2) = plt.subplots(1, 2, figsize=(11.5, 4.2)) x = np.linspace(0, .9995, 800) for nu, c, lab in ((1.5, BLUE, r"$\nu=1.5$"), (2.0, INK, r"$\nu=2$ (terminates: $1-3x^2=-2P_2$)"), (2.5, MAROON, r"$\nu=2.5$")): y = np.array([series(nu, xi) for xi in x]) if nu != 2.0 else 1 - 3*x**2 a1.plot(x, y, color=c, lw=2.4 if nu == 2 else 1.8, label=lab) a1.set_ylim(-6, 6); a1.axhline(0, color=INK, lw=.5) a1.set_xlabel(r"$x=\sin\varphi$ (right edge $=$ north pole)"); a1.set_ylabel(r"even solution $X(x)$, $X(0)=1$") a1.set_title("Only integer $\\nu$ stays finite at the pole", fontsize=11.5); a1.legend(frameon=False, fontsize=9.5, loc="lower left"); despine(a1) j = np.arange(2, 60, 2) for nu, c in ((1.5, BLUE), (2.5, MAROON)): a = np.zeros(200); a[0] = 1 for k in range(0, 198, 2): a[k+2] = (k*(k+1) - nu*(nu+1))/((k+1)*(k+2))*a[k] a2.plot(np.arange(0, 60, 2), np.abs(a[0:60:2]), "o-", ms=3.5, color=c, label=rf"$\nu={nu}$") a2.plot(np.arange(2, 60, 2), 1/np.arange(2, 60, 2)*1.0, color=INK, ls="--", label=r"$\propto 1/j$") a2.set_yscale("log"); a2.set_xlabel("power $j$ in $X=\\sum a_j x^j$"); a2.set_ylabel(r"$|a_j|$") a2.set_title(r"Coefficients only fall like $1/j$, so $\sum a_j x^j\sim -\ln(1-x)$", fontsize=11.5); a2.legend(frameon=False); despine(a2) fig.tight_layout(); fig.savefig("img/fig-grav-legendre-series.png", dpi=160, bbox_inches="tight"); plt.close(fig) # ------------------------------------------------------------------ spherical-harmonic patterns lon = np.radians(np.linspace(-180, 180, 361)); lat = np.radians(np.linspace(-90, 90, 181)) LO, LA = np.meshgrid(lon, lat) cases = [(2, 0, "zonal"), (4, 0, "zonal"), (2, 2, "sectoral"), (4, 4, "sectoral"), (3, 1, "tesseral"), (5, 3, "tesseral")] fig, axs = plt.subplots(2, 3, figsize=(11.5, 6.2)) for ax, (n, m, kind) in zip(axs.flat, cases): F = Pnm(n, m, np.sin(LA)) * np.cos(m*LO) F = F / np.abs(F).max() ax.pcolormesh(np.degrees(LO), np.degrees(LA), F, cmap="RdBu_r", vmin=-1, vmax=1, shading="auto") ax.contour(np.degrees(LO), np.degrees(LA), F, levels=[0], colors=INK, linewidths=1.3) ax.set_title(rf"$P_{{{n}{m}}}(\sin\varphi)\cos {m if m>1 else ''}\lambda$ — {kind}" if m else rf"$P_{{{n}}}(\sin\varphi)$ — {kind}", fontsize=11) ax.set_xticks([-180, -90, 0, 90, 180]); ax.set_yticks([-90, -45, 0, 45, 90]) ax.text(.02, .04, f"{n-m} nodal parallels,\n{2*m} nodal meridians", transform=ax.transAxes, fontsize=8.8, bbox=dict(fc=CREAM, ec="none", alpha=.8)) for ax in axs[1]: ax.set_xlabel("longitude $\\lambda$ (deg)") for ax in axs[:, 0]: ax.set_ylabel("latitude $\\varphi$ (deg)") fig.tight_layout(); fig.savefig("img/fig-grav-sphharm.png", dpi=160, bbox_inches="tight"); plt.close(fig) # ------------------------------------------------------------------ (n, m) grid fig, ax = plt.subplots(figsize=(4.8, 4.4)) for n in range(0, 8): for m in range(0, n+1): kind = "zonal" if m == 0 else ("sectoral" if m == n else "tesseral") c = {"zonal": BLUE, "sectoral": MAROON, "tesseral": GREEN}[kind] ax.plot(m, n, "o", color=c, ms=11) for k, c, lab in (("zonal", BLUE, "zonal ($m=0$)"), ("sectoral", MAROON, "sectoral ($m=n$)"), ("tesseral", GREEN, "tesseral ($0