"""Our analysis for https://factchecked.science/galactic-bending-waves/ : arithmetic checks of
Griv, E., and Yahalom, A. (2026), Generation of spiral bending waves in a rotating galactic stellar disk of finite
thickness, Astrophysics and Space Science 371, 118, https://doi.org/10.1007/s10509-026-04651-8 (CC BY 4.0),
using only the numbers printed in the paper (the caption of its Fig. 4, its Sects. 4-7) and in the sources cited
on the page, with standard unit conversions. Standard library only; reads no file.
Usage: python3 -I scale_checks.py
"""
import math
from datetime import date

KM_PER_KPC = 3.085677581e16
SEC_PER_YR = 365.25 * 86400

# --- Fig. 4 caption values (p. 12 of the paper) ---
r0 = 8.2          # kpc
Omega = 28.29     # km/s/kpc
kappa = 40.01     # km/s/kpc
c_r = 35.0        # km/s
c_z = 20.0        # km/s
print("1. Caption consistency")
print(f"   r0*Omega = {r0*Omega:.1f} km/s (caption: 232)")
print(f"   sqrt(2)*Omega = {math.sqrt(2)*Omega:.2f} km/s/kpc (caption kappa: 40.01)")

# --- time scales ---
Omega_s = Omega / KM_PER_KPC              # 1/s
Omega_yr = Omega_s * SEC_PER_YR           # 1/yr
inv_Omega_Myr = 1 / Omega_yr / 1e6
P_rot_Myr = 2 * math.pi * inv_Omega_Myr
print("\n2. Time scales at the caption's Omega")
print(f"   Omega = {Omega_yr:.3e} per yr (paper: 'Omega ~ 10^-8 yr', i.e. per yr)")
print(f"   1/Omega = {inv_Omega_Myr:.1f} Myr; rotation period 2*pi/Omega = {P_rot_Myr:.0f} Myr")
print(f"   10/Omega = {10*inv_Omega_Myr/1000:.2f} Gyr (paper: t_grow ~ 10/Omega ~ 1 Gyr, taking Omega ~ 1e-8 per yr)")
print(f"   10-20 rotation periods = {10*P_rot_Myr/1000:.1f}-{20*P_rot_Myr/1000:.1f} Gyr")
print(f"   e-folds in 10-20 rotation periods at Im(omega) = 0.1*Omega: {0.1*2*math.pi*10:.1f}-{0.1*2*math.pi*20:.1f}")
for g in (0.02, 0.01):
    print(f"   one e-fold at Im(omega) = {g}*Omega takes {inv_Omega_Myr/g/1000:.2f} Gyr")

# --- wavenumber <-> wavelength (x = k^2 c_r^2 / kappa^2, the axis of Fig. 4) ---
L = c_r / kappa   # kpc  (c_r/kappa)
print(f"\n3. Scales behind Fig. 4's x axis: c_r/kappa = {L:.3f} kpc; 2*pi*c_r/kappa = {2*math.pi*L:.2f} kpc")
print("   x       k (1/kpc)   1/k (kpc)   2*pi/k (kpc)")
for x in (0.05, 0.085, 0.1, 0.19, 0.25, 0.40, 0.5, 1.0):
    k = math.sqrt(x) / L
    print(f"   {x:<6}  {k:8.3f}   {1/k:8.2f}    {2*math.pi/k:8.1f}")
print("   Observed full wavelengths -> x and k*h_z:")
for lam in (2.0, 2.5, 3.0):
    k = 2 * math.pi / lam
    x = (k * L) ** 2
    print(f"   2*pi/k = {lam} kpc: k = {k:.2f}/kpc, x = {x:.2f}, k*h_z = {k*0.2:.2f}-{k*0.3:.2f} (h_z = 0.2-0.3 kpc)")
for inv_k in (2.0, 3.0):
    k = 1 / inv_k
    print(f"   1/k = {inv_k} kpc: x = {(k*L)**2:.3f}, full wavelength 2*pi/k = {2*math.pi*inv_k:.1f} kpc")


# --- Re(omega*) from eq. (48), with G rho0 h_z^2 = c_z^2/(4 pi) (h_z^2 ~ c_z^2/mu^2, mu^2 ~ 4 pi G rho0; p. 11) ---
def I1(x, terms=30):
    return sum((x / 2) ** (2 * j + 1) / (math.factorial(j) * math.factorial(j + 1)) for j in range(terms))


coef = (32 / math.pi) * (c_z ** 2 / (4 * math.pi)) / c_r ** 2   # (32/pi) G rho0 h_z^2 / c_r^2
print(f"\n4. Eq. (48): Re(omega*) = kappa*sqrt(1 - {coef:.4f} e^-x I1(x))")
for x in (0.0, 0.1, 0.19, 0.5, 1.0):
    re = kappa * math.sqrt(1 - coef * math.exp(-x) * I1(x))
    print(f"   x = {x:<4}: Re(omega*)/Omega = {re/Omega:.3f}; (Re(omega*) - kappa)/Omega = {(re-kappa)/Omega:+.4f}")
