--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()