# catSolver.py # Andrew Davison, Nov 2025, ad@coe.psu.ac.th import math from scipy.optimize import brentq def getCatA(H=None, d=None, L=None, tol=1e-10, maxIter=100): """ Calculates the catenary parameter 'a' for any valid combination of two of span (H), sag (d), and arc length (L), while avoiding numerical overflow. The function solves the catenary equations: d = a(cosh(H/(2a)) - 1) for (H, d) L = 2a sinh(H/(2a)) for (H, L) depending on which pair of parameters is supplied. * Use Brent’s method, brentq() from SciPy, for solving transcendental equations. Overflow-prevention checks: * Rejects physically impossible (H, L) combinations (L <= H) * Enforces a minimum safe a-value to prevent sinh(H/(2a)) overflow * Protects sinh/cosh evaluations by checking argument size -- note that both the d and L equation use arguments of the form h/2a * In root-finding, returns +inf when evaluation would overflow """ numParams = sum(x is not None for x in (H, d, L)) if numParams != 2: raise ValueError("Two of (H, d, L) must be provided") if H is not None and H <= 0: raise ValueError("H must be positive") if d is not None and d <= 0: raise ValueError("d must be positive") if L is not None and L <= 0: raise ValueError("L must be positive") # Overflow safety constants ARG_MAX = 700.0 # cosh/sinh overflow near ~710 def maxA(Hval): # calculates a such that H/(2a) == ARG_MAX return Hval / (2.0 * ARG_MAX) # Prevent impossible cases. if H is not None and L is not None: if L <= H * 1.0000001: # A catenary with span H has min arc length H raise ValueError( f"L={L} too close to span H={H}; forces a->0 and sinh(H/(2a)) overflow." ) # CASE 1: Closed-form (L,d) if L is not None and d is not None and H is None: a = (0.25 * L*L - d*d) / (2*d) if a <= 0: raise ValueError(f"Non-physical a obtained: {a}") if a < maxA(L): raise ValueError("Computed a too small; would overflow.") return a # CASE 2: Known H and d if H is not None and d is not None: ratio = d / H # Small sag branch (ratio < 0.1) if ratio < 0.1: aApprox = H*H / (8*d) if aApprox < maxA(H): aApprox = maxA(H) # Very small sag (ratio < 0.01) if ratio < 0.01: def smallSag(a): u = H / (2*a) if abs(u) > ARG_MAX: return float('inf') # safe fallback return 2*a * math.sinh(u/2)**2 - d # sinh-squared identity is more reliable here aLower = max(aApprox * 0.5, maxA(H)) aUpper = aApprox * 2.0 if aUpper <= aLower: aUpper = aLower * 1.1 return brentq(smallSag, aLower, aUpper, xtol=tol, maxiter=maxIter) # Moderate small sag: use Newton-Raphson a = aApprox for _ in range(maxIter): u = H / (2*a) if abs(u) > ARG_MAX: a = max(a, maxA(H)) u = H / (2*a) # Use sinh-squared identity for f and derivative f = 2*a * math.sinh(u/2)**2 - d df = 2 * math.sinh(u/2)**2 - (H/(2*a)) * math.sinh(u) if abs(df) < 1e-15: break aNew = a - f/df if aNew < maxA(H): aNew = maxA(H) if abs(aNew - a) < tol*abs(a): return aNew a = aNew return a # Larger sag branch def solveHd(a): # function used by brentq() for solving u = H / (2*a) # d = a(cosh(H/(2a)) - 1) for a if abs(u) > ARG_MAX: return float('inf') return a*(math.cosh(u)-1) - d aInit = H*H / (8*d) aInit = max(aInit, maxA(H)) aLower = max(aInit * 0.1, maxA(H)) aUpper = aInit * 10.0 # Ensure bracket safety for brentq() while solveHd(aLower) * solveHd(aUpper) > 0: aLower = max(aLower * 0.5, maxA(H)) aUpper *= 2.0 if aUpper > 1e12: raise ValueError("Cannot bracket solution safely (overflow risk).") return brentq(solveHd, aLower, aUpper, xtol=tol, maxiter=maxIter) # CASE 3: Known H and L if H is not None and L is not None: def solveHL(a): # function used by brentq() for solving u = H / (2*a) # L = 2a sinh(H/(2a)) for a if abs(u) > ARG_MAX: return float('inf') return 2*a * math.sinh(u) - L # Parabolic-based initial guess aInit = H * L / (L*L - H*H) if L > H else H aInit = max(aInit, maxA(H)) aLower = max(aInit * 0.1, maxA(H)) aUpper = aInit * 10.0 # Ensure bracket safety for brentq() while solveHL(aLower) * solveHL(aUpper) > 0: aLower = max(aLower * 0.5, maxA(H)) aUpper *= 2.0 if aUpper > 1e12: raise ValueError("Cannot bracket solution for H,L (overflow risk).") return brentq(solveHL, aLower, aUpper, xtol=tol, maxiter=maxIter) # ---------------------------------------- if __name__ == "__main__": print("Testing catenary parameter computation:\n") # Test 1: L and d (closed form) print("Test 1: L=110m, d=10m") a1 = getCatA(L=110, d=10) print(f" a = {a1:.6f} m\n") # Test 2: H and d with small sag print("Test 2: H=100m, d=0.5m (small sag, d/H = 0.005)") a2 = getCatA(H=100, d=0.5) print(f" a = {a2:.6f} m") print(f" Verification: d = {a2 * (math.cosh(100/(2*a2)) - 1):.6f} m\n") # Test 3: H and d with moderate sag print("Test 3: H=100m, d=10m (moderate sag)") a3 = getCatA(H=100, d=10) print(f" a = {a3:.6f} m") print(f" Verification: d = {a3 * (math.cosh(100/(2*a3)) - 1):.6f} m\n") # Test 4: Very small sag print("Test 4: H=1000m, d=0.1m (very small sag, d/H = 0.0001)") a5 = getCatA(H=1000, d=0.1) print(f" a = {a5:.6f} m") print(f" Approximation: H^2/(8h) = {1000**2/(8*0.1):.6f} m") print(f" Verification: d = {a5 * (math.cosh(1000/(2*a5)) - 1):.6f} m\n") # Test 5: H and L print("Test 5: H=100m, L=110m") a4 = getCatA(H=100, L=110) print(f" a = {a4:.6f} m") print(f" Verification: L = {2*a4*math.sinh(100/(2*a4)):.6f} m") # overflow / extreme cases -------------- # Test 6: L extremely close to H (forces a -> 0) print("\nTest 6: L extremely close to span H") try: H = 4 L = 4.0000000001 print("a =", getCatA(H=H, L=L)) except Exception as e: print("Caught:", e) # Test 7: L < H (physically impossible) print("\nTest 7: L < H (physically impossible)") try: print(getCatA(H=10, L=9.9)) except Exception as e: print("Caught:", e) # Test 8: Extremely small sag (d/H -> 0) print("\nTest 8: Extremely small sag") try: print(getCatA(H=10000, d=0.01)) # produces huge value except Exception as e: print("Caught:", e) # Test 9: Very large L (huge arc length) print("\nTest 9: Very large L relative to span") try: print(getCatA(H=10, L=10000)) # produces small value except Exception as e: print("Caught:", e) # Test 10: Very small H, huge L print("\nTest 10: Tiny H, huge L") try: print(getCatA(H=0.1, L=100)) # produces v.small value except Exception as e: print("Caught:", e) # Test 11: Sag ratio d/H -> 0 (numerical boundary) print("\nTest 11: d/H approaching zero") try: print(getCatA(H=5000, d=0.0001)) # produces huge value except Exception as e: print("Caught:", e) # Test 12: Use sinh/cosh argument near overflow boundary print("\nTest 12: Overflow boundary (u approx 700)") try: H = 1000 a_crit = H / (2 * 700.1) # produces u \approx 700.1 (unsafe) L = 2 * a_crit * math.sinh(H/(2*a_crit)) print(getCatA(H=H, L=L)) except Exception as e: print("Caught:", e)