print("   (the paper's text, p. 12, gives Re(omega*) as about Omega; Fig. 4's solid curves sit at 1.37-1.41)")

# --- drift coefficient, eq. (24), flat rotation curve (Omega ~ 1/r) at r0 ---
dOmega_dr = -Omega / r0   # km/s/kpc per kpc
D_cr2 = (c_r ** 2 / (r0 * kappa ** 2)) * dOmega_dr   # km/s/kpc  (eq. 24, second form, v_perp = c_r)
print("\n5. Eq. (24) for a flat rotation curve at r0 = 8.2 kpc")
print(f"   D c_r^2 = {D_cr2:.3f} km/s/kpc = {D_cr2/Omega:+.4f} Omega; 2 D c_r^2 = {2*D_cr2/Omega:+.4f} Omega (paper p. 11: '~ -0.1 Omega')")
print(f"   ratio of the paper's ~0.1 Omega to ours: {0.1/abs(2*D_cr2/Omega):.1f}")
print(f"   |L_Omega| = Omega/|dOmega/dr| = {Omega/abs(dOmega_dr):.1f} kpc (paper p. 8: '|L_Omega| ~ 3 kpc'); r_epic = c_r/kappa = {L:.2f} kpc (paper: '<~ 1 kpc')")
print(f"   with r_epic = 1 kpc and |L_Omega| = 3 kpc: |D v^2| = {1/(r0*3):.3f} Omega, 2|D v^2| = {2/(r0*3):.3f} Omega")
for x in (0.5, 1.0):
    re = kappa * math.sqrt(1 - coef * math.exp(-x) * I1(x))
    print(f"   m_crit (eq. 52) at x = {x}: {(re-kappa)/(2*D_cr2):.2f} with eq.-(24) D; {(re-kappa)/(-0.1*Omega):.2f} with 2Dc_r^2 = -0.1 Omega")

# --- observed ripple sizes in the sources the paper cites for '250-350 pc' and '2-3 kpc' (p. 2) ---
print("\n6. Ripple sizes in the cited sources (abstracts; Konietzka et al. and Schoenrich & Dehnen in the authors' arXiv versions)")
rows = [
    ("Schoenrich & Dehnen 2018", "W wave ~1 km/s (0.8 km/s in fits; 1 km/s ~ 35 pc)", "~2.5 kpc (2-2.5 kpc in R_g)"),
    ("Alves et al. 2020", "max ~160 pc", "average period ~2 kpc"),
    ("Narayan et al. 2020 (5 other galaxies)", "70-300 pc", "1-5 kpc"),
    ("Konietzka et al. 2024", "max ~220 pc", "mean wavelength ~2 kpc"),
    ("Poggio et al. 2025", "~150-200 pc (corrugation on the warp)", "length >= 10 kpc"),
    ("Sun et al. 2026 (not cited by the paper)", "~100-200 pc", "radial 3.9-7.9 kpc"),
]
for r in rows:
    print("   " + " | ".join(r))
print("   paper (p. 2): 'maximum out-of-plane amplitude of 250-350 pc and an average in-plane period of 2-3 kpc'")

# --- years since each report, at the episode date ---
ep = date(2026, 10, 27)
print("\n7. Years before 27 October 2026")
for label, d in [
    ("Burke 1957 / Kerr 1957 (warp; AJ 62, May 1957)", date(1957, 5, 1)),
    ("Varsavsky & Quiroga 1970 (IAU Symp. 38, 'a possible corrugation')", date(1970, 1, 1)),
    ("Quiroga 1974 (Ap&SS 27, April 1974, 'The corrugation of the galactic layer')", date(1974, 4, 1)),
    ("Widrow et al. 2012 (ApJL, online 25 Apr 2012; vertical waves)", date(2012, 4, 25)),
    ("Xu et al. 2015 (ApJ, online 11 Mar 2015; rings and radial waves)", date(2015, 3, 11)),
    ("Antoja et al. 2018 (Nature, online 19 Sep 2018; phase spiral)", date(2018, 9, 19)),
    ("Alves et al. 2020 (Nature, online 7 Jan 2020; Radcliffe Wave)", date(2020, 1, 7)),
    ("Konietzka et al. 2024 (Nature, online 20 Feb 2024; it oscillates)", date(2024, 2, 20)),
    ("Poggio et al. 2025 (A&A, online 14 Jul 2025; the great wave)", date(2025, 7, 14)),
    ("Sun et al. 2026 (Nat. Astron., online 7 Aug 2026; CO corrugations)", date(2026, 8, 7)),
]:
    print(f"   {label}: {(ep - d).days/365.25:.1f} years")