# Programm RKF5.py
from turtle import *
from random import *
from math import *

# globale Variable
G = 6.67408e-11         #Gravitationskonstante
mProAE=149597870700    #Astronomische Einheit
sekProTag=86400
# sekProJahr=365.2425*86400
# mProLj=9460730472580800

N=10 # Anzahl der Körper
r0 = [[0 for i in range(3)] for j in range(10)]
v0 = [[0 for i in range(3)] for j in range(10)]
a0 = [[0 for i in range(3)] for j in range(10)]
r1 = [[0 for i in range(3)] for j in range(10)]
v1 = [[0 for i in range(3)] for j in range(10)]
a1 = [[0 for i in range(3)] for j in range(10)]
r2 = [[0 for i in range(3)] for j in range(10)]
v2 = [[0 for i in range(3)] for j in range(10)]
a2 = [[0 for i in range(3)] for j in range(10)]
r3 = [[0 for i in range(3)] for j in range(10)]
v3 = [[0 for i in range(3)] for j in range(10)]
a3 = [[0 for i in range(3)] for j in range(10)]
r4 = [[0 for i in range(3)] for j in range(10)]
v4 = [[0 for i in range(3)] for j in range(10)]
a4 = [[0 for i in range(3)] for j in range(10)]
r5 = [[0 for i in range(3)] for j in range(10)]
v5 = [[0 for i in range(3)] for j in range(10)]
a5 = [[0 for i in range(3)] for j in range(10)]
r0_4 = [[0 for i in range(3)] for j in range(10)] # Ergebnis für RKF 4. Ordnung
m = [0 for i in range (10)]

def startwerte():
    # Orte und Geschwindigkeiten
    r0[0][0]=   -1063436181;    v0[0][0]=      9.29323348
    r0[0][1]=    -403178315;    v0[0][1]=    -12.81658147
    r0[0][2]=      27924311;    v0[0][2]=     -0.15918117
    r0[1][0]=  -20525703382;    v0[1][0]=  37004.00242937
    r0[1][1]=  -67316663267;    v0[1][1]= -11177.58979901
    r0[1][2]=   -3651721547;    v0[1][2]=  -4307.67352332
    r0[2][0]= -108522074771;    v0[2][0]=   1392.28564030
    r0[2][1]=   -5295518213;    v0[2][1]= -35152.38582183
    r0[2][2]=    6163650722;    v0[2][2]=   -560.19899547
    r0[3][0]=  -27570384049;    v0[3][0]= -29777.08241127
    r0[3][1]=  144289541085;    v0[3][1]=  -5492.11604817
    r0[3][2]=      27924311;    v0[3][2]=     -0.15918117
    r0[4][0]=  206972399720;    v0[4][0]=   1173.70663100
    r0[4][1]=   -2404931177;    v0[4][1]=  26284.74809179
    r0[4][2]=   -5126803613;    v0[4][2]=    522.10695479
    r0[5][0]=  597504078214;    v0[5][0]=  -7900.54555615
    r0[5][1]=  439201609603;    v0[5][1]=  11143.31683419
    r0[5][2]=  -15198533263;    v0[5][2]=    130.70502469
    r0[6][0]=  957322522683;    v0[6][0]=  -7422.72126397
    r0[6][1]=  982453470791;    v0[6][1]=   6723.09678597
    r0[6][2]=  -55183935225;    v0[6][2]=    178.08847902
    r0[7][0]= 2157914958773;    v0[7][0]=   4646.31018375
    r0[7][1]=-2055025724501;    v0[7][1]=   4614.84545075
    r0[7][2]=  -35599134212;    v0[7][2]=    -43.05247532
    r0[8][0]= 2513984445571;    v0[8][0]=   4475.19788359
    r0[8][1]=-3739120513758;    v0[8][1]=   3063.80412436
    r0[8][2]=   19061170936;    v0[8][2]=   -166.22362888
    r0[9][0]=-1478397309067;    v0[9][0]=   5253.44576775
    r0[9][1]=-4182980594505;    v0[9][1]=  -2675.42458857
    r0[9][2]=  875240628161;    v0[9][2]=  -1233.31598335
    # 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

def sq(x):
    return x*x;

def abstand(r, p, q):
    	return sqrt( sq(r[p][0]-r[q][0]) + sq(r[p][1]-r[q][1]) + sq(r[p][2]-r[q][2]) )
    
def beschleunigung(N, m, r, a):
    for p in range(N):    # alle Beschleunigungen löschen
        for k in range (3):
            a[p][k]=0;
    for p in range(N): # alle Planeten
        for q in range(p+1, N): # gegen alle anderen Planeten
            R=abstand(r, p, q)
            R3=R*R*R;
            for k in range(3):
                A=G*(r[q][k]-r[p][k])/R3
                a[p][k]+=m[q]*A
                a[q][k]-=m[p]*A

