The 3-body Problem

main.cpp

/* The 3-body Problem (actually n-body)
 * -------------------------------------
 *  
 * Copyright (C) 2025 Waldeck Schutzer <[email protected]> (@infinitymathart) 
 *  
 * This program is free software: you can redistribute it and/or modify  
 * it under the terms of the GNU General Public License as published by  
 * the Free Software Foundation, either version 3 of the License, or  
 * (at your option) any later version.  
 *  
 * This program 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 General Public License for more details.  
 *  
 * You should have received a copy of the GNU General Public License  
 * along with this program. If not, see <https://www.gnu.org/licenses/>.  
 */  


#include "ofMain.h"
#include "ofApp.h"

//========================================================================
int main( ){

	//Use ofGLFWWindowSettings for more options like multi-monitor fullscreen
	ofGLWindowSettings settings;
	settings.setSize(1024, 768);
	settings.windowMode = OF_WINDOW; //can also be OF_FULLSCREEN

	auto window = ofCreateWindow(settings);

	ofRunApp(window, make_shared<ofApp>());
	ofRunMainLoop();

}

ofApp.cpp

/* The 3-body Problem (actually n-body)
 * -------------------------------------
 *  
 * Copyright (C) 2025 Waldeck Schutzer <[email protected]> (@infinitymathart) 
 *  
 * This program is free software: you can redistribute it and/or modify  
 * it under the terms of the GNU General Public License as published by  
 * the Free Software Foundation, either version 3 of the License, or  
 * (at your option) any later version.  
 *  
 * This program 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 General Public License for more details.  
 *  
 * You should have received a copy of the GNU General Public License  
 * along with this program. If not, see <https://www.gnu.org/licenses/>.  
 */  

#include "ofApp.h"

//--------------------------------------------------------------
void ofApp::setup(){
    ofSetFrameRate(60);
    ofSetVerticalSync(true);
    ofSetBackgroundColor(0);
    
    // Initialize the three-body system
    threeBodySystem.setup();
}

//--------------------------------------------------------------
void ofApp::update(){
    threeBodySystem.update();
}

//--------------------------------------------------------------
void ofApp::draw(){
    ofEnableDepthTest();
    threeBodySystem.draw();
    
    // Optional: Display FPS in corner
    ofSetColor(255);
    string fps = "FPS: " + ofToString(ofGetFrameRate(), 0);
    ofDrawBitmapString(fps, 10, 20);
}

//--------------------------------------------------------------
void ofApp::keyPressed(int key){
    if(key == 'f') {
        ofToggleFullscreen();
    }
}

//--------------------------------------------------------------
void ofApp::keyReleased(int key){

}

//--------------------------------------------------------------
void ofApp::mouseMoved(int x, int y){

}

//--------------------------------------------------------------
void ofApp::mouseDragged(int x, int y, int button){

}

//--------------------------------------------------------------
void ofApp::mousePressed(int x, int y, int button){

}

//--------------------------------------------------------------
void ofApp::mouseReleased(int x, int y, int button){

}

//--------------------------------------------------------------
void ofApp::mouseEntered(int x, int y){

}

//--------------------------------------------------------------
void ofApp::mouseExited(int x, int y){

}

//--------------------------------------------------------------
void ofApp::windowResized(int w, int h){

}

//--------------------------------------------------------------
void ofApp::gotMessage(ofMessage msg){

}

//--------------------------------------------------------------
void ofApp::dragEvent(ofDragInfo dragInfo){ 

}

ofApp.h

/* The 3-body Problem (actually n-body)
 * -------------------------------------
 *  
 * Copyright (C) 2025 Waldeck Schutzer <[email protected]> (@infinitymathart) 
 *  
 * This program is free software: you can redistribute it and/or modify  
 * it under the terms of the GNU General Public License as published by  
 * the Free Software Foundation, either version 3 of the License, or  
 * (at your option) any later version.  
 *  
 * This program 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 General Public License for more details.  
 *  
 * You should have received a copy of the GNU General Public License  
 * along with this program. If not, see <https://www.gnu.org/licenses/>.  
 */  

#pragma once

#include "ofMain.h"
#include "ThreeBodySystem.h"

class ofApp : public ofBaseApp {
public:
    void setup();
    void update();
    void draw();
    
