[SOLVED] Hanging rope (catenary)


#1

I know this topic has been discussed here before and to my belief it never got to a true mathematical approach of a hanging rope.

What I am looking for is a semi realistic simulation of a hanging rope between two points (of different height) and a given rope length. I cannot use physics or MassFX or anything like that, because I will be using it is a live max viewing session and I need things to behave quick and effortless.

Now this rope can be devoid of any physical parameters. I don’t need it to stretch or have elasticity as it is called, but I would like to know the actual lowest point and the way the rope or cable should hang. So basically, I would like to generate a curve with N points, and I would like to calculate the height of said points.

Behaviour of a rope can be described using a Catenary, a mathematical description of a rope. I have giving the formulas for that shape a lot of thought this past week, but the math is beyond me at the moment.

Now I have found two example codes that I have both tried to convert to maxscript with varying result:
http://members.chello.nl/j.beentjes3/Ruud/catfiles/catenary.pdf#page=7&zoom=100,0,0
https://nl.mathworks.com/matlabcentral/fileexchange/38550-catenary-hanging-rope-between-two-points
(see functions and then the catenary function)

I have made a dummy that has a scale_script controller that controls the position of the spline points. These points are calculated using the logic, but somehow my logic isn’t sound.

While the below image looks good it fails to behave as it should when the right point is lower than the left point, and the lenght also isn’t always close to the actual target length.

below is my Maxscript snippet of my logic from the second solution that is giving me the best result so far. I have racked my brain on this for some days now, so if someone could help me with the final stretch that would be very much appreciated.

My question would be, can someone help me with proper logic for a catenary, or perhaps point me in the right direction as I am also not sure what I am doing wrong here.

fn g s d r_length h= 2.0*sinh(s*d/2.0)/s - sqrt(r_length^2.0-h^2.0)

fn dg s d = 2.0*cosh(s*d/2.0)*d/(2.0*s) - 2*sinh(s*d/2.0)/(s^2.0)

fn catenary a b r_length N sagInit = (
	-- given two points a=[ax ay] and b=[bx by] in the vertical plane,
	-- rope length r_length, and the number of intermediate points N,
	-- outputs the coordinates X and Y of the hanging rope from a to b
	-- the optional input sagInit initializes the sag parameter for the
	-- root-finding procedure.
	maxIter	= 100		-- maximum number of iterations
	minGrad	= 1e-10		-- minimum norm of gradient
	minVal	= 1e-8		-- minimum norm of sag function
	stepDec	= 0.5			-- factor for decreasing stepsize
	minStep	= 1e-9		-- minimum step size
	minHoriz	= 1e-3		-- minumum horizontal distance
	if saginit == undefined then (
		sag = 1
	)
	else (
		sag = sagInit
	)
	
	pointA = a
	pointB = b  	
	
	d = distance pointB [pointA.x,pointA.y,pointB.z]
	h = distance pointA [pointA.x,pointA.y,pointB.z]
	
	PointArray = #()
	
	Xdist = pointB.x - pointA.x
	Ydist = pointB.y - pointA.z
	Zdist = pointB.z - pointA.z
	
	for i=0 to N do 
	(
		P = [PointA.x + ((Xdist/N)*i),pointA.y+((Ydist/N)*i),pointA.z]
		append PointArray P
	)
	
	-- rope is stretched: straight line
	if r_length <= sqrt(d^2+h^2) then (
		for i=0 to N do
		(
			P = PointArray[i+1]
			P.z = PointA.z + ((Zdist/N)*i)
			PointArray[i+1] = P
		)
	)
	
	-- find rope sag
	else 
	(
		for i= 1 to maxIter do (
			val = g sag d r_length h
			grad = dg sag h
			if (abs val) < minVal or (abs grad) < minGrad then (
				break
			)
			
			search	= -(g sag d r_length h )/(dg sag h)
			alpha		= 1
			sag_new  = sag + alpha*search
			
			while sag_new < 0 or (abs(g sag_new d r_length h)) > (abs val) do 
			(
				alpha		= stepDec*alpha
				if alpha < minStep then 
				(
					break
				)
				sag_new = sag + alpha*search			
			)
			sag = sag_new
		)
		
		Xstep = d/N
		
		-- get location of rope minimum and vertical bias
		x_left	= 0.5*(log((r_length+h)/(r_length-h))/sag-d)
		
		for i=0 to N do (
			X = Xstep*i

			x_min		= -x_left
			bias		= pointA.z - cosh(x_left*sag)/sag
			Z			= cosh((X-x_min)*sag)/sag + bias
			
			P = PointArray[i+1]
			P.z = pointA.z + Z
			PointArray[i+1] = P
		)
	)
	PointArray
)