def rkf5(tsteps, dt):
    for t in range(tsteps):
        beschleunigung(N, m, r0, a0);						
        for i in range(N):
            for k in range(3):
                r1[i][k] = r0[i][k] + \
                    1/4*dt * v0[i][k]	
                v1[i][k] = v0[i][k] + \
                    1/4*dt * a0[i][k]
        beschleunigung(N, m, r1, a1);
        for i in range(N):
            for k in range(3):
                r2[i][k] = r0[i][k] + \
                    3/32 * dt*v0[i][k] + \
                    9/32 * dt*v1[i][k]
                v2[i][k] = v0[i][k] + \
                    3/32 * dt*a0[i][k] + \
                    9/32 * dt*a1[i][k]
        beschleunigung(N, m, r2, a2);						
        for i in range(N):
            for k in range(3):
                r3[i][k] = r0[i][k] + \
                    1932/2197 * dt*v0[i][k] + \
                    -7200/2197 * dt*v1[i][k] + \
                    7296/2197 * dt*v2[i][k] 
                v3[i][k] = v0[i][k] + \
                    1932/2197 * dt*a0[i][k] + \
                    -7200/2197 * dt*a1[i][k] + \
                    7296/2197 *  dt*a2[i][k]
        beschleunigung(N, m, r3, a3);						
        for i in range(N):
            for k in range(3):
                r4[i][k] = r0[i][k] + \
                    439/216 *   dt*v0[i][k] + \
                    -8 *        dt*v1[i][k] + \
                    3680/513  * dt*v2[i][k] + \
                    -845/4104 * dt*v3[i][k]
                v4[i][k] = v0[i][k] + \
                    439/216 *   dt*a0[i][k] + \
                    -8 *        dt*a1[i][k] + \
                    3680/513 *  dt*a2[i][k] + \
                    -845/4104 * dt*a3[i][k] 
        beschleunigung(N, m, r4, a4);						
        for i in range(N):
            for k in range(3):
                r5[i][k] = r0[i][k] + \
                    -8/27 *      dt*v0[i][k] + \
                    2 *          dt*v1[i][k] + \
                    -3544/2565 * dt*v2[i][k] + \
                    1859/4104 *  dt*v3[i][k] + \
                    -11/40 *     dt*v4[i][k] 
                v5[i][k] = v0[i][k] + \
                    -8/27 *      dt*a0[i][k] + \
                    2 *          dt*a1[i][k] + \
                    -3544/2565 * dt*a2[i][k] + \
                    1859/4104 *  dt*a3[i][k] + \
                    -11/40 *     dt*a4[i][k] 
        beschleunigung(N, m, r5, a5);						             
        for i in range(N):
            for k in range(3):
                r0_4[i][k] = r0[i][k]
                # Ergebnisse 4. Ordnung zur Fehlerberechnung
                r0_4[i][k]+= dt* ( \
                    25/216 *       v0[i][k] + \
                    0 *            v1[i][k] + \
                    1408/2565 *    v2[i][k] + \
                    2197.0/4104 *  v3[i][k] + \
                    -1/5 *         v4[i][k] + \
                    0 *            v5[i][k] \
                )
                r0[i][k]+= dt* ( \
                    16/135 *      v0[i][k] + \
                    0*            v1[i][k] + \
                    6656/12825 *  v2[i][k] + \
                    28561/56430 * v3[i][k] + \
                    -9/50 *       v4[i][k] + \
                    2/55 *        v5[i][k] \
                )
                v0[i][k]+= dt* ( \
                    16/135 *      a0[i][k] + \
                    0*            a1[i][k] + \
                    6656/12825 *  a2[i][k] + \
                    28561/56430 * a3[i][k] + \
                    -9/50 *       a4[i][k] + \
                    2/55 *        a5[i][k] \
                )
        # Fehler ausrechnen. Damit kann man nur etwas anfangen,
        # wenn man die for-Schleife durch eine while-Schleife ersetzt
        err=0;
        for i in range(N): # welcher Körper gibt den größten Fehler?
            localErr=0; 
            for k in range(3):
                localErr+=sq(r0_4[i][k]-r0[i][k])
            localErr=sqrt(localErr)
            if (localErr>err):
                err=localErr;



rMin = [0 for i in range(N)]
rMax = [0 for i in range(N)]

# Hauptprogramm
startwerte()
for i in range(2000):
    # minimale und maximalen Sonnenabstand löschen
    for j in range(N):
        rMin[j]=7e77
        rMax[j]=0
    # je 100 Jahre rechnen
    for j in range(36524):     # 36524 Schleifen ...
        rkf5(2, 0.5*sekProTag) # .. zu 1 Tagen = 36524 Tage 
    # minimalen und maximalen Sonnenabstand ermitteln
        for j in range(N):
            r=sqrt( sq(r0[j][0]) + sq(r0[j][1]) + sq(r0[j][2]))
            if (rMin[j]>r):
                rMin[j]=r;
            if (rMax[j]<r):
                rMax[j]=r;
    # minimalen und maximalen Sonnenabstand ausgeben
    print("Minimum, t=",i/10," : ",end='')
    for j in range(N):
        print (rMin[j]/mProAE," ",end='')
    print("")
    print("Maximum, t=",i/10," : ",end='')
    for j in range(N):
        print (rMax[j]/mProAE," ",end='')
    print("")

name = input("Fertig?")

/*
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 /**/
}