N-body problem: Difference between revisions

Content added Content deleted
Line 428:
===Further work===
And as in every experiment report/journal paper/whatever, time for some philosophizing, a more rigorous implementation will require the bodies to be realistic, this means they will not be point particles and hence the dynamics will be much more complex. Tidal effects, collisions, orbits will require much more computation. When I say realistic, I don't mean rigid bodies because even rigid bodies are an idealization of reality. And all of this, even with classical mechanics, employing the Newtonian version of Gravitation will be a worthy task for any competent Physicist, let alone a Programmer. Introduce the General Theory of Relativity and everything changes. I can go on and on, but this is a coding, not a rambling site. I will end by saying that the Inverse Square Law, a name often used for Newtonian Gravitation, is just one class of forces studied in Dynamics. The force relation can take any form, the N-body problem can therefore be solved in (countably/uncountably ?) infinite ways.
 
=={{header|C++}}==
{{trans|Java}}
<lang cpp>#include <algorithm>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
 
class Vector {
private:
double px, py, pz;
 
public:
Vector() : px(0.0), py(0.0), pz(0.0) {
// empty
}
 
Vector(double x, double y, double z) : px(x), py(y), pz(z) {
// empty
}
 
double mod() const {
return sqrt(px*px + py * py + pz * pz);
}
 
Vector operator+(const Vector& rhs) const {
return Vector(px + rhs.px, py + rhs.py, pz + rhs.pz);
}
 
Vector operator-(const Vector& rhs) const {
return Vector(px - rhs.px, py - rhs.py, pz - rhs.pz);
}
 
Vector operator*(double s) const {
return Vector(px*s, py*s, pz*s);
}
 
bool operator==(const Vector& rhs) const {
return px == rhs.px
&& py == rhs.py
&& pz == rhs.pz;
}
 
friend std::istream& operator>>(std::istream&, Vector&);
friend std::ostream& operator<<(std::ostream&, Vector&);
};
 
std::istream& operator>>(std::istream& in, Vector& v) {
return in >> v.px >> v.py >> v.pz;
}
 
std::ostream& operator<<(std::ostream& out, Vector& v) {
auto precision = out.precision();
auto width = out.width();
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.px << " ";
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.py << " ";
out << std::fixed << std::setw(width) << std::setprecision(precision) << v.pz;
return out;
}
 
const Vector ORIGIN{ 0.0, 0.0, 0.0 };
 
class NBody {
private:
double gc;
int bodies;
int timeSteps;
std::vector<double> masses;
std::vector<Vector> positions;
std::vector<Vector> velocities;
std::vector<Vector> accelerations;
 
void resolveCollisions() {
for (int i = 0; i < bodies; ++i) {
for (int j = i + 1; j < bodies; ++j) {
if (positions[i] == positions[j]) {
std::swap(velocities[i], velocities[j]);
}
}
}
}
 
void computeAccelerations() {
for (int i = 0; i < bodies; ++i) {
accelerations[i] = ORIGIN;
for (int j = 0; j < bodies; ++j) {
if (i != j) {
double temp = gc * masses[j] / pow((positions[i] - positions[j]).mod(), 3);
accelerations[i] = accelerations[i] + (positions[j] - positions[i]) * temp;
}
}
}
}
 
void computeVelocities() {
for (int i = 0; i < bodies; ++i) {
velocities[i] = velocities[i] + accelerations[i];
}
}
 
void computePositions() {
for (int i = 0; i < bodies; ++i) {
positions[i] = positions[i] + velocities[i] + accelerations[i] * 0.5;
}
}
 
public:
NBody(std::string& fileName) {
using namespace std;
 
ifstream ifs(fileName);
if (!ifs.is_open()) {
throw runtime_error("Could not open file.");
}
 
ifs >> gc >> bodies >> timeSteps;
 
masses.resize(bodies);
positions.resize(bodies);
fill(positions.begin(), positions.end(), ORIGIN);
velocities.resize(bodies);
fill(velocities.begin(), velocities.end(), ORIGIN);
accelerations.resize(bodies);
fill(accelerations.begin(), accelerations.end(), ORIGIN);
 
for (int i = 0; i < bodies; ++i) {
ifs >> masses[i] >> positions[i] >> velocities[i];
}
 
cout << "Contents of " << fileName << '\n';
cout << gc << ' ' << bodies << ' ' << timeSteps << '\n';
for (int i = 0; i < bodies; ++i) {
cout << masses[i] << '\n';
cout << positions[i] << '\n';
cout << velocities[i] << '\n';
}
cout << "\nBody : x y z | vx vy vz\n";
}
 
int getTimeSteps() {
return timeSteps;
}
 
void simulate() {
computeAccelerations();
computePositions();
computeVelocities();
resolveCollisions();
}
 
friend std::ostream& operator<<(std::ostream&, NBody&);
};
 
