Forum der Vereinigung der Sternfreunde

Forum of the German Amateur Astronomy Association
Aktuelle Zeit: 22. Oktober 2020, 14:56:29 PM

Alle Zeiten sind UTC+02:00




Ein neues Thema erstellen  Auf das Thema antworten  [ 4 Beiträge ] 
Autor Nachricht
 Betreff des Beitrags: Heft 73/74 - Das Mehrkörperproblem
BeitragVerfasst: 04. April 2020, 21:00:28 PM 
Offline
Meister
Benutzeravatar

Registriert: 20. Januar 2013, 20:03:54 PM
Beiträge: 958
Wohnort: Leipzig
Heute habe ich das VdS-Journal 73 aus dem Kasten geholt. Ich beginne in diesem Heft eine zweiteilige Serie zum Mehrkörperproblem. Ab drei Körpern lässt sich die Bewegung im Gravitationsfeld nicht mehr geschlossen lösen. Die Serie erläutert, wie man theoretisch herangehen kann und gibt am Schluss (in Heft 74) ein sehr leistungsfähiges Verfahren an.

Die Programme zur Serie sind auf unserer Webseite hinterlegt. Da der gesamte Satz Startwerte hinterlegt ist, wird das Programm recht lang. Deshalb klebe ich es nicht hier hinein.

_________________
Uwe Pilz, Fachgruppen Kometen und Astrophysik/Algorithmen.
Oft benutzte Instrumente: Swarovski SLC 7x50B, Fujinon 16x70 FMT-SX-2, TMB Apo 105/650, Skywatcher Evostar 120/900 ED, Ninja Dobson 320/1440


Nach oben
   
BeitragVerfasst: 29. Juni 2020, 14:28:00 PM 
Offline
Meister
Benutzeravatar

Registriert: 20. Januar 2013, 20:03:54 PM
Beiträge: 958
Wohnort: Leipzig
Liebe Sternfreunde,

im Heft 74 habe ich nun eine hochwertige Route zur Berechnung des Mehrkörperproblems veröffentlicht. Es ist sehr viel Programmcode dabei, und wir haben kurz überlegt, ob wir das wirklich abdrucken sollen. Dieses Programm halte ich aber für so fundamental (und damit auch nach Jahrzehnten noch nutzbar), dass ich mich nicht allein auf die im Internet gespeicherten Version verlassen möchte.
Simuliert wird die Bewegung von mehreren, möglicherweise auch vielen Punktmassen. Damit kann man verschiedenen Fragestellungen nachgehen. Im Aufsatz selbst ist die Positionsberechnung im Sonnensystem angesprochen und eine Art "Störungsrechnung". Mir fallen sofort eine ganze Reihe von astronomischen Fragen, die man mit diesem Algorithmus beantworten kann:

- Man kann dieses Programm um weitere Körper ergänzen, z.B. Kometen oder Asteroiden und die sog. gestörte Bahn berechnen.
- man kann nach zyklischen (periodischen) Bahnen im Sonnensystem suchen, wie sie bei einigen Asteroiden auftauchen
- Es ist möglich, die Stabilität von Mehrkörpersystemen zu untersuchen, z.B. Mehrfachsterne
- Die dynamische Entwicklung von Stern- und Galaxienhaufen lässt sich simulieren

Bei dieser Fülle von Anwendungen findet man natürlich verschiedene Ansätze und auch ausführbare Programme anderer Autoren:

- Schon in den 80er Jahren habe ich mit dem sog. Heidelberger Programm von Joachim Schubart experimentiert. Es verfolgt ähnliche Ziele wie das abgedruckte Programm. Es ist jedoch ein Mehrschrittverfahren, welches mehrere zeitlich zurückliegende Ergebnisse erfordert. Diese sind am Anfang der Rechnung nicht vorhanden, so dass eine zusätzliche Anlaufrechnung benötigt wird.
Die Genauigkeit ist sehr hoch, aber die Schrittweitensteuerung ist schwierig, weil stets eine neue Anlaufrechnung nötig ist. Das Programm enthält Startwerde des Planetensystems für die Epoche 2'430'000,5 (6.1.1941).


- Im Buch "Astronomie mit dem Personalcomputer" von Oliver Montenbruck und Thomas Pfleger ist das Programm DE enthalten (in C++). DE heißt differential equations und ist die C++-Portierung des frei verfügbaren Fortran-Programms ODE von Lawrence Shampine und Marilyn Gordon. So wie das Programm im Buch abgedruckt ist, berechnet es einen Kleinkörper in Gravitationsfeld der Planeten. Die Planeten werden aber nicht wie bei mir numerisch simuliert, sondern mit der normalen Planetentheorie. Man kann dies aber auf ein Mehrkörperprogramm umbauen.
Das Programm verwendet ebenfalls ein Mehrschrittverfahren und ist sehr genau. Die Anlaufrechnung ist integriert, so dass die Schrittweitensteuerung einfach vorgenommen werden kann. Allerdings ist der Programmcode recht komplex.


