Atomic simulation

atomic.cpp

// Required dependencies:
// * SFML
// * Some font, place it in res/font.ttf

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
using namespace std;

#include <SFML/Graphics.hpp>
using namespace sf;

// Constants
const int sx=599,sy=599,fq=120;
const double dtr = 57.2957795;
const double MAX_SPAWN_AREA = 10000, MAX_RADIUS = MAX_SPAWN_AREA / 1.41;
const int MAX_BODIES = 300;
const double FITNESS_PERIOD = 120;
const double MIN_VEL_TO_GET_ENERGY = 20;

double mpl = 1.0;
sf::Vector2f camPos;
double last_am = 0, last_bm = 0, last_density = 0;
std::ofstream stat_file;

struct PeriodicStats
{
    size_t deaths = 0;
    double death_rate = 25;

    size_t collisions = 0;
    double collision_rate = 10;

    size_t vel_overruns = 0;
    double vel_overrun_rate = 5;

    size_t pos_overruns = 0;
    double pos_overrun_rate = 10;

    size_t energy_overruns = 0;
    double energy_overrun_rate = 10;
} stats;

struct Body
{
    Body(double px, double py, double r);

    double vecX, vecY, posX, posY, fX, fY, radius, density, mass, am, bm;
    double energy = 10;
    bool dead = false;

    void update();
    void moveUpdate();
};

std::vector<Body*> bodies;

double dist(Body& _1, Body& _2)
{
    double x = _2.posX - _1.posX;
    double y = _2.posY - _1.posY;
    return sqrt(x*x+y*y);
}

double angle(Body& _s, Body& _e)
{
    double x = _s.posX - _e.posX;
    double y = _s.posY - _e.posY;
    return atan2(x,y) * dtr;
}

Body::Body(double px, double py, double r)
{
    while(density < 0.1 || isnan(density) || !isfinite(density) || density >= 1000) density = last_density + ((double)rand() / RAND_MAX - 0.5) * stats.death_rate / 10.0;
    last_density = density;
    while(am <= 1 || isnan(am) || !isfinite(am) || am <= 0 || am >= 1000) am = last_am + ((double)rand() / RAND_MAX - 0.5) * stats.death_rate / 10.0;
    last_am = am;
    while(bm <= 1 || isnan(bm) || !isfinite(bm) || bm <= 0 || bm >= 1000) bm = last_bm + ((double)rand() / RAND_MAX - 0.5) * stats.death_rate / 10.0;
    last_bm = bm;

    posX = px;
    posY = py;
    radius = r / 2.f;
    mass = density * radius * radius * 3.14;
    fX=0; fY=0;
    vecX = (double(rand()) / RAND_MAX - 0.5) * 10;
    vecY = (double(rand()) / RAND_MAX - 0.5) * 10;

    std::cout << "new body [" << px << "," << py << "]" << std::endl;
}

void Body::update()
{
    fX = 0.f;
    fY = 0.f;

    // interaction with other bodies
    // TODO: don't interact with everything !!!
    for(unsigned i = 0; i < bodies.size(); i++)
    {
        Body* body = bodies[i];

        if(body != this)
        {
            double r = dist(*body, *this);

            if(2 * r < body->radius + radius)
            {
                std::cout << "body dead: collided; r_this=" << radius << ", r_body=" << body->radius << ", r_diff=" << r << std::endl;
                body->dead = true;
                dead = true;
                stats.collisions++;
                return;
            }

            double deg = angle(*body, *this);

            // interaction
            double G = 1e6;
            double A = G*mass*body->mass / am;
            double B = G*mass*body->mass / bm;
            double U = (B-A)/pow(r,2);
            double F = U / r;
            fX += F * sin(deg/dtr);
            fY += F * cos(deg/dtr);
        }
    }

    //cout << endl;
    //cout << "p= " << posX << ", " << posY << endl;
    //cout << "v= " << vecX << ", " << vecY << endl;
    //cout << "F= " << fX << ", " << fY << endl;
}
void Body::moveUpdate()
{
    posX += vecX / fq;
    posY += vecY / fq;

    vecX += fX / mass / fq;
    vecY += fY / mass / fq;

    // make something to keep particles moving if all particles are spawned
    energy += vecX*vecX+vecY*vecY - MIN_VEL_TO_GET_ENERGY;

    if(vecX*vecX+vecY*vecY > 1e8)
    {
        std::cout << "body dead: vec=[" << vecX << "," << vecY << "]" << std::endl;
        stats.vel_overruns++;
        dead = true;
    }
    if(posX*posX+posY*posY > MAX_RADIUS*MAX_RADIUS)
    {
        std::cout << "body dead: pos=[" << posX << "," << posY << "]" << std::endl;
        stats.pos_overruns++;
        dead = true;
    }
    if(bodies.size() == MAX_BODIES && energy <= 0)
    {
        std::cout << "body dead: energy overrun" << std::endl;
        stats.energy_overruns++;
        dead = true;
    }
}