    void keyPressed(int key);
    void keyReleased(int key);
    void mouseMoved(int x, int y);
    void mouseDragged(int x, int y, int button);
    void mousePressed(int x, int y, int button);
    void mouseReleased(int x, int y, int button);
    void mouseEntered(int x, int y);
    void mouseExited(int x, int y);
    void windowResized(int w, int h);
    void dragEvent(ofDragInfo dragInfo);
    void gotMessage(ofMessage msg);

private:
    ThreeBodySystem threeBodySystem;
};

ThreeBodySystem.cpp

/* The 3-body Problem (actually n-body)
 * -------------------------------------
 *  
 * Copyright (C) 2025 Waldeck Schutzer <[email protected]> (@infinitymathart) 
 *  
 * This program is free software: you can redistribute it and/or modify  
 * it under the terms of the GNU General Public License as published by  
 * the Free Software Foundation, either version 3 of the License, or  
 * (at your option) any later version.  
 *  
 * This program 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 General Public License for more details.  
 *  
 * You should have received a copy of the GNU General Public License  
 * along with this program. If not, see <https://www.gnu.org/licenses/>.  
 */  

#include "ThreeBodySystem.h"

// IAS15 constants
const int G_IAS15 = 7;  // Number of Gauss-Radau spacings
const double H_IAS15[8] = {0.0, 0.0562625605369221464656521910318, 0.180240691736892364987579942780,
                          0.352624717113169637373907769648, 0.547153626330555383001448554766,
                          0.734210177215410531523210605558, 0.885320946839095768090359771030,
                          0.977520613561287501891174488626};

void ThreeBodySystem::setup() {
    // Initialize bodies
    bodies.resize(NUM_BODIES);
    
    // Figure-8 initial conditions
    // These values are scaled and adjusted for better visualization
    const double scale = 3.0;  // Spatial scale
    const double vScale = 0.35; // Velocity scale
    
    // Body 1 (Red)
    bodies[0].position = glm::dvec3(0.97000436, -0.24308753, 0.0) * scale;
    bodies[0].velocity = glm::dvec3(-0.93240737, -0.86473146, 0.0) * vScale;
    bodies[0].color.set(255, 50, 50);
    bodies[0].sphere.setRadius(0.05);
    
    // Body 2 (Yellow)
    bodies[1].position = glm::dvec3(-0.97000436, 0.24308753, 0.0) * scale;
    bodies[1].velocity = glm::dvec3(-0.93240737, -0.86473146, 0.0) * vScale;
    bodies[1].color.set(255, 255, 50);
    bodies[1].sphere.setRadius(0.05);
    
    // Body 4 (Blue)
    if (NUM_BODIES>=3)
    {
        // Body 3 (Blue)
        bodies[2].position = glm::dvec3(0.0, 0.0, 0.0) * scale;
        bodies[2].velocity = glm::dvec3(1.86481474, 1.72946292, 0.0) * vScale;
        bodies[2].color.set(50, 50, 255);
        bodies[2].sphere.setRadius(0.05);
    }

    // Body 4 (Blue)
    if (NUM_BODIES>=4)
    {
      bodies[3].position = glm::dvec3(0.0, 0.0, 3.0) * scale;
      bodies[3].velocity = glm::dvec3(1.44481474, -0.86473146, -0.4) * vScale;
      bodies[3].color.set(66, 66, 86);
      bodies[3].sphere.setRadius(0.035);
      bodies[3].mass = 0.035;
    }


    if (NUM_BODIES>=5)
    {
      for(int k=4; k<NUM_BODIES; k++)
      {
        bodies[k].position = glm::dvec3(ofRandom(-1.0,1.0), ofRandom(-1.0,1.0), ofRandom(-1.0,1.0)) * scale;
        bodies[k].velocity = glm::dvec3(ofRandom(-1.0,1.0), ofRandom(-1.0,1.0), ofRandom(-1.0,1.0)) * vScale;
        bodies[k].color.set(ofColor::fromHsb(ofMap(k,4,NUM_BODIES-1,0,255),200,255));
        bodies[k].sphere.setRadius(0.008);
        bodies[3].mass = 0.0025;

      }
    }

    // Initialize stars
    int numStars = ofRandom(1000, 2000);
    stars.resize(numStars);
    for(auto& star : stars) {
        star.position = glm::dvec3(
            ofRandom(-100, 100),
            ofRandom(-100, 100),
            ofRandom(-100, 100)
        );
        star.size = ofRandom(1, 4);
    }
    
    // Setup camera
    cam.setDistance(5);
    cam.setNearClip(0.1);
    cam.setFarClip(1000);
    
    // Initialize center of mass tracking
    calculateCenterOfMass();
    smoothedCenterOfMass = centerOfMass;
}