- Sverre Aarseth hat fast sein ganzes Berufsleben der Mehrkörpersimulation gewidmet. Er zielte insbesondere auf Viel-Körper-Systeme ab, wie Offene Sternhaufen oder gar Kugelsternhaufen. Hierbei steht die exakte Position der einzelnen Körper nicht im Vordergrund, sondern statistische Größen, wie Änderung der Sterndichte über dem Radius in Abhängigkeit der Zeit, Rate an Doppel- und Mehrfachsternen und die Gesamtlebensdauer.
Das Programm ist ziemlich aufwendig und nach so vielen Jahrzehnten sehr ausgereift. Es rechnet stets in mehreren Zeitebenen, damit die Rechengeschwindigkeit nicht durch enge Systeme gebremst wird. Das Programm lässt sich mit dem Gnu-Fortran übersetzen.
Sverre Aarseth hat sogar ein Buch über seine Programme publiziert, es steht bei mir im Regal.


Ich habe versucht, ein inhaltlich einigermaßen überschaubares, kurzes Programm zu schreiben, was dennoch leistungsfähig ist. Es gehört zur Klasse der Runge-Kutta-Verfahren und verwendet keine zurückliegenden Werte. Damit ist eine Schrittweitensteuerung einfach und wurde in der Variante von Erwin Fehlberg berücksichtigt. Die Genauigkeit entspricht (bezogen auf die Rechengeschwindigkeit) durchaus den Ansätzen von Schubart und Shampine/Gordon.
Ich habe Experimente mit Vielkörpersystemen gemacht. Doppel- und Mehrfachsterne muss man dann berücksichtigen, was ich auf einfache Weise getan habe. Dennoch reicht man mit diesem recht einfachen Programm nicht an die Qualitäten von Sverre Aarseth heran. Die Rechenlast bei Vielkörpersystemen ist sehr hoch, so dass sich eine Unetrstützung von Mehrprozessormaschinen oder Grafikkarten lohnt. Aarseth oder seine Nachfolger haben hierfür Lösungen geschrieben, welche alle öffentlich nutzbar sind.

Ich werde künftig auf dieses Programm zurückgreifen, um weitere fesselnde Probleme zu berechnen.

_________________
Uwe Pilz, Fachgruppen Kometen und Astrophysik/Algorithmen.
Oft benutzte Instrumente: Swarovski SLC 7x50B, Fujinon 16x70 FMT-SX-2, TMB Apo 105/650, Skywatcher Evostar 120/900 ED, Ninja Dobson 320/1440


Nach oben
   
BeitragVerfasst: 29. Juni 2020, 14:30:17 PM 
Offline
Meister
Benutzeravatar

Registriert: 20. Januar 2013, 20:03:54 PM
Beiträge: 958
Wohnort: Leipzig
Das Python-Programm.
Code:
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?")


_________________
Uwe Pilz, Fachgruppen Kometen und Astrophysik/Algorithmen.
Oft benutzte Instrumente: Swarovski SLC 7x50B, Fujinon 16x70 FMT-SX-2, TMB Apo 105/650, Skywatcher Evostar 120/900 ED, Ninja Dobson 320/1440


Nach oben
   
BeitragVerfasst: 29. Juni 2020, 14:50:47 PM 
Offline
Meister
Benutzeravatar

Registriert: 20. Januar 2013, 20:03:54 PM
Beiträge: 958
Wohnort: Leipzig
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 /**/
}

_________________
Uwe Pilz, Fachgruppen Kometen und Astrophysik/Algorithmen.
Oft benutzte Instrumente: Swarovski SLC 7x50B, Fujinon 16x70 FMT-SX-2, TMB Apo 105/650, Skywatcher Evostar 120/900 ED, Ninja Dobson 320/1440


Nach oben
   
Beiträge der letzten Zeit anzeigen:  Sortiere nach  
Ein neues Thema erstellen  Auf das Thema antworten  [ 4 Beiträge ] 

Alle Zeiten sind UTC+02:00


Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 1 Gast


Du darfst keine neuen Themen in diesem Forum erstellen.
Du darfst keine Antworten zu Themen in diesem Forum erstellen.
Du darfst deine Beiträge in diesem Forum nicht ändern.
Du darfst deine Beiträge in diesem Forum nicht löschen.
Du darfst keine Dateianhänge in diesem Forum erstellen.

Suche nach:
Gehe zu:  
cron
Powered by phpBB® Forum Software © phpBB Limited
Deutsche Übersetzung durch phpBB.de