void render(sf::RenderTarget& t, sf::Font& font)
{
    sf::CircleShape cs(MAX_RADIUS, 100);
    cs.setOrigin(MAX_RADIUS, MAX_RADIUS);

    cs.setOutlineThickness(1.f / mpl);
    cs.setOutlineColor(sf::Color::Red);
    cs.setFillColor(sf::Color::Transparent);
    t.draw(cs);

    sf::RectangleShape rs({MAX_SPAWN_AREA, MAX_SPAWN_AREA});
    rs.setOrigin(MAX_SPAWN_AREA/2,MAX_SPAWN_AREA/2);
    rs.setOutlineThickness(1.f / mpl);
    rs.setOutlineColor(sf::Color::Green);
    rs.setFillColor(sf::Color::Transparent);
    t.draw(rs);

    for(Body* body: bodies)
    {
        sf::CircleShape csb(body->radius);
        csb.setOrigin(body->radius, body->radius);
        csb.setPosition(body->posX, body->posY);
        csb.setOutlineColor(sf::Color::White);
        float spdmpl = std::min(1.0, (body->vecX * body->vecX + body->vecY * body->vecY) / 250000.0);
        sf::Color color(255 * spdmpl, 255 * spdmpl, 255 * spdmpl);
        csb.setFillColor(sf::Color::Red * color + sf::Color::Green * (sf::Color::White - color));
        csb.setOutlineThickness(1.f / mpl);
        t.draw(csb);


        std::ostringstream oss;
        oss << "mass: " << body->mass << std::endl;
        oss << "am: " << body->am << std::endl;
        oss << "bm: " << body->bm << std::endl;
        oss << "e: " << body->energy << std::endl;

        sf::Text text(oss.str(), font, 15);
        text.setPosition(body->posX, body->posY);
        t.draw(text);
    }
}

void renderGUI(sf::RenderTarget& t, sf::Font& font)
{
    std::ostringstream oss;
    oss << "--- environment ---" << std::endl;
    oss << "* bodies: " << bodies.size() << std::endl;
    oss << "* mpl: " << mpl << std::endl;
    oss << "* density: " << last_density << std::endl;
    oss << "* am: " << last_am << std::endl;
    oss << "* bm: " << last_bm << std::endl;
    oss << "--- stats ---" << std::endl;
    oss << "* dr: " << stats.death_rate << " (" << stats.deaths << ")" << std::endl;
    oss << "* cr: " << stats.collision_rate << " (" << stats.collisions << ")" << std::endl;
    oss << "* vor: " << stats.vel_overrun_rate << " (" << stats.vel_overruns << ")" << std::endl;
    oss << "* por: " << stats.pos_overrun_rate << " (" << stats.pos_overruns << ")" << std::endl;
    oss << "* eor: " << stats.energy_overrun_rate << " (" << stats.energy_overruns << ")" << std::endl;

    sf::Text text(oss.str(), font, 15);
    text.setPosition(5, 5);
    t.draw(text);
}