void ThreeBodySystem::updatePhysicsIAS15() {
    static constexpr double MAX_FORCE = 10.0; // Force limiting
    
    // Arrays for storing intermediate values
    glm::dvec3 k[G_IAS15+1][NUM_BODIES];
    glm::dvec3 r[NUM_BODIES];
    glm::dvec3 v[NUM_BODIES];
    
    // Store initial positions and velocities
    for(int i = 0; i < NUM_BODIES; i++) {
        r[i] = bodies[i].position;
        v[i] = bodies[i].velocity;
    }
    
    // Predictor step
    for(int i = 0; i < NUM_BODIES; i++) {
        k[0][i] = calculateAcceleration(bodies, i);
        
        // Limit acceleration magnitude
        double accMag = k[0][i].length();
        if(accMag > MAX_FORCE) {
            k[0][i] *= MAX_FORCE / accMag;
        }
        
        bodies[i].position = r[i] + v[i] * DT + k[0][i] * (DT * DT / 2);
        bodies[i].velocity = v[i] + k[0][i] * DT;
        
        // Limit velocity magnitude
        double velMag = bodies[i].velocity.length();
        static constexpr double MAX_VELOCITY = 5.0;
        if(velMag > MAX_VELOCITY) {
            bodies[i].velocity *= MAX_VELOCITY / velMag;
        }
    }
    
    // Corrector iterations
    for(int j = 1; j <= G_IAS15; j++) {
        // Calculate new accelerations
        for(int i = 0; i < NUM_BODIES; i++) {
            k[j][i] = calculateAcceleration(bodies, i);
            
            // Limit acceleration magnitude
            double accMag = k[j][i].length();
            if(accMag > MAX_FORCE) {
                k[j][i] *= MAX_FORCE / accMag;
            }
        }
        
        // Update positions and velocities
        for(int i = 0; i < NUM_BODIES; i++) {
            glm::dvec3 dr(0,0,0), dv(0,0,0);
            
            for(int l = 0; l <= j; l++) {
                double c = H_IAS15[l] * DT;
                dr += k[l][i] * (c * c / 2);
                dv += k[l][i] * c;
            }
            
            bodies[i].position = r[i] + v[i] * DT + dr;
            bodies[i].velocity = v[i] + dv;
            
            // Limit velocity magnitude
            double velMag = bodies[i].velocity.length();
            static constexpr double MAX_VELOCITY = 5.0;
            if(velMag > MAX_VELOCITY) {
                bodies[i].velocity *= MAX_VELOCITY / velMag;
            }
        }
    }
}

glm::dvec3 ThreeBodySystem::calculateAcceleration(const vector<Body>& bodies, int bodyIndex) {
    glm::dvec3 acc(0,0,0);
    static constexpr double SOFTENING2 = SOFTENING; // Softening parameter squared
    
    for(int j = 0; j < NUM_BODIES; j++) {
        if(j != bodyIndex) {
            glm::dvec3 r = bodies[j].position - bodies[bodyIndex].position;
            double distance2 = glm::length2(r) + SOFTENING2; // Softened distance squared
            double distance = sqrt(distance2);
            
            // Calculate force with smoothing near zero
            double forceMagnitude = G * bodies[j].mass / distance2;
            
            // Smooth transition for very close encounters
/*
            if(distance < MIN_DISTANCE) {
                forceMagnitude *= (distance / MIN_DISTANCE);
            }
 */           
            acc += glm::normalize(r) * forceMagnitude;
        }
    }
    
    return acc;
}


void ThreeBodySystem::updateTrails() {
    for(auto& body : bodies) {
        // Add current position to trail
        body.trail.push_back(body.position);
        
        // Remove oldest position if trail is too long
        while(body.trail.size() > TRAIL_LENGTH) {
            body.trail.erase(body.trail.begin());
        }
    }
}

