Files
2025-09-29 00:52:08 +02:00

136 lines
2.3 KiB
Plaintext
Executable File

--Catenary solve
--Translated by Neil Gregory from C++ code supplied by Bernie Freidin
struct CatenaryCurve
(
points = #(),
slack = 0.3,
maxIter = 75,
fn asinh_rad x =
(
log(x + sqrt(1.0 + x*x))
),
fn solveCatenary x0 dx m b =
(
--print "solveCatenary"
while cosh (radToDeg(x0 + dx)) < m * (x0 + dx) + b do
(
dx *= 2.0
)
x1 = x0 + dx
for iter = 0 to maxIter do
(
x2 = 0.5 * (x1 + x0)
y2 = cosh (radToDeg x2)
if y2 < m*x2 + b then
(
x0 = x2
)
else
(
x1 = x2
)
)
0.5*(x1 + x0)
),
fn solveSlack &x0 &x1 xmid m b =
(
--print "solveSlack"
x0 = solveCatenary xmid -1.0 m b
x1 = solveCatenary xmid 1.0 m b
dx = x1 - x0
dy = cosh (radToDeg x1) - cosh (radToDeg x0)
arc = sinh (radToDeg x1) - sinh (radToDeg x0) --arc length along cosh(x) between [x0,x1]
dist = sqrt (dx*dx + dy*dy)
arc / dist - 1.0
),
fn solveHangingWire &x0 &x1 m slack =
(
--print "solvehangingWire"
xmid = asinh_rad m
b0 = sqrt (1.0 + m*m) - m*xmid --m*x + b0 is tangent to cosh(x) at x = xmid
db = 1.0
if slack <= 0.0 then
(
x0 = xmid
x1 = xmid
return b0
)
while (solveSlack &x0 &x1 xmid m b0 + db) < slack do
(
db *= 2.0
)
b1 = b0 + db
for iter = 0 to maxIter do
(
b2 = 0.5*(b1 + b0)
if (solveSlack &x0 &x1 xmid m b2) < slack then
(
b0 = b2
)
else
(
b1 = b2
)
)
--print "."
0.5*(b1 + b0)
),
fn simulateHangingWire p0 p1 numSegs =
(
--clear points array
points = #()
v = p1 - p0
m = v.z / length [v.x, v.y, 0]
x0 = undefined
x1 = undefined
solveHangingWire &x0 &x1 m slack
y0 = cosh (radToDeg x0)
y1 = cosh (radToDeg x1)
points[1] = p0
points[numSegs] = p1
for i = 1 to numSegs do
(
t = i as Float / numSegs as Float
x = x0 + (x1 - x0)*t
--x = asinh_rad (sinh (radToDeg x0) + t * sinh (radToDeg x1) - sinh (radToDeg x0))
y = cosh (radToDeg x)
pxy = p0 + (p1 - p0)*t
--pxy = p0 + (p1 - p0) * ((x - x0)/(x1 - x0))
pz = p0 + (p1 - p0) * ((y - y0)/(y1 - y0))
points[i] = [pxy.x, pxy.y, pz.z]
)
)
)
global catenaryCurve = CatenaryCurve()