a = catenary [0,0,0] [20,0,5] 28 10 1

new_spline = line()
addNewSpline new_spline
for i=1 to a.count do (
	addKnot new_spline 1 #smooth #curve a[i]
)
updateshape new_spline

#2

I have found another very interesting example Javascript that does basically what I need, and is even visualized but again I fail to convert this to proper maxscript… :frowning: I must be overlooking something.

The example can be found here: https://codepen.io/sketchpunk/pen/ajGrqZ

As I compare my code to the code of the Js I cannot see what is different, but somehow it just doesn’t do the same thing. :frowning: Do any of you see what I am missing?

My current new code:

fn asinh x = 
(
	log(x+(x*x+1)^0.5)
)

fn catenary a x = 
(
	a*(cosh(x/a))
)

fn getA posA posB ropeLength = 
(
	--Solving A comes from : http://rhin.crai.archi.fr/rld/plugin_details.php?id=990
	zDelta = posB.z - posA.z
	vecLen = distance posA posB
	
	if (zDelta > ropeLength) or (vecLen > ropeLength) then
	(
		exit
	)
	--Swap verts, low end needs to be on the left side
	if (zDelta < 0.0) then 
	(
		tmp 	= posA
		posA		= posB
		posB		= tmp
		zDelta	*= -1.0
	)
	
	max_tries = 100
	diff = 0.0
	a = 100.0
	aTmp = 0.0
	zRopeDelta = 0.5 * sqrt((ropeLength*ropeLength) - (zDelta*zDelta))
	vecLenHalf = 0.5 * vecLen
	
	for i=0 to max_tries do 
	(
		aTmp	= vecLenHalf / (asinh (zRopeDelta/a))
		diff = abs( (aTmp - a) / a )
		a = aTmp
		if diff < 0.0001 then 
		(			
			exit
		)
	)
	--return a
	a
)

fn getPoints posA posB ropeLength segCnt = 
(
	dist = distance posA posB
	A = getA posA posB ropeLength
	format "a = %\n" a
	distHalf 	= dist * 0.5
	segInc = dist / segCnt
	offset = catenary A (-distHalf)
	format "offset = %\n" offset
	pnt = [0,0,0]
	xpos = undefined
	c = undefined
	
	pntArray = #(posA)
	for i=1 to segCnt do 
	(
		pnt = posA+(((posB - posA)/segCnt)*i)
		xpos	= i * segInc - distHalf
		--format "xpos = %\n" xpos 
		c = catenary A xpos
		--format "c = %\n" c
		--format "offset = %\n" offset
		pnt.z -= offset - c
		append pntArray pnt
	)
	--return the resulting pntArray
	pntArray
)

posA = [0,0,0]
posB = [200,0,50]
cableLength = 280
segCnt = 20

pointArray = getPoints posA posB cableLength segCnt
new_spline = line()
addNewSpline new_spline
for i=1 to segCnt do (
	addKnot new_spline 1 #corner #line pointArray[i]
)
updateshape new_spline

#3

EDIT: turns out this was discussed back in 2015 and it is all on me :frowning: : [Max 2010] The Hyperbolic Trig functions are wrong

My goodness!!! Turns out the default function cosh x is not working properly in max 2017. Not sure about other versions.

I thought I would test all the default functions once more, and voila there is was:

cosh 0.5 returns 1.00004 instead of the actual 1.12763.

I have made my own altCosh function that actually works and now I have a proper hanging rope curve:

fn altCosh x = ((e^x + e^ -x)/2)

I tried searching if anyone else ran into this issue, but it seems not, so let this be a warning to all of you who are hoping to use this function in the future :wink:

Look at those curves:


#4

I’m glad you found the solution yourself. I wanted to explain the problem this morning, but my California morning has just come …

As you can see from the MXS help, all trigonometric functions are in degrees. Including hyperbolic. So the corresponding:


fn sh x = (exp(x) - exp(-x))/2

sh x  == sinh (radtodeg x)

Another major problem with all mxs hyperbolic functions is that they are in floats, so it is very easy to reach their floating point limit during computation. Therefore, wherever possible, try to work with normalized values.


#5


#6

Well … I was playing around with ‘catenaty’ about five years ago … looking for a fast realtime solution (a vertex shader).

But … finally found and stayed with a parabolic solution (sort of ‘parabolic’).
I know that a parabolic curve looks “sharper” than a “catenary”, but it solves 10 times faster!


#7

That’s is looking real nice and thanks for the replies! That does clear things up some more. I am having rounding errors during calculation, so the actual target length is never truly maintained. This is probably due to what you mentioned about the hyperbolic functions being expressed in floats, but at the moment I am not confident how to ‘work in normalized’ values. I was glad I got this far :relaxed:

