From 2c05db08262b869be2cdf2efed5facf91e87bb8e Mon Sep 17 00:00:00 2001 From: "armindamon.riess" <armindamon.riess@uzh.ch> Date: Tue, 25 Oct 2022 13:39:24 +0200 Subject: [PATCH] added some comments --- integrator.hpp | 4 ++++ main.cpp | 8 +++++++- parameters.txt | 2 +- plot.py | 1 + timeplot.py | 1 + vector.hpp | 1 + 6 files changed, 15 insertions(+), 2 deletions(-) diff --git a/integrator.hpp b/integrator.hpp index dae7469..018f7ea 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 9e2f1ce..86da633 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 39dc639..9585428 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 e9c4433..f24500b 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 2159878..00d50d1 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 a7441fe..b2f6a69 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; -- GitLab