std::ostream& operator<<(std::ostream& out, NBody& nb) {
for (int i = 0; i < nb.bodies; ++i) {
out << "Body " << i + 1 << " : ";
out << std::setprecision(6) << std::setw(9) << nb.positions[i];
out << " | ";
out << std::setprecision(6) << std::setw(9) << nb.velocities[i];
out << '\n';
}
return out;
}
 
int main() {
std::string fileName = "nbody.txt";
NBody nb(fileName);
 
for (int i = 0; i < nb.getTimeSteps(); ++i) {
std::cout << "\nCycle " << i + 1 << '\n';
nb.simulate();
std::cout << nb;
}
 
return 0;
}</lang>
{{out}}
<pre>Contents of nbody.txt
0.01 3 20
1
0.000000 0.000000 0.000000
0.010000 0.000000 0.000000
0.100000
1.000000 1.000000 0.000000
0.000000 0.000000 0.020000
0.001000
0.000000 1.000000 1.000000
0.010000 -0.010000 -0.010000
 
Body : x y z | vx vy vz
 
Cycle 1
Body 1 : 0.010177 0.000179 0.000002 | 0.010354 0.000357 0.000004
Body 2 : 0.998230 0.998232 0.020002 | -0.003539 -0.003536 0.020004
Body 3 : 0.010177 0.988232 0.988055 | 0.010354 -0.013536 -0.013889
 
Cycle 2
Body 1 : 0.020709 0.000718 0.000011 | 0.010710 0.000721 0.000014
Body 2 : 0.992907 0.992896 0.039971 | -0.007109 -0.007138 0.019935
Body 3 : 0.020717 0.972888 0.972173 | 0.010727 -0.017153 -0.017876
 
Cycle 3
Body 1 : 0.031600 0.001625 0.000034 | 0.011072 0.001094 0.000033
Body 2 : 0.983985 0.983910 0.059834 | -0.010735 -0.010835 0.019790
Body 3 : 0.031643 0.953868 0.952235 | 0.011125 -0.020886 -0.021999
 
Cycle 4
Body 1 : 0.042858 0.002913 0.000081 | 0.011443 0.001481 0.000060
Body 2 : 0.971393 0.971163 0.079509 | -0.014448 -0.014659 0.019561
Body 3 : 0.042981 0.931039 0.928087 | 0.011552 -0.024772 -0.026299
 
Cycle 5
Body 1 : 0.054492 0.004595 0.000160 | 0.011826 0.001884 0.000097
Body 2 : 0.955030 0.954509 0.098909 | -0.018278 -0.018649 0.019238
Body 3 : 0.054766 0.904225 0.899522 | 0.012018 -0.028857 -0.030829
 
Cycle 6
Body 1 : 0.066517 0.006691 0.000281 | 0.012224 0.002308 0.000145
Body 2 : 0.934759 0.933760 0.117931 | -0.022265 -0.022849 0.018806
Body 3 : 0.067040 0.873197 0.866280 | 0.012530 -0.033199 -0.035655
 
Cycle 7
Body 1 : 0.078950 0.009225 0.000456 | 0.012642 0.002759 0.000206
Body 2 : 0.910400 0.908677 0.136456 | -0.026454 -0.027316 0.018244
Body 3 : 0.079856 0.837662 0.828023 | 0.013101 -0.037871 -0.040861
 
Cycle 8
Body 1 : 0.091815 0.012227 0.000702 | 0.013086 0.003245 0.000284
Body 2 : 0.881722 0.878958 0.154340 | -0.030902 -0.032122 0.017523
Body 3 : 0.093281 0.797239 0.784313 | 0.013749 -0.042975 -0.046559
 
