pendulum c++ code

[ARFA Essosnam]

hello,
Can anyone tell me where I can find the c++ code of the point kinematics example of pendulum?

THANKS.

[Dimitar Stanev]

Hi,

Simbody is the underlying physics engine, while OpenSim implements additional features related to biomechanics. There is a Simbody example with a double pendulum.

https://github.com/simbody/simbody

Code: Select all

#include "Simbody.h"
using namespace SimTK;
int main() {
    // Define the system.
    MultibodySystem system;
    SimbodyMatterSubsystem matter(system);
    GeneralForceSubsystem forces(system);
    Force::Gravity gravity(forces, matter, -YAxis, 9.8);

    // Describe mass and visualization properties for a generic body.
    Body::Rigid bodyInfo(MassProperties(1.0, Vec3(0), UnitInertia(1)));
    bodyInfo.addDecoration(Transform(), DecorativeSphere(0.1));

    // Create the moving (mobilized) bodies of the pendulum.
    MobilizedBody::Pin pendulum1(matter.Ground(), Transform(Vec3(0)),
            bodyInfo, Transform(Vec3(0, 1, 0)));
    MobilizedBody::Pin pendulum2(pendulum1, Transform(Vec3(0)),
            bodyInfo, Transform(Vec3(0, 1, 0)));

    // Set up visualization.
    Visualizer viz(system);
    system.addEventReporter(new Visualizer::Reporter(viz, 0.01));

    // Initialize the system and state.
    State state = system.realizeTopology();
    pendulum2.setRate(state, 5.0);

    // Simulate for 20 seconds.
    RungeKuttaMersonIntegrator integ(system);
    TimeStepper ts(system, integ);
    ts.initialize(state);
    ts.stepTo(20.0);
}

Of course you can port this example very easily through OpenSim API, by changing the corresponding body, mobilizers and call the appropriate model functions to build the model.

[ARFA Essosnam]

thanks jim.

[Namrata Kaundal]

You can write a simple verlet integrator, which will allow you to do many more fun things easily (this is the home-made version of using a 2D Physics engine).




#include <iostream>
#include <cmath>

float const degrees = std::atan(1.0f)/45.0f;

struct Vector2D {
float x, y;
Vector2D(float x, float y) : x(x), y(y) {
}
};

std::ostream &operator<<(std::ostream &os, Vector2D v) {
return os << '(' << v.x << ',' << v.y << ')';
}

Vector2D operator+(Vector2D v, Vector2D w) {
return Vector2D(v.x+w.x, v.y+w.y);
}

Vector2D operator-(Vector2D v, Vector2D w) {
return Vector2D(v.x-w.x, v.y-w.y);
}

Vector2D operator*(Vector2D v, float s) {
return Vector2D(v.x*s, v.y*s);
}

float length(Vector2D v) {
return std::sqrt(v.x*v.x+v.y*v.y);
}

int main() {
float const dt = 0.01f; // (s) (100 Hz simulation)
float const dt_squared = dt*dt;
float length_of_string = 0.1f; // (m)
float gravity = 9.81f; // (m/s^2)
float friction = 1.5f; // (s^-1), I believe, but you better ask a Physicist if you care
float friction_coef = std::exp(-friction*dt);
float alpha = 15.0f*degrees;
Vector2D position(length_of_string*std::sin(alpha), -length_of_string*std::cos(alpha));
Vector2D previous_position = position;
Vector2D acceleration(0.0f, -gravity);

// For 10 seconds
for (float t = 0; t < 10.0; t += dt) {
// Verlet integration
Vector2D next_position = position + (position - previous_position) * friction_coef + acceleration * dt_squared;

// Constrain the point to be at the desired length from the origin
next_position = next_position * (length_of_string / length(next_position));

// Update
previous_position = position;
position = next_position;

// Output
std::cout << t << ' ' << position << '\n';
}