diff --git a/integrator.hpp b/integrator.hpp index dae746958b14077a18458d22231aabe07112f973..018f7ea592fa47eb2bd75c02e7bcf5f39469289c 100644 --- a/integrator.hpp +++ b/integrator.hpp @@ -14,12 +14,14 @@ public: Integrator() = delete; Integrator(double dt, F f) : dt_(dt), f_(f) {} ~Integrator() = default; + // save data to vector, four entries correspond to one time step void logData(std::vector<double>& data, const Planet& planet) { data.push_back(planet.r.x); data.push_back(planet.r.y); data.push_back(planet.energy()); data.push_back(planet.angularMomentum()); } + // perform integration with different algorithms void ExplicitEulerStep(Planet& planet) { logData(EEdat_, planet); Planet p0 = planet; @@ -54,6 +56,7 @@ public: planet.r = planet.r + dt_*planet.v; planet.v = planet.v + 0.5*dt_*f_(planet.r); } + // write data to file void print(std::ofstream& file) { for (unsigned i = 0; i < EEdat_.size(); ++i) { file << EEdat_[i] << " "; @@ -78,6 +81,7 @@ public: private: double dt_; F f_; + // vectors containing x position, y position, energy, and angular momentum std::vector<double> EEdat_; std::vector<double> RK2dat_; std::vector<double> RK4dat_; diff --git a/main.cpp b/main.cpp index 9e2f1ce5d74c0cac7d8a57b42da0c0ec4d28bdf8..86da63341426c21da26084c1a77ac9c34a2555e0 100644 --- a/main.cpp +++ b/main.cpp @@ -7,6 +7,7 @@ #include "planet.hpp" #include "integrator.hpp" +// force at a given position struct f { Vector operator()(const Vector pos) { double r = pos.norm(); @@ -37,12 +38,14 @@ int main(int argc, char *argv[]) { return 1; } + // read parameters std::string dummy; double dt, e; unsigned N; getline(file, dummy); file >> dt >> N >> e; + // initialize integrator Integrator<f> integrator(dt, f()); // initial conditions @@ -50,13 +53,14 @@ int main(int argc, char *argv[]) { Vector v0{0, std::sqrt(1+e)}; // initialize planets + // in hindsight, it would have been better to use a vector of planets and loop over them instead of repeating the same code for each one Planet EEplanet{r0, v0}; Planet RK2planet{r0, v0}; Planet RK4planet{r0, v0}; Planet SIplanet{r0, v0}; Planet LFplanet{r0, v0}; - // run simulations and measure execution times + // run simulations and measure computation times auto start = std::chrono::high_resolution_clock::now(); for (unsigned i = 0; i < N; ++i) { integrator.ExplicitEulerStep(EEplanet); @@ -92,10 +96,12 @@ int main(int argc, char *argv[]) { stop = std::chrono::high_resolution_clock::now(); auto LFtime = std::chrono::duration_cast<std::chrono::microseconds>(stop - start); + // save results std::ofstream datafile("./data.txt"); integrator.print(datafile); datafile.close(); + // print computation times std::cout << "Explicit Euler: " << EEtime.count()/1000.0 << " milliseconds" << std::endl << "RK2: " << RK2time.count()/1000.0 << " milliseconds" << std::endl << "RK4: " << RK4time.count()/1000.0 << " milliseconds" << std::endl diff --git a/parameters.txt b/parameters.txt index 39dc63916842d89bf1ebb1584016d031b48cd6c1..9585428eb227956a49ccecdd33ce1e0e5a561464 100644 --- a/parameters.txt +++ b/parameters.txt @@ -1,2 +1,2 @@ Time Step Number of Steps Excentricity -0.0005 1000 0.9 \ No newline at end of file +0.0005 1000000 0.9 \ No newline at end of file diff --git a/plot.py b/plot.py index e9c4433766856f8298dddb4720b2cc1c01bf56a8..f24500b014f60d6f884885a9f5c51434f97c3c35 100644 --- a/plot.py +++ b/plot.py @@ -4,6 +4,7 @@ import matplotlib.pyplot as plt # Read data from file data = np.loadtxt('build/data.txt') +# read parameters params = open('parameters.txt') params.readline() parameters = params.readline().split() diff --git a/timeplot.py b/timeplot.py index 215987856883aba3f5bccee9154990ba12f4047c..00d50d10f2094c32be7075fa6dcf7a461f11bb47 100644 --- a/timeplot.py +++ b/timeplot.py @@ -1,6 +1,7 @@ import matplotlib.pyplot as plt import numpy as np +# measured computation times data = [[0.129, 1.40, 14.8, 138], [0.154, 1.71, 19.0, 161], [0.241, 3.49, 27.1, 242], diff --git a/vector.hpp b/vector.hpp index a7441fe9a7f387cb1cf0c2ef1e63fa3903f7f9ee..b2f6a69d501c898c90d88f26f5e07888c761131d 100644 --- a/vector.hpp +++ b/vector.hpp @@ -3,6 +3,7 @@ #include <cmath> +// I probably should have used an existing library, but this was easier to implement. struct Vector { double x; double y;