# Rope¶

Now we’ll try to simulate a rope. Here you’ll see the advantage of PyMbs in
that it is fully scriptable using Python. If you want to try it first, or look
at the complete source code, see
`Rope.py`

.

We start with importing `time`

(you’ll see why in a moment) and PyMbs, and
set up our initial frame of reference as before:

```
import time
from PyMbs.Input import *
world = MbsSystem([0, 0, -1])
```

To make our script a bit more flexible, we add the length of the rope, the rotational damping, the mass density, the radius and the number of segments used to simulate it as variables. Additionally, we choose which method is to be used to derive the equations of motion:

```
L = 3.0 # Total Length
d = 0.2 # Rotational Damping
rho = 500.0 # 500 kg/m^3
r = 0.1 # radius
segments = 15
# Choose between "Explicit", "OrderN" and "Recursive"
method = "Recursive"
```

Now we initialize the lists that will hold our bodies and joints. Also, we calculate the length, mass and inertia of each rope segment:

```
bodies = [None]*segments
joints = [None]*segments
# Parameters
l = L/segments
m = rho*pi*r**2*l
I = 1/12*m*l**2*diag([1, 1, 1])
```

The interesting part is creating all the bodies and joints. We loop over the
number of segments, create a body and an additional frame at the end for each
segment. We then connect the first segment to the world and the other segments
to each other. A simple line is used for visualisation. In the end, a joint
load is added to each joint, with its value being calculated from multiplying
the angular velocity (`s[1]`

) with the rotational damping `d`

:

```
for i in xrange(0, segments):
# Create Body and Frame
bodies[i] = world.addBody(m, cg=[l/2, 0, 0], inertia=I)
bodies[i].addFrame(name='end', p=[l, 0, 0])
# Create Joints
if (i==0):
joints[i] = world.addJoint(world, bodies[i], 'Ry')
else:
joints[i] = world.addJoint(bodies[i-1].end, bodies[i], 'Ry')
world.addVisualisation.Line(bodies[i], l)
# Add Damping
s = world.addSensor.Joint(symbol='state%i_'%i, joint=joints[i])
D = world.addExpression(symbol_str='T%i_'%i, exp=-d*s[1])
world.addLoad.Joint(joint=joints[i], symbol=D)
```

## Method of Equation Generation¶

Now we get to the issue of which method to choose for generating the equations
of motion. The choice is between *Explicit*, *Recursive* and *OrderN*. The
explicit generation scheme is comparatively slow and should not be used for
larger systems. On my system, using eight segments as a benchmark, the
*Explicit* scheme took 12.9 seconds to generate the equations, *OrderN* took
1.06 seconds and the *Recursive* one 0.05 seconds.

You can test these for yourself using the `time`

module:

```
t = time.clock()
if method == "OrderN":
world.genEquations.OrderN()
elif method == "Recursive":
world.genEquations.Recursive()
elif method == "Explicit":
world.genEquations.Explicit()
print "Time needed for generating equations: %5.2f s" % (time.clock() - t)
```

Not only is the time to generate the equations different, the complexity
of the equations differs as well. Again for eight segments, the python code
generated that is used for the simulation has 26.6 million characters (!) when
using *Explicit*, but only 172 thousand when using *OrderN* and 34 thousand
for *Recursive*. So in this case you should use the *Recursive* scheme; for
other models *OrderN* might be beneficial.

Now we can take a look at our rope swinging:

```
world.show('Rope')
```

This is what it should look like:

If you hit simulate with more than say six segments, the simulation will be
quite slow, because the time integration has to evaluate a python function
to obtain the accelerations for each time step. However, you can use the
*Compile Model in Fortran* or *Compile Model in C* buttons to automatically
export the equations to Fortran or C, compile them and reimport them for
faster integration. This way, even 20 segments appear smooth on my machine.

In the end, the complete source looks like this:

```
# -*- coding: utf-8 -*-
'''
This file is part of PyMbs.
PyMbs is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation, either version 3 of
the License, or (at your option) any later version.
PyMbs is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with PyMbs.
If not, see <http://www.gnu.org/licenses/>.
Copyright 2011, 2012 Carsten Knoll, Christian Schubert,
Jens Frenkel, Sebastian Voigt
'''
# Warning: The source code of the examples is quoted in the documentation. If
# you change this file, you'll have to change the corresponding file in the
# documentation (see doc/examples).
import time
from PyMbs.Input import *
world = MbsSystem([0, 0, -1])
L = 3.0 # Total Length
d = 0.2 # Rotational Damping
rho = 500.0 # 500 kg/m^3
r = 0.1 # radius
segments = 15
# Choose between "Explicit", "OrderN" and "Recursive"
method = "Recursive"
# Initialize Lists
bodies = [None]*segments
joints = [None]*segments
# Parameters
l = L/segments
m = rho*pi*r**2*l
I = 1/12*m*l**2*diag([1, 1, 1])
# Create Bodies and connect them
for i in xrange(0, segments):
# Create Body and Frame
bodies[i] = world.addBody(m, cg=[l/2, 0, 0], inertia=I)
bodies[i].addFrame(name='end', p=[l, 0, 0])
# Create Joints
if (i==0):
joints[i] = world.addJoint(world, bodies[i], 'Ry')
else:
joints[i] = world.addJoint(bodies[i-1].end, bodies[i], 'Ry')
world.addVisualisation.Line(bodies[i], l)
# Add Damping
s = world.addSensor.Joint(symbol='state%i_'%i, joint=joints[i])
D = world.addExpression(symbol_str='T%i_'%i, exp=-d*s[1])
world.addLoad.Joint(joint=joints[i], symbol=D)
t = time.clock()
if method == "OrderN":
world.genEquations.OrderN()
elif method == "Recursive":
world.genEquations.Recursive()
elif method == "Explicit":
world.genEquations.Explicit()
print "Time needed for generating equations: %5.2f s" % (time.clock() - t)
world.show('Rope')
```