I have hooked up my catenary calculation to a scale script on one of the dummies and can update the spline’s vertex positions that way. I must say, I am not disappointed about the speed of the update at all. I’ve done some quick tests with 20+ of them and things are still looking real sunny.

I’ll make a nice GIF when the workday is over and I have access to screen capture software :wink:

In the meantime, I am very interested in your parabolic approach. Are you able to retain a nice line length? (up until the points are too far spaced)


#8

But as you can see there is a difference in the target distance and the actual measured distance (curvedistance $) of the cables. Now I did find 2 other script files (C# and python) that I will also convert and test for accuracy, but I’m guessing they will give similar results.

How is the precision on your end @denisT?


#9

The value of the precision is well known. The smaller the deflection value (the difference in chain length and distance between supports), the higher the precision. Those from 0.3% to 5.0%.
The parabolic curve is “sharper” than the Catenary. But you can make some “correction” for the version of the Bezier curve …


#10

several people asked me to post my solution of catenary. here is it (pure MXS version) …


--fn _sinh x = (exp(x) - exp(-x))/2
--fn _cosh x = (exp(x) + exp(-x))/2

fn _cosh x = cosh (radtodeg x)
fn _sinh x = sinh (radtodeg x)

fn makePointTab count start end = 
(
	step = (end - start)/(count-1)
	for k=0 to count-1 collect (start + step*k)
)

fn compPoint3Tab tab xtab ytab ztab = 
(
	for p in tab do 
	(
		append xtab p.x
		append ytab p.y
		append ztab p.z
	)
)

fn combPoint3Tab xtab ytab ztab scale:1 = 
(
	n = amin #(xtab.count, ytab.count, ztab.count)
	for k=1 to n collect ([xtab[k], ytab[k], ztab[k]] * scale)
)
fn cutPointTab tab = 
(
	tab.count -= 1
	tab
)
	
fn catenary pt0 pt1 extent count &sag:1 max_iterations:100 =
(
	pt0 = pt0 / extent
	pt1 = pt1 / extent
	
	len = 1.0
	
	fn interp s len d h = 
	( 
		2 * _sinh(s*d/2)/s - sqrt(len*len - h*h)
	)
	fn deriv s len d h = 
	(
		2 * _cosh(s*d/2)*d/(2*s) - 2 * _sinh(s*d/2)/(s*s)
	)

	max_iterations			-- maximum number of iterations
	min_grad	 	= 1e-10	-- minimum norm of gradient
	min_val	 		= 1e-8	-- minimum norm of sag function
	step_size 		= 0.5	-- factor for decreasing step_size
	min_step	 	= 1e-9	-- minimum step size
		
	if pt0.x > pt1.x do swap pt1 pt0
	
	d = pt1.x - pt0.x
	h = pt1.z - pt0.z
	
	pts = makePointTab count pt0 pt1
	
	ptx = #()
	pty = #()
	ptz = #()
	
	compPoint3Tab pts ptx pty ptz
	
	if len > sqrt(d*d + h*h) do -- rope is not stretched: straight line
	(
		-- find rope sag
		end = off
		for i=1 to max_iterations while not end do
		(
			val	= interp sag len d h
			grad = deriv sag len d h
			
			end = abs(val) < min_val or abs(grad) < min_grad
			
			if not end do
			(
				search	= -(interp sag len d h)/(deriv sag len d h)
				
				alpha = 1
				_sag = sag + alpha*search
				
				out = false
				while (_sag < 0 or abs(interp _sag len d h) > abs(val)) and not out do
				(
					alpha *= step_size
					if not (out = alpha < min_step) do
					(
						_sag = sag + alpha*search
					)
				)
				
				sag = _sag
			)
		)
		
		-- get location of rope minimum and vertical bias
		x_left = 0.5 * (log((len + h)/(len - h))/sag - d)
		x_min = pt0.x - x_left
		
		_bias = pt0.z - _cosh(x_left * sag)/sag
		
		ptz	= #()
		for x in ptx do append ptz (_cosh((x - x_min)*sag)/sag + _bias)
	)
	
	pts = combPoint3Tab ptx pty ptz scale:extent
)

fn setSplineKnots sp pp type:#corner =
(
	for k=1 to pp.count do 
	(
		if numknots sp < k then addknot sp 1 type #line pp[k] 
		else setKnotPoint sp 1 k pp[k] 
	)
	updateshape sp
)	

fn applyCatenary dd extent sp: count:20 = 
(
	sag = 1 
	pp = catenary dd[1].pos dd[2].pos extent count sag:&sag max_iterations:10
	
	if sp == unsupplied do sp = shapes[1]
	if not iskindof sp splineshape do sp = splineshape vertexTicks:on wirecolor:orange
	if (numsplines sp == 0) do addnewspline sp
	
	setSplineKnots sp pp type:#corner --#auto
)


/*
delete objects
gc()
d0 = dummy boxsize:[10,10,10] pos:[-10,0,40]
d1 = dummy boxsize:[5,5,5] pos:[10,0,40]


(
	dd = for obj in helpers where iskindof obj Dummy collect obj  
	--new_spline.vertexTicks = off

	extent = 100 

	deleteAllChangeHandlers id:#catenary
	when transform dd change id:#catenary do applyCatenary dd extent --(distance dd[1] dd[2] * 3) 
	applyCatenary dd extent
)

global new_spline = shapes[1]
*/

it is open to code optimizations and improvements


#11

I am using extent (extra scale) to minimize the disadvantages of mxs float calculations


#12

Yes, real effective! Thanks again, hopefully I can play with it today a little


#13

I noticed there was a small issue with the code above. The catenary worked splendid on the x plane, but was looking off in the y plane. (moving one of the helpers to the same x value as the other)

I added a small portion of code to the catenary function to account for this. Basically I’m just looking at the caterary in the x plane and then translate the calculated positions to the plane between the two points.
(see top and bottom of the function)

Thanks to you I now have a perfect and quick catenary:

fn catenary pt0 pt1 extent count &sag:1 max_iterations:100 =
(
	--build local coordsys for flat approach
	xMax = distance pt0 [pt1.x,pt1.y,pt0.z]
	zMax = pt1.z - pt0.z
	
	oldpt0 = pt0
	oldpt1 = pt1
	
	
	pt0 = [0,0,0]
	pt1 = [xMax,0,zMax] / extent
	
	len = 1.0
	
	fn interp s len d h = 
	( 
		2 * _sinh(s*d/2)/s - sqrt(len*len - h*h)
	)
	fn deriv s len d h = 
	(
		2 * _cosh(s*d/2)*d/(2*s) - 2 * _sinh(s*d/2)/(s*s)
	)

	max_iterations			-- maximum number of iterations
	min_grad	 	= 1e-10	-- minimum norm of gradient
	min_val	 		= 1e-8	-- minimum norm of sag function
	step_size 		= 0.5	-- factor for decreasing step_size
	min_step	 	= 1e-9	-- minimum step size
		
	if pt0.x > pt1.x do swap pt1 pt0
	
	d = pt1.x - pt0.x
	h = pt1.z - pt0.z
	
	pts = makePointTab count pt0 pt1
	
	ptx = #()
	pty = #()
	ptz = #()
	
	compPoint3Tab pts ptx pty ptz
	
	if len > sqrt(d*d + h*h) do -- rope is not stretched: straight line
	(
		-- find rope sag
		end = off
		for i=1 to max_iterations while not end do
		(
			val	= interp sag len d h
			grad = deriv sag len d h
			
			end = abs(val) < min_val or abs(grad) < min_grad
			
			if not end do
			(
				search	= -(interp sag len d h)/(deriv sag len d h)
				
				alpha = 1
				_sag = sag + alpha*search
				
				out = false
				while (_sag < 0 or abs(interp _sag len d h) > abs(val)) and not out do
				(
					alpha *= step_size
					if not (out = alpha < min_step) do
					(
						_sag = sag + alpha*search
					)
				)
				
				sag = _sag
			)
		)
		
		-- get location of rope minimum and vertical bias
		x_left = 0.5 * (log((len + h)/(len - h))/sag - d)
		x_min = pt0.x - x_left
		
		_bias = pt0.z - _cosh(x_left * sag)/sag
		
		ptz	= #()
		for x in ptx do append ptz (_cosh((x - x_min)*sag)/sag + _bias)
	)
	
	--scale back to original values
	pts = combPoint3Tab ptx pty ptz scale:extent

	--recalculate in actual 3d plane
	for i=1 to pts.count do 
	(
		pt = pts[i]
		horFactor = pt.x / xMax
		
		pt.x = oldPt0.x + ((oldPt1.x - oldpt0.x) * horFactor)
		pt.y = oldPt0.y + ((oldPt1.y - oldpt0.y) * horFactor)
		pt.z = oldPt0.z + pt.z
		
		pts[i] = pt
	)
	pts
)

I also tried and successfully got a parabola to work but the sharpness at the bottom of the curve was bothering me. Now I have the ‘real’ thing at least. I also figured out what I was doing wrong in my own catenary function. I was superclose but the normalizing and properly setting the bias was where I was off. Thanks again man