Cycle 9
Body 1 : 0.105140 0.015737 0.001035 | 0.013564 0.003775 0.000383
Body 2 : 0.848429 0.844216 0.171401 | -0.035684 -0.037362 0.016600
Body 3 : 0.107405 0.751427 0.734579 | 0.014498 -0.048649 -0.052908
 
Cycle 10
Body 1 : 0.118964 0.019805 0.001481 | 0.014085 0.004362 0.000509
Body 2 : 0.810137 0.803953 0.187408 | -0.040900 -0.043166 0.015414
Body 3 : 0.122346 0.699554 0.678056 | 0.015384 -0.055097 -0.060138
 
Cycle 11
Body 1 : 0.133337 0.024498 0.002071 | 0.014662 0.005025 0.000672
Body 2 : 0.766343 0.757509 0.202050 | -0.046687 -0.049720 0.013868
Body 3 : 0.138268 0.640690 0.613685 | 0.016460 -0.062633 -0.068603
 
Cycle 12
Body 1 : 0.148327 0.029907 0.002851 | 0.015317 0.005792 0.000888
Body 2 : 0.716377 0.703998 0.214889 | -0.053246 -0.057302 0.011810
Body 3 : 0.155406 0.573482 0.539941 | 0.017816 -0.071782 -0.078886
 
Cycle 13
Body 1 : 0.164025 0.036157 0.003887 | 0.016079 0.006709 0.001184
Body 2 : 0.659310 0.642172 0.225282 | -0.060887 -0.066351 0.008976
Body 3 : 0.174112 0.495836 0.454475 | 0.019596 -0.083511 -0.092045
 
Cycle 14
Body 1 : 0.180564 0.043437 0.005286 | 0.017000 0.007852 0.001613
Body 2 : 0.593807 0.570186 0.232208 | -0.070119 -0.077621 0.004875
Body 3 : 0.194929 0.404136 0.353320 | 0.022038 -0.099890 -0.110265
 
Cycle 15
Body 1 : 0.198150 0.052049 0.007234 | 0.018171 0.009372 0.002283
Body 2 : 0.517817 0.485100 0.233878 | -0.081861 -0.092550 -0.001535
Body 3 : 0.218605 0.290860 0.228583 | 0.025314 -0.126661 -0.139210
 
Cycle 16
Body 1 : 0.217126 0.062542 0.010117 | 0.019781 0.011614 0.003484
Body 2 : 0.427899 0.381659 0.226654 | -0.097974 -0.114332 -0.012913
Body 3 : 0.244268 0.131956 0.057562 | 0.026013 -0.191148 -0.202831
 
Cycle 17
Body 1 : 0.238346 0.076539 0.015221 | 0.022658 0.016380 0.006723
Body 2 : 0.317489 0.248502 0.200967 | -0.122846 -0.151982 -0.038461
Body 3 : 0.075592 -0.559591 -0.487315 | -0.363366 -1.191945 -0.886924
 
Cycle 18
Body 1 : 0.263123 0.097523 0.026918 | 0.026898 0.025587 0.016672
Body 2 : 0.173428 0.050424 0.112716 | -0.165275 -0.244174 -0.138041
Body 3 : -0.286241 -1.745597 -1.369528 | -0.360299 -1.180066 -0.877501
 
Cycle 19
Body 1 : 0.270854 0.113045 0.061923 | -0.011436 0.005457 0.053339
Body 2 : 0.199821 -0.093105 -0.208666 | 0.218061 -0.042882 -0.504723
Body 3 : -0.646318 -2.924909 -2.246453 | -0.359856 -1.178559 -0.876350
 
Cycle 20
Body 1 : 0.258572 0.116046 0.112038 | -0.013129 0.000544 0.046890
Body 2 : 0.426346 -0.111425 -0.681150 | 0.234987 0.006241 -0.440245
Body 3 : -1.006089 -4.103186 -3.122591 | -0.359686 -1.177995 -0.875924</pre>
 
=={{header|C#|C sharp}}==
Retrieved from "https://rosettacode.org/wiki/N-body_problem"