void ThreeBodySystem::draw() {
    ofEnableDepthTest();
    ofBackground(0);
    
    cam.begin();
    //cam.setPosition(smoothedCenterOfMass + glm::dvec3(0, 0, 10));
    cam.lookAt(smoothedCenterOfMass);
    
    // Draw stars
    ofPushStyle();
    ofSetColor(255, 255, 255, 100);
    for(const auto& star : stars) {
        ofDrawCircle(glm::vec3(star.position), star.size * 0.01);
    }
    ofPopStyle();
    
    // Draw bodies and trails
    for(const auto& body : bodies) {
        // Draw trail
        if(body.trail.size() > 1) {
            ofMesh trail;
            trail.setMode(OF_PRIMITIVE_LINE_STRIP);
            
            for(size_t i = 0; i < body.trail.size(); i++) {
                double alpha = (double)i / body.trail.size();
                ofColor trailColor = body.color;
                trailColor.a = alpha * 255;
                trail.addColor(trailColor);
                trail.addVertex(body.trail[i]);
            }
            
            ofEnableAlphaBlending();
            trail.draw();
            ofDisableAlphaBlending();
        }

        // Add glow effect
        const float glowsize = pow(1.40,1.0/10);
        for(int k=0; k<10; k++)
        {
            ofPushStyle();
            ofEnableBlendMode(OF_BLENDMODE_ADD);
            ofSetColor(body.color, 150-10*k);
            double glowSize = body.sphere.getRadius() * pow(glowsize,k);
            ofDrawSphere(glm::vec3(body.position), glowSize);
            ofDisableBlendMode();
            ofPopStyle();
        }

        // Draw sphere
        ofPushMatrix();
        ofTranslate(glm::vec3(body.position));
        ofSetColor(body.color);
        ofEnableBlendMode(OF_BLENDMODE_DISABLED);
        body.sphere.draw();
                
        ofPopMatrix();
    }
    
    cam.end();
    ofDisableDepthTest();
}

void ThreeBodySystem::update() {
    // Run physics 5 times per frame for faster animation
    double temp_dt = DT;
    int oversamp = 240;
    DT = temp_dt / oversamp;
    for(int i = 0; i < 8*oversamp; i++) {
        updatePhysicsIAS15();
    }
    DT = temp_dt;
    updateTrails();
    calculateCenterOfMass();
    
    // Smooth camera movement
    smoothedCenterOfMass = smoothedCenterOfMass * SMOOTH_FACTOR + 
                          centerOfMass * (1.0f - SMOOTH_FACTOR);
}

void ThreeBodySystem::calculateCenterOfMass() {
    centerOfMass = glm::dvec3(0, 0, 0);
    double totalMass = 0;
    
    for(const auto& body : bodies) {
        centerOfMass += body.position * body.mass;
        totalMass += body.mass;
    }
    
    if(totalMass > 0) {
        centerOfMass /= totalMass;
    }
}

ThreeBodySystem.h

#pragma once
#include "ofMain.h"

class ThreeBodySystem {
public:
    void setup();
    void update();
    void draw();

private:
    // Bodies
    struct Body {
        glm::dvec3 position;
        glm::dvec3 velocity;
        ofColor color;
        ofSpherePrimitive sphere;
        vector<glm::dvec3> trail;
        double mass = 1.0f;
        void adjustMass(double factor) { mass *= factor; }
    };
    
// Add these constants to ThreeBodySystem.h in the private section:
static constexpr double MAX_VELOCITY = 10.0f;
static constexpr double MIN_DISTANCE = 0.0001f;
static constexpr double SOFTENING = 0.25f;

    vector<Body> bodies;
    static constexpr int NUM_BODIES = 3;
    static constexpr int TRAIL_LENGTH = 240; // 3 seconds at 60 FPS
    
    // Camera and movement
    ofEasyCam cam;
    glm::dvec3 centerOfMass;
    glm::dvec3 smoothedCenterOfMass;
    static constexpr double SMOOTH_FACTOR = 0.95f;
    
    // Stars background
    struct Star {
        glm::dvec3 position;
        double size;
    };
    vector<Star> stars;
    
    // Physics
    static constexpr double G = 1.0f; // Gravitational constant
    double DT = 0.001f; // Time step
    
    void updatePhysics();
    void updateTrails();
    glm::dvec3 calculateForce(const Body& b1, const Body& b2);
    void calculateCenterOfMass();
    // Add to private section of ThreeBodySystem.h:
    void updatePhysicsIAS15();
    glm::dvec3 calculateAcceleration(const vector<Body>& bodies, int bodyIndex);

};