In C++ rechnet es viel schneller. Übersetz werden kann mit Gnu-C oder Visual Studio.
Code:
/*
RKF5nbody - Kartesische Koordinaten
*/
#include <cstdio>
#include <string>
#include<vector>
#include<ctime>
#include<cmath>
#include<iostream>
#include<fstream>
// Visual Studio hat kein pi.
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#define MAIN
using namespace std;
// meine eigenen Bibliotheken
#include "getopt.h"
#include "tools.h"
#include <limits>
#include <iostream>
#include <iomanip>
/* Globale Variablen */
const double G = 6.67408e-11; // Gravitationskonstante
const double mProAE=149597870700; // Astronomische Einheit
const double sekProTag=86400;
const double sekProJahr=365.2425*86400;
const double Lj=9460730472580800; // in Metern
int N; // Anzahl der Körper
double S=1; // allg. Steuerfktor, Benutzung abhängig von der Aufgabe
// für Butcher-Tableau etwas weiter schematisiert
#define maxOrdnung 5 // Ordnung der Fehlberg-Methode
#define maxBodies 10000 // max. Anzahhl vorn Körpern
double r[maxOrdnung+1][maxBodies][3];
double v[maxOrdnung+1][maxBodies][3];
double a[maxOrdnung+1][maxBodies][3];
double m[maxBodies];
// für Protokoll
vector<double> rMin(maxBodies,myNAN);
vector<double> rMax(maxBodies,0);
bool IDE=false;
void heliozentrisch2baryzentrisch() { // ab zeile 366: Geht nur für die Sonne, nicht für nbody allg. S. coCo.
int i,j,k;
double gmasse=0;
for (i=1;i<N;i++) gmasse+=m[i];
gmasse+=m[0]; // Sonne zuletzt
for (k=0;k<3;k++) { // alle Koordinaten
for (j = 1; j < N; j++) { // Sonne korrigieren
r[0][0][k]+=m[j]*r[0][j][k];
v[0][0][k]+=m[j]*v[0][j][k];
}
r[0][0][k]=-r[0][0][k]/gmasse;
v[0][0][k]=-v[0][0][k]/gmasse;
for (j = 1; j < N; j++) { // Planeten ins Baryzentrum rücken
r[0][j][k]+=r[0][0][k];
v[0][j][k]+=v[0][0][k];
}
}
}
double init1() {
// für 2451545,0, aus xyzPlanet (Bahnelemente für 2000), heliozentrisch
const double kGauss = 0.01720209895; // Gravitationskonstante nach Gauß, 3. Keplersches Gesetz
// Die Sonnemasse aus der Gaußschen Gravitationskonstante ausrechnen,
// damit das Produkt aus G und MSun wieder stimmt
double MSun=sq(kGauss)/G *mProAE*mProAE*mProAE/sq(sekProTag);
N=10; // 9 Planeten und Sonne
// moderne Massen
m[0]=MSun;
m[1]=MSun/6023600;
m[2]=MSun/ 408523.5;
m[3]=MSun/ 328900.5;
m[4]=MSun/3098710;
m[5]=MSun/ 1047.355;
m[6]=MSun/ 3498.5;
m[7]=MSun/ 22869;
m[8]=MSun/ 19314;
m[9]=MSun/300000;
r[0][0][0]=0;
r[0][0][1]=0;
r[0][0][2]=0;
v[0][0][0]=0;
v[0][0][1]=0;
v[0][0][2]=0;
// Merkur
r[0][1][0]= -0.1300972207003186;
r[0][1][1]= -0.4472890198141531;
r[0][1][2]= -0.0245969133239213;
v[0][1][0]= 0.0213662324174032;
v[0][1][1]= -0.0064481960971869;
v[0][1][2]= -0.0024877977034037;
// Venus
r[0][2][0]= -0.7183166316917509;
r[0][2][1]= -0.0327032722771552;
r[0][2][2]= 0.0410147977488960;
v[0][2][0]= 0.0007987449513147;
v[0][2][1]= -0.0202947994390621;
v[0][2][2]= -0.0003234500580008;
// Erde
r[0][3][0]= -0.1771880023702050;
r[0][3][1]= 0.9672110887942617;
r[0][3][2]= -0.0000000000000000;
v[0][3][0]= -0.0172030714318618;
v[0][3][1]= -0.0031645602421189;
v[0][3][2]= 0.0000000000000000;
// Mars
r[0][4][0]= 1.3906336696396515;
r[0][4][1]= -0.0133808914031599;
r[0][4][2]= -0.0344572279098808;
v[0][4][0]= 0.0006725050100979;
v[0][4][1]= 0.0151881144907859;
v[0][4][2]= 0.0003016339332627;
// Jupiter
r[0][5][0]= 4.0011766985384991;
r[0][5][1]= 2.9385765042116430;
r[0][5][2]= -0.1017825822273280;
v[0][5][0]= -0.0045683141626712;
v[0][5][1]= 0.0064432061940609;
v[0][5][2]= 0.0000755804032107;
// Saturn
r[0][6][0]= 6.4064144387907085;
r[0][6][1]= 6.5699908996536998;
r[0][6][2]= -0.3690684852593678;
v[0][6][0]= -0.0042923475419510;
v[0][6][1]= 0.0038903154986327;
v[0][6][2]= 0.0001029466380039;
// Uranus
r[0][7][0]= 14.4318791761688789;
r[0][7][1]=-13.7343034133537500;
r[0][7][2]= -0.2381521766122678;
v[0][7][0]= 0.0026781013835866;
v[0][7][1]= 0.0026726984663166;
v[0][7][2]= -0.0000247729502902;
// Neptun
r[0][8][0]= 16.8120566822566566;
r[0][8][1]=-24.9917817543021776;
r[0][8][2]= 0.1272293952807222;
v[0][8][0]= 0.0025792757608357;
v[0][8][1]= 0.0017768971425853;
v[0][8][2]= -0.0000959102440067;
// Pluto
r[0][9][0]= -9.8753669819861276;
r[0][9][1]=-27.9588031341511325;
r[0][9][2]= 5.8504355694029950;
v[0][9][0]= 0.0030287515246100;
v[0][9][1]= -0.0015377848009245;
v[0][9][2]= -0.0007122076484756;
// Orte umrechnen in Meter
for (int i = 0; i < N; i++) {
for (int k=0;k<3;k++) r[0][i][k] *=mProAE;
}
// Geschwindigkeiten umrechnen in m/s
for (int i = 0; i < N; i++) {
for (int k=0;k<3;k++) v[0][i][k] *=mProAE/sekProTag;
}
heliozentrisch2baryzentrisch() ;
return(2451545.0*sekProTag);
}
double init2() { // baryzentrisch
// für 2451545,0, aus xyzPlanet (Bahnelemente für 2000), heliozentrisch
const double kGauss = 0.01720209895; // Gravitationskonstante nach Gauß, 3. Keplersches Gesetz
double MSun=sq(kGauss)/G *mProAE*mProAE*mProAE/sq(sekProTag);
N=10; // 9 Planeten und Sonne
r[0][0][0]= -1063436181.5;
r[0][0][1]= -403178315.5;
r[0][0][2]= 27924311.6;
r[0][1][0]= -20525703382.2;
r[0][1][1]= -67316663267.2;
r[0][1][2]= -3651721547.4;
r[0][2][0]= -108522074771.0;
r[0][2][1]= -5295518213.1;
r[0][2][2]= 6163650722.0;
r[0][3][0]= -27570384049.7;
r[0][3][1]= 144289541085.5;
r[0][3][2]= 27924311.6;
r[0][4][0]= 206972399720.3;
r[0][4][1]= -2404931177.5;
r[0][4][2]= -5126803613.9;
r[0][5][0]= 597504078214.3;
r[0][5][1]= 439201609603.6;
r[0][5][2]= -15198533263.9;
r[0][6][0]= 957322522683.3;
r[0][6][1]= 982453470791.1;
r[0][6][2]= -55183935225.7;
r[0][7][0]= 2157914958773.0;
r[0][7][1]=-2055025724501.0;
r[0][7][2]= -35599134212.1;
r[0][8][0]= 2513984445571.8;
r[0][8][1]=-3739120513758.2;
r[0][8][2]= 19061170936.1;
r[0][9][0]=-1478397309067.7;
r[0][9][1]=-4182980594505.0;
r[0][9][2]= 875240628161.8;
v[0][0][0]= 9.29323348;
v[0][0][1]= -12.81658147;
v[0][0][2]= -0.15918117;
v[0][1][0]= 37004.00242937;
v[0][1][1]= -11177.58979901;
v[0][1][2]= -4307.67352332;
v[0][2][0]= 1392.28564030;
v[0][2][1]= -35152.38582183;
v[0][2][2]= -560.19899547;
v[0][3][0]= -29777.08241127;
v[0][3][1]= -5492.11604817;
v[0][3][2]= -0.15918117;
v[0][4][0]= 1173.70663100;
v[0][4][1]= 26284.74809179;
v[0][4][2]= 522.10695479;
v[0][5][0]= -7900.54555615;
v[0][5][1]= 11143.31683419;
v[0][5][2]= 130.70502469;
v[0][6][0]= -7422.72126397;
v[0][6][1]= 6723.09678597;
v[0][6][2]= 178.08847902;
v[0][7][0]= 4646.31018375;
v[0][7][1]= 4614.84545075;
v[0][7][2]= -43.05247532;
v[0][8][0]= 4475.19788359;
v[0][8][1]= 3063.80412436;
v[0][8][2]= -166.22362888;
v[0][9][0]= 5253.44576775;
v[0][9][1]= -2675.42458857;
v[0][9][2]= -1233.31598335;
// moderne Massen
m[0]=1.9884754159665367161e+30;
m[1]=3.3011412045397049632e+23;
m[2]=4.8674688627864414779e+24;
m[3]=6.0458266739227718677e+24;
m[4]=6.4171071702951773445e+23;
m[5]=1.8985686953960564018e+27;
m[6]=5.6837942431514554871e+26;
m[7]=8.6950693776139604302e+25;
m[8]=1.0295513181974406248e+26;
m[9]=6.6282513865551224955e+24;
return(2451545.0*sekProTag);
}
// Abstand des Körpers q von Körper p
double abstand (double r[][3], int p, int q) {
return sqrt( sq(r[p][0]-r[q][0]) + sq(r[p][1]-r[q][1]) + sq(r[p][2]-r[q][2]) );
}
int anzBeschl=0;
void beschleunigung(int N, double m[], double r[][3], double a[][3]) {
anzBeschl++;
int p,q,k;
double R3, A;
// alle Beschleunigungen löschen
for (p = 0; p < N; p++) for (k = 0; k < 3; k++) a[p][k]=0;
for (p=0; p<N; p++) { // alle Planeten p
for (q=p+1; q<N; q++) { // gegen alle andern Planeten q
double R=abstand(r, p, q); R3=R*R*R;
for (k=0; k<3; k++) {
A=G*(r[q][k]-r[p][k])/R3;
a[p][k]+=m[q]*A;
a[q][k]-=m[p]*A;
}
}
} // Beschleunigungen löschen
}
double energie () {
/* Lokale Variablen */
int i,j;
double d;
/* Berechnung der kinetischen Energie */
double Ekin = 0;
for (i = 0;i < N;i ++) {
double v2=sq(v[0][i][0]) + sq(v[0][i][1]) + sq(v[0][i][2]);
Ekin += 0.5*m[i] * v2;
}
/* Berechnung der potentiellen Energie */
double Epot = 0;
for (i = 0;i < N;i ++) {
/* Beitrag pro Abstandspaar zur potentiellen Energie */
for (j = i + 1;j < N;j ++) {
d=abstand(r[0], i,j);
/* Potentielle Energie pro Abstandspaar */
Epot -= G * m[i] * m[j] / d;
}
}
return Ekin+Epot;
}
vector<double> impuls() {
double ImpSonne;
vector<double> Imp(5);
for (int i = 0; i < N; i++) {
Imp[0]+=m[i]*v[0][i][0];
Imp[1]+=m[i]*v[0][i][1];
Imp[2]+=m[i]*v[0][i][2];
if (0==i) ImpSonne=sqrt(sq(Imp[0]) + sq(Imp[1]) +sq(Imp[2]) );
}
Imp[3]=sqrt(sq(Imp[0]) + sq(Imp[1]) +sq(Imp[2]) );
Imp[4]=Imp[3]/ImpSonne;
return Imp;
}
vector<double> drehimpuls() {
vector<double> Dimp(4);
for (int i=0;i<N;i++) {
double a1=r[0][i][0];
double a2=r[0][i][1];
double a3=r[0][i][2];
double b1=v[0][i][0];
double b2=v[0][i][1];
double b3=v[0][i][2];
// Kreuzprodukt
Dimp[0]+=a2*b3-a3*b2;
Dimp[1]+=a3*b1-a1*b3;
Dimp[2]+=a1*b2-a2*b1;
}
Dimp[3]=sqrt(sq(Dimp[0]) + sq(Dimp[1]) +sq(Dimp[2]) );
return Dimp;
}
void ausgabe(double time) {
cout<<std::setprecision(20);
cout<<"zeit: "<<time<<"\n";
for (int p = 0; p < N; p++) {
cout<<abstand(r[0],p,0)/mProAE<<"\t";
for (int k=0;k<3;k++) cout<<(r[0][p][k]-r[0][0][k])/mProAE<<" ";
cout<<"\n";
}
}
void ausgabeHelioz() {
cout<<std::setprecision(20);
for (int p = 0; p < N; p++) {
double x=r[0][p][0]-r[0][0][0];
double y=r[0][p][1]-r[0][0][1];
double z=r[0][p][2]-r[0][0][2];
double r=sqrt(x*x+y*y+z*z);
double l=atan2(y,x); if (l<0) l+=2*M_PI;
printf("r=%lf, l=%lf, b=%lf\n",r/mProAE, 180*l/M_PI, 180*atan(z/r)/M_PI);
}
}
void printRMinMax(double time, double dt) {
vector<double> result=impuls();double Imp=result[3];double ImpR=result[4];
result=drehimpuls();double Dimp=result[3];
cout<<std::setprecision(20);
if (time < 1e-10) return; // keine MyNAN ausgeben: maxR ist so gesetzt
cout<<time/sekProTag/365.25<<" Jahre, SW:"<<dt/sekProTag<<" Tage. Energie: "<<energie()<<" Impuls: "<<Imp<<" "<<ImpR<<" D-Imp: "<< Dimp <<"\n";
cout << "rMin ";
for (int i = 1; i < N; i++) printf("%1.4lf ",rMin[i]/mProAE);
cout<<"\nrMax ";
for (int i = 1; i < N; i++) printf("%1.4lf ",rMax[i]/mProAE);
cout<<"\n";
for (int i = 1; i < N; i++) { // wieder löschen
rMin[i]=myNAN;
rMax[i]=0;
}
}
void protokoll(string I, double t, double p, double dt) {
if ("99"==I) return; // kein Protokoll beim letzten Schritt
if (I=="3") printRMinMax(t, dt); // Ausgabe der Sonnenabstände
else ausgabe(t/p); // Datei eingelesen
}
// Runge- Kutta-Fehlberg mit Butcher-Tableau
double rkf45(double tEnd, double dt, double p, string I) {
int ordnung=5;
double rErg1[maxBodies][3];
double rErg2[maxBodies][3];
double vErg1[maxBodies][3];
double vErg2[maxBodies][3];
bool swAenderbar=true;
if (dt < 0) { // Schrittweite<0 gegeben: Bleibt konstant
swAenderbar=false;
dt=fabs(dt);
}
// Rechnung in negative Richtung behandeln
if (tEnd < 0) {
tEnd=-tEnd;
dt=-dt;
}
// das Tableau ...
// ... Stützstellen für die Zeit, wird nicht benötigt
double sst[6]={0, 1.0/4, 3.0/8, 12.0/13, 1.0, 1.0/2};
// .. Koeffizienten , mit 0 aufgefüllt
double c[6][5]={
{0, 0, 0, 0, 0},
{1.0/4, 0, 0, 0, 0}, // 0.25
{3.0/32, 9.0/32, 0, 0, 0}, // 0.375
{1932.0/2197, -7200.0/2197, 7296.0/2197, 0, 0}, // 0.923
{439.0/216, -8.0, 3680.0/513, -845.0/4104, 0}, // 0.9998
{-8.0/27, 2.0, -3544.0/2565, 1859.0/4104, -11.0/40} // 0.5
};
// ... und die beiden Ergebnisse
double c1[6]={25.0/216, 0, 1408.0/2565, 2197.0/4104, -1.0/5, 0}; // Summe ist 1
double c2[6]={16.0/135, 0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55}; // Summe ist 1
int i,j,k, schritt;
double rErrMin, rErrMax;
// Beginn der Rechnung
double t=0;
double progressTime=0;
protokoll(I, t, p, dt);
while (t<tEnd) { // For-Schleife geht hier nicht
t+=fabs(dt);
progressTime+=fabs(dt);
for (schritt=0;schritt<ordnung+1;schritt++) {
for (i = 0; i < N;i ++) { // alle Körper
for (k = 0; k < 3;k ++) { // alle Koordinaten
r[schritt][i][k]=r[0][i][k]; // Beginn immer mit dem alten Ergebnis
v[schritt][i][k]=v[0][i][k];
for (j=0; j<schritt; j++) { // alle k_i
r[schritt][i][k] += c[schritt][j] * dt * v[j][i][k] ;
v[schritt][i][k] += c[schritt][j] * dt * a[j][i][k] ;
}
}
}
beschleunigung (N, m, r[schritt], a[schritt]);
} // Schritt des Fehlberg-Verfahrens
// Fehlberg-Ergebnisse berechnen
for (i = 0; i < N;i ++) {
for (k = 0;k < 3;k ++) {
// Beginn immer mit dem alten Ergebnis
rErg1[i][k]=r[0][i][k];
vErg1[i][k]=v[0][i][k];
rErg2[i][k]=r[0][i][k];
vErg2[i][k]=v[0][i][k];
for (j=0;j<ordnung+1; j++) {
rErg1[i][k] += dt * c1[j] * v[j][i][k];
vErg1[i][k] += dt * c1[j] * a[j][i][k];
rErg2[i][k] += dt * c2[j] * v[j][i][k];
vErg2[i][k] += dt * c2[j] * a[j][i][k];
}
}
}
if (swAenderbar) {
// entscheiden, ob die Schrittweite verändert wird.
double rErr=0;
double vErr=0;
// Fehlermaximum und Abstände suchen
for (i = 0; i < N;i ++) {
double localRErr=0;
double localVErr=0;
for (k = 0;k < 3;k ++) {
localRErr+=sq(rErg1[i][k]-rErg2[i][k]);
localVErr+=sq(vErg1[i][k]-vErg2[i][k]);
}
localRErr=sqrt(localRErr);
localVErr=sqrt(localVErr);
if (localRErr>rErr) rErr=localRErr;
if (localVErr>rErr) rErr=localVErr;
}
if (t<1.1*fabs(dt)) { // in der ersten Runde werden die Fehlergrenzen aus der SW berechnet
rErrMin=rErr/2;
rErrMax=rErr;
}
/* Schrittweitensteuerung à la piu58:
- es wird ein absoluter Fehler in r ausgewertet
- wenn der Fehler größer als die rErrMax ist,
dann wird die Schrittweite um 20% verringert
- wenn der Fehler kleiner ist als rErrMin,
dann wird die Schrittweite um 20% erhöht
*/
if (rErr<rErrMin) {
dt*=1.2;
// cerr<<dt<<" ^\n";
}
if (rErr>rErrMax) {
dt*=0.8;
// cerr<<dt<<" v\n";
}
} // Schrittweite darf geändert werden
//
for (i = 0; i < N;i ++) { // das genauere Ergebnis
for (k = 0;k < 3;k ++) {
r[0][i][k]=rErg2[i][k];
v[0][i][k]=vErg2[i][k];
}
}
// minimalen und maximalen Abstand suchen
for (i = 0; i < N;i ++) { // das genauere Ergebnis
double myR=sqrt( sq(r[0][i][0]) + sq(r[0][i][1]) + sq(r[0][i][2]));
if (rMin[i]>myR) rMin[i]=myR;
if (rMax[i]<myR) rMax[i]=myR;
}
if (progressTime > p) {
protokoll(I, t, p, dt); progressTime=0;
}
} // Zeitschritt
return sgn(dt)*t;
} // RKF45
void rk4(int tsteps, double dt) {
int i,k; // Butcher-Tableau
for (int t = 0;t < tsteps;t++) { // (tim) a0 a1 a2 a3
beschleunigung (N, m, r[0], a[0]); // %
for (i = 0; i < N;i ++) { // 0.5 0.5
for (k = 0;k < 3;k ++) {
r[1][i][k] = r[0][i][k] + 0.5 * dt * v[0][i][k];
v[1][i][k] = v[0][i][k] + 0.5 * dt * a[0][i][k];
}
}
beschleunigung (N, m, r[1], a[1]); //k1
for (i = 0; i < N;i ++) { // 0.5 0 0.5
for (k = 0;k < 3;k ++) {
r[2][i][k] = r[0][i][k] + 0.5 * dt * v[1][i][k];
v[2][i][k] = v[0][i][k] + 0.5 * dt * a[1][i][k];
}
}
beschleunigung (N, m, r[2], a[2]);
for (i = 0; i < N;i ++) {
for (k = 0;k < 3;k ++) {
r[3][i][k] = r[0][i][k] + 1 * v[2][i][k] * dt; // 1 0 0 1
v[3][i][k] = v[0][i][k] + 1 * a[2][i][k] * dt;
}
}
beschleunigung (N, m, r[3], a[3]);
for (i = 0; i < N;i ++) { // Ende 1/6 1/3 1/3 1/6
for (k = 0;k < 3;k ++) {
r[0][i][k] += (v[0][i][k] + 2 * v[1][i][k] + 2 * v[2][i][k] + v[3][i][k]) * dt / 6;
v[0][i][k] += (a[0][i][k] + 2 * a[1][i][k] + 2 * a[2][i][k] + a[3][i][k]) * dt / 6;
}
}
}
} // RK4
double init(string fname) {
ifstream mystream(fname);
if (!mystream.good()) error("Datei nicht zum lesen bereit");
string z;
getline(mystream, z); split(z.c_str());
N=toIntE(0); double startzeit=toDoubleE(1);
int i=0;
while (mystream.good()) {
getline(mystream, z); split(z.c_str());
m[i]=toDoubleE(0);
r[0][i][0]=toDoubleE(1);
r[0][i][1]=toDoubleE(2);
r[0][i][2]=toDoubleE(3);
v[0][i][0]=toDoubleE(4);
v[0][i][1]=toDoubleE(5);
v[0][i][2]=toDoubleE(6);
i++;
if (i==N) break;
}
if (i<N) error("zu wenige Zeilen");
return startzeit;
}
double T = -21544.5*sekProTag;
double p = 10 * sekProJahr; // alle 10a
double dt = 1 * sekProTag;
double t0 = 0;
int main(int argc, char **argv) {
// init1(); // Sonnensystem in astronomischen Koordinaten
init2(); // Sonnensystem in SI-Einheiten
double tEnd=rkf45(T, dt, p, I);
double restzeit=(T-tEnd);
tEnd=rkf45(0.99999*restzeit, restzeit, p, "99"); // ein Schritt zur Zeitkorrektur, 99: ohne Protokoll /**/
}
Registrieren
Anmelden
FAQ
Suche