void tick()
{
    static size_t deathRateClock = 0;

    for(Body* body: bodies)
    {
        body->update();
    }
    for(Body* body: bodies)
    {
        body->moveUpdate();
    }

    // dead update
    decltype(bodies)::iterator lastIt = bodies.begin();
    for(auto it = lastIt; !bodies.empty() && it != bodies.end(); it++)
    {
        if((*it)->dead)
        {
            bodies.erase(it);
            it = lastIt;
            stats.deaths++;
            continue;
        }
        lastIt = it;
    }

    // body spawn
    if(bodies.size() < MAX_BODIES)
        bodies.push_back(new Body((double(rand()) / RAND_MAX - 0.5) * MAX_SPAWN_AREA, (double(rand()) / RAND_MAX - 0.5) * MAX_SPAWN_AREA, 50.f));

    // death rate
    if(deathRateClock > FITNESS_PERIOD)
    {
        deathRateClock = 0;

        stats.death_rate = (double)stats.deaths / FITNESS_PERIOD * 60;
        stats.deaths = 0;
        stats.collision_rate = (double)stats.collisions / FITNESS_PERIOD * 60;
        stats.collisions = 0;
        stats.vel_overrun_rate = (double)stats.vel_overruns / FITNESS_PERIOD * 60;
        stats.vel_overruns = 0;
        stats.pos_overrun_rate = (double)stats.pos_overruns / FITNESS_PERIOD * 60;
        stats.pos_overruns = 0;
        stats.energy_overrun_rate = (double)stats.energy_overruns / FITNESS_PERIOD * 60;
        stats.energy_overruns = 0;

        stat_file << bodies.size() << " " << stats.death_rate << " " << stats.collision_rate << " " << stats.vel_overrun_rate << " " << stats.pos_overrun_rate << " " << stats.energy_overrun_rate << std::endl;
    }
    deathRateClock++;
}

int main()
{
    sf::RenderWindow wnd(sf::VideoMode(sx,sy), "Atomic Physics");

    std::cout << "waiting 10 seconds" << std::endl;
    sleep(sf::seconds(10));

    //wnd.setFramerateLimit(fq);
    sf::View view(sf::Vector2f(sx / 2, sy / 2), sf::Vector2f(sx, sy));
    sf::View guiView(sf::FloatRect(0, 0, sx, sy));

    srand(time(nullptr));

    stat_file.open("stats.txt");
    if(stat_file.fail())
    {
        std::cout << "failed to open stat file for writing" << std::endl;
        return 1;
    }
    stat_file << "bodies.size()" << " " << "stats.death_rate" << " " << "stats.collision_rate" << " " << "stats.vel_overrun_rate" << " " << "stats.pos_overrun_rate" << " " << "stats.energy_overrun_rate" << std::endl;

    sf::Font font;
    if(!font.loadFromFile("res/font.ttf"))
        return 1;

    while(wnd.isOpen())
    {
        sf::Event e;
        while(wnd.pollEvent(e))
        {
            if(e.type == sf::Event::Closed)
                wnd.close();
            else if(e.type == sf::Event::Resized)
            {
                view.setSize(sf::Vector2f(e.size.width / mpl, e.size.height / mpl));
                guiView = sf::View(sf::FloatRect({}, sf::Vector2f(e.size.width, e.size.height)));
            }
            else if(e.type == sf::Event::MouseButtonPressed && e.mouseButton.button == sf::Mouse::Left)
            {
                sf::Vector2f pos = wnd.mapPixelToCoords(sf::Vector2i(e.mouseButton.x, e.mouseButton.y), view);
                bodies.push_back(new Body(pos.x, pos.y, 50.f));
            }
            else if(e.type == sf::Event::MouseWheelScrolled && e.mouseWheelScroll.wheel == sf::Mouse::VerticalWheel)
            {
                float d = pow(2, -e.mouseWheelScroll.delta);
                view.zoom(d);
                mpl /= d;
            }
        }

        if(sf::Keyboard::isKeyPressed(sf::Keyboard::Up)) {camPos.y--; view.move(0.f, -6.f / mpl);}
        if(sf::Keyboard::isKeyPressed(sf::Keyboard::Down)) {camPos.y++; view.move(0.f, 6.f / mpl);}
        if(sf::Keyboard::isKeyPressed(sf::Keyboard::Right)) {camPos.x++; view.move(6.f / mpl, 0.f);}
        if(sf::Keyboard::isKeyPressed(sf::Keyboard::Left)) {camPos.x--; view.move(-6.f / mpl, 0.f);}

        tick();
        wnd.clear();
        wnd.setView(view);
        render(wnd, font);
        wnd.setView(guiView);
        renderGUI(wnd, font);
        wnd.display();
    }

    for(Body* body: bodies)
        delete body;
}