////
Copyright (c) 2024 The C++ Alliance, Inc. (https://cppalliance.org)

Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

Official repository: https://github.com/boostorg/website-v2-docs
////
= Real-Time Simulation
:navtitle: Simulation

Creating a real-time simulation of objects or processes involves various aspects, including accurate timing, physical modeling, collisions and deformation, concurrent programming for real-time response, data storage, networking for multi-vehicle simulation, and usually a graphic interface though sometimes logging results is enough.

[square]
* <<Libraries>>
* <<Simulate Movement in 3D Space>>
* <<Add Collision Detection>>
* <<Deformation on Impact>>
* <<Blowing in the Wind>>
* <<See Also>>

== Libraries

Here are some libraries that may be helpful:

[circle]

* boost:geometry[]: For spatial computations and geometric algorithms, which you will likely need for modeling the physical behavior and interactions of your vehicles.

* boost:units[]: Helps with calculations involving units of measurement. It provides classes and functions that can enforce the correct usage of units and conversions between them, which could be helpful in a physical simulation.
 
* boost:algorithm[] : Provides a variety of utilities for numerical and string processing.

* boost:graph[]: In case you need to represent roads or pathways as a graph, this library provides a flexible and powerful way to represent and manipulate graphs. It also includes a number of graph algorithms.

* boost:chrono[]: Timing is critical in real-time applications. This library can help you measure time intervals, which could be useful for controlling the timing of your simulation.

* boost:thread[]: To achieve real-time performance, you might need to make use of multi-threading. This library provides classes and functions for multi-threading, synchronization, and inter-thread communication. 

* boost:interprocess[]: If you need to share data between different processes in real-time, this library can be useful. It supports shared memory, memory-mapped files, semaphores, and more.

* boost:mpi[] or boost:asio[]: For distributed simulations that run across multiple systems, you might need a library for network communication. boost:mpi[] provides a pass:[C++] interface for the Message Passing Interface (MPI) standard for distributed computing. boost:asio[] can also handle networking tasks and it is a bit lower-level.

* boost:serialization[]: To save the state of the simulation or to communicate complex data structures over a network, you might find this library helpful.

* boost:log[]: Supports severity levels, which you can use to categorize and filter your log messages. This can help you control the amount of log output and focus on what's important.

Note:: The code in this tutorial was written and tested using Microsoft Visual Studio (Visual C++ 2022, Console App project) with Boost version 1.88.0.

== Simulate Movement in 3D Space

The following sample simulates the movement of a geometric shape (a 3D box) in space using boost:geometry[] for shape representation, boost:chrono[] for timing, and boost:algorithm[] for numerical adjustments, such as scaling, normalization and smoothing movement (avoiding unrealistic motion).

The simulation itself runs for just a few seconds.

[source,cpp]
----
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/chrono.hpp>
#include <boost/algorithm/clamp.hpp>
#include <boost/thread/thread.hpp>

namespace bg = boost::geometry;
using namespace boost::chrono;

// Define a 3D box with min and max corners
using point3d = bg::model::point<double, 3, bg::cs::cartesian>;
using box3d = bg::model::box<point3d>;

const int xcoord = 0;
const int ycoord = 1;
const int zcoord = 2;
const double increment = 0.2;       // Guideline time increment in seconds

// Retrieve a point3d co-ordinate
double getCoord(int i, point3d p)
{
    switch (i)
    {
    case 0: return bg::get<xcoord>(p);
    case 1: return bg::get<ycoord>(p);
    case 2: return bg::get<zcoord>(p);
    }
}

// Set a point3d co-ordinate
void setCoord(int i, point3d& p, double v)
{
    switch (i)
    {
    case 0: bg::set<xcoord>(p, v);
        break;

    case 1: bg::set<ycoord>(p, v);
        break;

    case 2: bg::set<zcoord>(p, v);
        break;

    }
}

// Update the position of the box
void move_box(box3d& box, const point3d& velocity, double delta_time) {

    // Temporary variables for use in the movement calculations
    point3d min_corner = box.min_corner();
    point3d max_corner = box.max_corner();

    // Update the box position using velocity * delta time
    for (int i = 0; i < 3; ++i) {

        double new_min = getCoord(i,min_corner) + (getCoord(i,velocity) * delta_time);
        double new_max = getCoord(i,max_corner) + (getCoord(i,velocity) * delta_time);

        // Example clamp to avoid excessive movement (Boost.Algorithm)
        new_min = boost::algorithm::clamp(new_min, -100.0, 100.0);
        new_max = boost::algorithm::clamp(new_max, -100.0, 100.0);

        setCoord(i, min_corner, new_min);
        setCoord(i, max_corner, new_max);
    }

    // Update the box with new corners
    box = box3d(min_corner, max_corner);
}

int main() {

    // Set the sleep duration based on the guideline increment in seconds
    boost::chrono::duration<double> sleep_duration(increment);

    // Create a 3D box (min corner and max corner)
    box3d box(point3d(0.0, 0.0, 0.0), point3d(1.0, 1.0, 1.0));

    // Define velocity (units per second)
    point3d velocity(0.5, 0.3, -0.2);

    // Start timing
    steady_clock::time_point start_time = steady_clock::now();
    steady_clock::time_point current_time;
    steady_clock::time_point previous_time = start_time;
    double elapsed_seconds, increment_seconds;

    // Run simulation for 50 increments
    for (int i = 0; i < 50; ++i) {

        // Measure elapsed time - both from the start, and from the previous move
        current_time = steady_clock::now();
        elapsed_seconds = duration<double>(current_time - start_time).count();
        increment_seconds = duration<double>(current_time - previous_time).count();
        previous_time = current_time;

        // Move the box every increment, noting the increment will vary by tiny fractions of a second each time.
        move_box(box, velocity, increment_seconds);

        // Print times and the updated box position
        point3d min_corner = box.min_corner();
        point3d max_corner = box.max_corner();
        std::cout << "Time: " << elapsed_seconds << " sec | "
            << "Inc: " << increment_seconds << " sec | "
            << "Box Position: Min("
            << getCoord(xcoord,min_corner) << ", "
            << getCoord(ycoord,min_corner) << ", "
            << getCoord(zcoord,min_corner) << ") "
            << " Max("
            << getCoord(xcoord, max_corner) << ", "
            << getCoord(ycoord, max_corner) << ", "
            << getCoord(zcoord, max_corner) << ")\n";

        boost::this_thread::sleep_for(sleep_duration);       
    }

    return 0;
}

----

An example of the output:

[source,text]
----
Time: 0.209428 sec | Inc: 0.209424 sec | Box Position: Min(0.104714, 0.0628285, -0.0418857)  Max(1.10471, 1.06283, 0.958114)
Time: 0.429642 sec | Inc: 0.220214 sec | Box Position: Min(0.214821, 0.128893, -0.0859284)  Max(1.21482, 1.12889, 0.914072)
Time: 0.648115 sec | Inc: 0.218473 sec | Box Position: Min(0.324058, 0.194435, -0.129623)  Max(1.32406, 1.19443, 0.870377)
......

----

Note:: The time increment varies slightly on each loop, and this value is used when calculating movement.

== Add Collision Detection

Most 3D simulations require collision detection, which usually has a significant impact on the performance of a simulation, particularly in three dimensions. We'll introduce a bounding volume (a larger 3D box representing the environment), and detect when our moving box collides with its boundaries.

Collision detection is handled by checking if the box's min/max corners exceed the bounds. Some "bounce" mechanics are added to invert velocity after impact. In this example, the box moves continuously, rebounding off the walls, without consequences!

[source,cpp]
----
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/chrono.hpp>
#include <boost/algorithm/clamp.hpp>
#include <boost/thread/thread.hpp>

namespace bg = boost::geometry;
using namespace boost::chrono;

// Define a 3D box with min and max corners
using point3d = bg::model::point<double, 3, bg::cs::cartesian>;
using box3d = bg::model::box<point3d>;

const int xcoord = 0;
const int ycoord = 1;
const int zcoord = 2;
const double increment = 0.2;       // Guideline time increment in seconds

// Retrieve a point3d co-ordinate
double getCoord(int i, point3d p)
{
    switch (i)
    {
    case 0: return bg::get<xcoord>(p);
    case 1: return bg::get<ycoord>(p);
    case 2: return bg::get<zcoord>(p);
    }
}

// Set a point3d co-ordinate
void setCoord(int i, point3d& p, double v)
{
    switch (i)
    {
    case 0: bg::set<xcoord>(p, v);
        break;

    case 1: bg::set<ycoord>(p, v);
        break;

    case 2: bg::set<zcoord>(p, v);
        break;

    }
}

// Update the position of the box
void move_box(box3d& box, const point3d& velocity, double delta_time) {

    // Temporary variables for use in the movement calculations
    point3d min_corner = box.min_corner();
    point3d max_corner = box.max_corner();

    // Update the box position using velocity * delta time
    for (int i = 0; i < 3; ++i) {

        double new_min = getCoord(i,min_corner) + (getCoord(i,velocity) * delta_time);
        double new_max = getCoord(i,max_corner) + (getCoord(i,velocity) * delta_time);

        // Example clamp to avoid excessive movement (Boost.Algorithm)
        new_min = boost::algorithm::clamp(new_min, -100.0, 100.0);
        new_max = boost::algorithm::clamp(new_max, -100.0, 100.0);

        setCoord(i, min_corner, new_min);
        setCoord(i, max_corner, new_max);
    }

    // Update the box with new corners
    box = box3d(min_corner, max_corner);
}

// Function to check and handle collisions with the bounding box
void handle_collision(box3d& box, point3d& velocity, const box3d& bounds) {
    for (int i = 0; i < 3; ++i) {
        double min_pos = getCoord(i, box.min_corner());
        double max_pos = getCoord(i, box.max_corner());
        double bound_min = getCoord(i, bounds.min_corner());
        double bound_max = getCoord(i, bounds.max_corner());

        // If box collides with environment limits, reverse velocity
        if (min_pos <= bound_min || max_pos >= bound_max) {
            setCoord(i, velocity, -getCoord(i, velocity));              // Reverse direction
            double impact_force = std::abs(getCoord(i,velocity));       // Impact force = velocity along the impact axis
           std::cout << "\nCollision with impact force: " << impact_force << "\n"
                << "New velocity : ("
                << getCoord(xcoord, velocity) << ", "
                << getCoord(ycoord, velocity) << ", "
                << getCoord(zcoord, velocity) << ")\n\n";
        }
    }
}


int main() {

    // Set the sleep duration based on the guideline increment in seconds
    boost::chrono::duration<double> sleep_duration(increment);

    // Define the 3D simulation space (bounding box)
    box3d bounds(point3d(-5.0, -5.0, -5.0), point3d(5.0, 5.0, 5.0));

    // Create a 3D box (min corner and max corner)
    box3d box(point3d(0.0, 0.0, 0.0), point3d(1.0, 1.0, 1.0));

    // Define velocity (units per second)
    point3d velocity(0.5, 0.3, -0.2);

    // Start timing
    steady_clock::time_point start_time = steady_clock::now();
    steady_clock::time_point current_time;
    steady_clock::time_point previous_time = start_time;
    double elapsed_seconds, increment_seconds;

    // Run simulation for 80 increments
    for (int i = 0; i < 80; ++i) {

        // Measure elapsed time - both from the start, and from the previous move
        current_time = steady_clock::now();
        elapsed_seconds = duration<double>(current_time - start_time).count();
        increment_seconds = duration<double>(current_time - previous_time).count();
        previous_time = current_time;

        // Move the box every increment, noting the increment will vary by tiny fractions of a second each time.
        move_box(box, velocity, increment_seconds);

        // Check for collision
        handle_collision(box, velocity, bounds);

        // Print times and the updated box position
        point3d min_corner = box.min_corner();
        point3d max_corner = box.max_corner();
        std::cout << "Time: " << elapsed_seconds << " sec | "
            << "Inc: " << increment_seconds << " sec | "
            << "Box Position: Min("
            << getCoord(xcoord,min_corner) << ", "
            << getCoord(ycoord,min_corner) << ", "
            << getCoord(zcoord,min_corner) << ") "
            << " Max("
            << getCoord(xcoord, max_corner) << ", "
            << getCoord(ycoord, max_corner) << ", "
            << getCoord(zcoord, max_corner) << ")\n";

        boost::this_thread::sleep_for(sleep_duration);       
    }

    return 0;
}

----

Note:: An unlikely event perhaps, but a collision between the box and two or three sides of the bounding area will be handled by this code.


An example of the output:

[source,text]
----
Time: 7.76744 sec | Inc: 0.2183 sec | Box Position: Min(3.88372, 2.33023, -1.55349)  Max(4.88372, 3.33023, -0.553488)
Time: 7.98306 sec | Inc: 0.215622 sec | Box Position: Min(3.99153, 2.39492, -1.59661)  Max(4.99153, 3.39492, -0.596613)

Collision with impact force: 0.5
New velocity : (-0.5, 0.3, -0.2)

Time: 8.19834 sec | Inc: 0.21528 sec | Box Position: Min(4.09917, 2.4595, -1.63967)  Max(5.09917, 3.4595, -0.639669)
Time: 8.41678 sec | Inc: 0.218435 sec | Box Position: Min(3.98995, 2.52503, -1.68336)  Max(4.98995, 3.52503, -0.683356)
......

----

== Deformation on Impact

Collisions rarely have no consequences. To simulate deformation on impact, we can modify the shape of the box when it collides with a boundary.

image::deformation.png[cube deforming in 3D space]

The box will deform along the axis of impact, based on the velocity.

[source,cpp]
----
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/chrono.hpp>
#include <boost/algorithm/clamp.hpp>
#include <boost/thread/thread.hpp>

namespace bg = boost::geometry;
using namespace boost::chrono;

// Define a 3D box with min and max corners
using point3d = bg::model::point<double, 3, bg::cs::cartesian>;
using box3d = bg::model::box<point3d>;

const int xcoord = 0;
const int ycoord = 1;
const int zcoord = 2;
const double increment = 0.2;       // Guideline time increment in seconds

// Retrieve a point3d co-ordinate
double getCoord(int i, point3d p)
{
    switch (i)
    {
    case 0: return bg::get<xcoord>(p);
    case 1: return bg::get<ycoord>(p);
    case 2: return bg::get<zcoord>(p);
    }
}

// Set a point3d co-ordinate
void setCoord(int i, point3d& p, double v)
{
    switch (i)
    {
    case 0: bg::set<xcoord>(p, v);
        break;

    case 1: bg::set<ycoord>(p, v);
        break;

    case 2: bg::set<zcoord>(p, v);
        break;

    }
}

// To verify deformation, calculate the volume of the box
double volume(const box3d& box) {
    double dx = getCoord(xcoord, box.max_corner()) - getCoord(xcoord, box.min_corner());
    double dy = getCoord(ycoord, box.max_corner()) - getCoord(ycoord, box.min_corner());
    double dz = getCoord(zcoord, box.max_corner()) - getCoord(zcoord, box.min_corner());
    return dx * dy * dz;
}

// Function to deform the box upon collision
void deform_box(box3d& box, int axis, double impact_force) {
    // Get the box corners
    point3d min_corner = box.min_corner();
    point3d max_corner = box.max_corner();

    // Deformation ratio (scales based on impact)
    double deformation = 1.0 - (impact_force * 0.2);
    deformation = boost::algorithm::clamp(deformation, 0.7, 1.0);  // Prevent over-deformation

    // Scale the box on the axis of impact
    double center = (getCoord(axis,min_corner) + getCoord(axis,max_corner)) / 2.0;
    setCoord(axis,min_corner, center - (center - getCoord(axis,min_corner)) * deformation);
    setCoord(axis,max_corner, center + (getCoord(axis,max_corner) - center) * deformation);

    // Update the box
    box = box3d(min_corner, max_corner);
}

// Update the position of the box
void move_box(box3d& box, const point3d& velocity, double delta_time) {

    // Temporary variables for use in the movement calculations
    point3d min_corner = box.min_corner();
    point3d max_corner = box.max_corner();

    // Update the box position using velocity * delta time
    for (int i = 0; i < 3; ++i) {

        double new_min = getCoord(i,min_corner) + (getCoord(i,velocity) * delta_time);
        double new_max = getCoord(i,max_corner) + (getCoord(i,velocity) * delta_time);

        // Example clamp to avoid excessive movement (Boost.Algorithm)
        new_min = boost::algorithm::clamp(new_min, -100.0, 100.0);
        new_max = boost::algorithm::clamp(new_max, -100.0, 100.0);

        setCoord(i, min_corner, new_min);
        setCoord(i, max_corner, new_max);
    }

    // Update the box with new corners
    box = box3d(min_corner, max_corner);
}

// Function to check and handle collisions with the bounding box
void handle_collision(box3d& box, point3d& velocity, const box3d& bounds) {
    for (int i = 0; i < 3; ++i) {
        double min_pos = getCoord(i, box.min_corner());
        double max_pos = getCoord(i, box.max_corner());
        double bound_min = getCoord(i, bounds.min_corner());
        double bound_max = getCoord(i, bounds.max_corner());

        // If box collides with environment limits, reverse velocity
        if (min_pos <= bound_min || max_pos >= bound_max) {
            setCoord(i, velocity, -getCoord(i, velocity));          // Reverse direction
            double impact_force = std::abs(getCoord(i,velocity));   // Higher velocity = more deformation
            deform_box(box, i, impact_force);                       // Apply deformation
            std::cout << "\nCollision with impact force: " << impact_force << "\n"
                << "New velocity : ("
                << getCoord(xcoord, velocity) << ", "
                << getCoord(ycoord, velocity) << ", "
                << getCoord(zcoord, velocity) << "), Volume of box: "
                << volume(box) << "\n\n";
        }
    }
}

int main() {

    // Set the sleep duration based on the guideline increment in seconds
    boost::chrono::duration<double> sleep_duration(increment);

    // Define the 3D simulation space (bounding box)
    box3d bounds(point3d(-5.0, -5.0, -5.0), point3d(5.0, 5.0, 5.0));

    // Create a 3D box (min corner and max corner)
    box3d box(point3d(0.0, 0.0, 0.0), point3d(1.0, 1.0, 1.0));

    // Define velocity (units per second)
    point3d velocity(0.5, 0.3, -0.2);

    // Start timing
    steady_clock::time_point start_time = steady_clock::now();
    steady_clock::time_point current_time;
    steady_clock::time_point previous_time = start_time;
    double elapsed_seconds, increment_seconds;

    // Run simulation for 80 increments
    for (int i = 0; i < 80; ++i) {

        // Measure elapsed time - both from the start, and from the previous move
        current_time = steady_clock::now();
        elapsed_seconds = duration<double>(current_time - start_time).count();
        increment_seconds = duration<double>(current_time - previous_time).count();
        previous_time = current_time;

        // Move the box every increment, noting the increment will vary by tiny fractions of a second each time.
        move_box(box, velocity, increment_seconds);

        // Check for collision
        handle_collision(box, velocity, bounds);

        // Print times and the updated box position
        point3d min_corner = box.min_corner();
        point3d max_corner = box.max_corner();
        std::cout << "Time: " << elapsed_seconds << " sec | "
            << "Inc: " << increment_seconds << " sec | "
            << "Box Position: Min("
            << getCoord(xcoord,min_corner) << ", "
            << getCoord(ycoord,min_corner) << ", "
            << getCoord(zcoord,min_corner) << ") "
            << " Max("
            << getCoord(xcoord, max_corner) << ", "
            << getCoord(ycoord, max_corner) << ", "
            << getCoord(zcoord, max_corner) << ")\n";

        boost::this_thread::sleep_for(sleep_duration);       
    }

    return 0;
}

----

An example of the output:

[source,text]
----
Time: 7.40482 sec | Inc: 0.20316 sec | Box Position: Min(3.70241, 2.22145, -1.48096)  Max(4.70241, 3.22145, -0.480964)
Time: 7.62483 sec | Inc: 0.22001 sec | Box Position: Min(3.81242, 2.28745, -1.52497)  Max(4.81242, 3.28745, -0.524966)
Time: 7.8286 sec | Inc: 0.203764 sec | Box Position: Min(3.9143, 2.34858, -1.56572)  Max(4.9143, 3.34858, -0.565719)

Collision with impact force:0.5
New velocity : (-0.5, 0.3, -0.2), Volume of box: 0.9

Time: 8.03252 sec | Inc: 0.203922 sec | Box Position: Min(4.06626, 2.40976, -1.6065)  Max(4.96626, 3.40976, -0.606503)
Time: 8.25219 sec | Inc: 0.219674 sec | Box Position: Min(3.95642, 2.47566, -1.65044)  Max(4.85642, 3.47566, -0.650438)

----

== Blowing in the Wind

Let's look at something more natural and how it might be simulated. Say you have a real-time simulation that is out in the open, and you need to realistically model the effects of wind on your objects.

This is one of those cases where a little mathematics goes a long way. Real wind is not random. It has structure across multiple time scales. Meteorologists often think of wind as having four components:

. Mean wind (the forecast strength and direction)
. Long-period variation (30 to 120 second changes)
. Gusts (2 to 10 second bursts)
. Turbulence (small, rapid fluctuations)

Avoid the temptation of using a random number every frame, you can model each component separately.

Say the weather forecast is for a 10 mile an hour wind, from the NNW, gusting to 20 mph. Let's look at the four components.

. Mean wind : this should be converted to a vector (speed x direction), and direction is fixed at NNW (337.5 degrees).
. Long-period variation : Real wind slowly strengthens and weakens. A sine wave works surprisingly well, something like:
+
[source,cpp]
----
double long_period(double t)
{
    return 1.5 *
        sin(2.0 * pi * t / 90.0);
}
----

. Gusts : Gusting to 20 mph according to our forecast. We do not want sudden jumps, instead smoother pulses.
+
[source,cpp]
----
double gust(double t)
{
    return exp(-x*x);
}
----
+
To be natural, each gust might last 4 to 8 seconds.

. Turbulence : We need to avoid obvious repetition, which we can manage with several sine waves, though the frequencies are unrelated:
+
[source,cpp]
----
double turbulence(double t)
{
    return
        0.4*sin(5.3*t)
      + 0.3*sin(8.9*t)
      + 0.2*sin(13.7*t);
}
----

Total wind speed will now look something like this:

[source,cpp]
----
double wind_speed(double t)
{
    return mean_wind + long_period(t) + gust_strength(t) + turbulence(t);
}
----

Wind direction is of course not _fixed_, but will wander left and right, perhaps plus or minus 5 degrees:

[source,cpp]
----
double wind_direction(double t)
{
    return mean_direction + 4*sin(t/40) + 2*sin(t/11) + 1*sin(7*t);
}
----

Using the constants available in boost:math[], and the units from boost:units[], a first cut at a wind simulator might look like the following. Here we have three types of object - leaves, paper, tumbleweed - each with certain physical properties that determine how quickly the wind will move them. 

[source,cpp]
----
#include <boost/math/constants/constants.hpp>

#include <boost/units/systems/si.hpp>
#include <boost/units/io.hpp>

#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>

namespace si = boost::units::si;

using boost::units::quantity;

constexpr double pi =
boost::math::constants::pi<double>();
//------------------------------------------------------------
// Small-scale turbulence.
//
// Combines several sine waves of different frequencies to
// produce smooth, non-repeating fluctuations.
//------------------------------------------------------------

double turbulence(double t)
{
    return

        0.30 * std::sin(7.3 * t)

        + 0.20 * std::sin(11.7 * t + 0.8)

        + 0.12 * std::sin(17.9 * t + 2.1)

        + 0.08 * std::sin(29.4 * t);
}
//------------------------------------------------------------
// Wind speed as a function of time.
//
// Returns metres per second.
//------------------------------------------------------------

quantity<si::velocity>
wind_speed(quantity<si::time> t)
{
    // Mean wind = 4.47 m/s (≈10 mph)

    quantity<si::velocity> mean =
        4.47 * si::meters_per_second;

    // Smooth variation

    double variation =
        0.9 *
        std::sin(
            2.0 * pi *
            t.value() / 60.0);

    // Gaussian gust

    double x =
        (t.value() - 30.0) / 5.0;

    double gust =
        4.5 * std::exp(-x * x);

    double turbulent =
        turbulence(t.value());

    return
        mean
        + variation * si::meters_per_second
        + gust * si::meters_per_second
        + turbulent * si::meters_per_second;
}

struct Object
{
    std::string name;

    quantity<si::mass> mass;

    double drag;

    double friction;

    quantity<si::length> position;

    quantity<si::velocity> velocity;
};

int main()
{
    std::vector<Object> objects =
    {
        {
            "Leaf",

            0.005 * si::kilograms,

            0.90,

            0.02,

            0.0 * si::meters,

            0.0 * si::meters_per_second
        },

        {
            "Paper",

            0.010 * si::kilograms,

            0.80,

            0.05,

            0.0 * si::meters,

            0.0 * si::meters_per_second
        },

        {
            "Tumbleweed",

            0.300 * si::kilograms,

            0.60,

            0.08,

            0.0 * si::meters,

            0.0 * si::meters_per_second
        }
    };

    quantity<si::time> dt =
        1.0 * si::seconds;

    std::cout
        << std::fixed
        << std::setprecision(2);

    for (quantity<si::time> t = 0.0 * si::seconds;
        t <= 60.0 * si::seconds;
        t += dt)
    {
        auto wind = wind_speed(t);

        std::cout
            << "\nTime = "
            << t
            << "   Wind = "
            << wind
            << "\n";

        for (auto& o : objects)
        {
            //------------------------------------------------
            // Very simple physics
            //------------------------------------------------

            quantity<si::force> force =
                wind * o.drag *
                o.mass /
                (1.0 * si::seconds);

            // Simple damping

            force -=
                o.friction *
                o.velocity *
                o.mass /
                (1.0 * si::seconds);

            quantity<si::acceleration>
                acceleration =
                force / o.mass;

            o.velocity +=
                acceleration * dt;

            o.position +=
                o.velocity * dt;

            std::cout
                << std::setw(12)
                << o.name
                << "   Position = "
                << o.position
                << "   Velocity = "
                << o.velocity
                << '\n';
        }
    }
}

----

Run this program and you will get a lot of numbers!

[source,text]
----
...

Time = 58.00 s   Wind = 4.63 m s^-1
        Leaf   Position = 6073.07 m   Velocity = 155.59 m s^-1
       Paper   Position = 3533.71 m   Velocity = 70.33 m s^-1
  Tumbleweed   Position = 1911.14 m   Velocity = 32.32 m s^-1

Time = 59.00 s   Wind = 4.37 m s^-1
        Leaf   Position = 6229.48 m   Velocity = 156.41 m s^-1
       Paper   Position = 3604.02 m   Velocity = 70.31 m s^-1
  Tumbleweed   Position = 1943.50 m   Velocity = 32.36 m s^-1

Time = 60.00 s   Wind = 4.06 m s^-1
        Leaf   Position = 6386.41 m   Velocity = 156.94 m s^-1
       Paper   Position = 3674.07 m   Velocity = 70.04 m s^-1
  Tumbleweed   Position = 1975.71 m   Velocity = 32.20 m s^-1
----

For extra realism, treat gusts like moving weather cells instead of global events. Imagine a gust as a "blob" of stronger wind drifting across the landscape. The gust has a position, radius, speed, and intensity. The result is that waves of motion travel naturally through your environment, much like you see when watching a real forest in a breeze or wind gusts crossing a lake. It looks dramatically more realistic than having every object receive the same gust at the same instant, and the underlying mathematics is still quite simple: each object samples the wind field at its own position.

== Next Steps

It is quite difficult to assess the "naturalness" of a physical process just by looking at numbers. To really assess how well our wind simulator works, we would need to have graphical output (not necessarily a full-featured 3D rendering) where we could visualize the different objects blowing across the terrain, answering the question "does it look right?".

It is good practice when designing a simulation of real-world activity to clearly define what is to be simulated and what is not. All simulations are simplifications to an extent, though they do tend to be large and challenging programs to write. A complex simulation might have several processes running on different threads. For a sample of multi-threading code, refer to xref:task-parallel-computation.adoc[].

== See Also

* https://www.boost.org/doc/libs/latest/libs/libraries.htm#Algorithms[Category: Algorithms]
* https://www.boost.org/doc/libs/latest/libs/libraries.htm#Image-processing[Category: Image processing]
* https://www.boost.org/doc/libs/latest/libs/libraries.htm#Math[Category: Math and numerics]
* https://www.boost.org/doc/libs/latest/libs/libraries.htm#State[Category: State Machines]
