Ich untersuche einen bestimmten Fall des eingeschränkten Drei-Körper-Problems. Es wurde festgestellt, dass einige Objekte einem Hufeisenbahnmuster folgen, und ich versuche, etwas durch einen Integrationscode in C zu sortieren. Ich folge einigen Ratschlägen im Artikel Familien periodischer Hufeisenbahnen im eingeschränkten Drei-Körper-Problem . Das gibt mir ideale Anfangsbedingungen und die Gleichungen im Schwerpunktsystem. (m ist die Masse der Erde und die sich daraus ergebende Position der Sonne im Massenschwerpunkt-Bezugssystem (x, y) sind die Koordinaten des dritten Körpers, die als masselos angenommen werden (wie es das eingeschränkte Problem erfordert).
Die Positionen von "Sonne" und "Erde" sind im gleichen Bezugssystem auf (m, 0) und (m-1,0) festgelegt. (rotierendes Bezugssystem, vorausgesetzt, die Erde hat eine kreisförmige Umlaufbahn.)
Aus all dem habe ich die Gleichungen berechnet, um das System zu beschreiben:
Ich habe den Algorithmus von Runge-Kutta 4 verwendet, um diese Gleichungen zu integrieren. (Ich weiß, dass der Code ziemlich umwerfend ist, aber ich kann einfach keine Zeiger verwenden und ich verwende überall Strukturen).
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define dt 0.0001
#define N 100
typedef struct{
long double x,y;
}vec;
typedef struct{
vec k1,k2,k3,k4;
}runge;
typedef struct{
runge r,v;
}big;
double dS,dE,m;
double accx(double,double,double);
double accy(double,double,double);
void rad(vec);
big rungekutta(vec,vec);
vec moto(vec,runge);
double jacobi(vec);
int main(){
vec r,v;
big f;
double J,t;
int i,Num;
FILE* s1;
s1=fopen("HorseShoe.dat","w");
Num=(int)N/dt;
scanf("%Lf",&r.x);
scanf("%Lf",&r.y);
scanf("%Lf",&v.x);
scanf("%Lf",&v.y);
scanf("%lf",&m);
for(i=0;i<Num;i++){
t=(i+1)*dt;
rad(r);
f=rungekutta(r,v);
r=moto(r,f.r);
v=moto(v,f.v);
J=jacobi(r);
fprintf(s1,"%lf\t%Lf\t%Lf\t%Lf\t%Lf\t%lf\n",t,r.x,r.y,v.x,v.y,J);
}
return 0;
}
void rad(vec r){
dS=pow(r.x-m,2)+pow(r.y,2);
dE=pow(r.x-m+1,2)+pow(r.y,2);
}
double jacobi(vec r){
return pow(r.x,2)+pow(r.y,2)+2*(1-m)/dS+2*m/dE+m*(1-m);
}
double accx(double x,double y,double v){
return x-(x-m)*(1-m)/pow(pow(x-m,2)+pow(y,2),1.5)-m*(x-m+1)/pow(pow(x-m+1,2)+pow(y,2),1.5)+2*v;
}
double accy(double x,double y,double v){
return y-(1-m)*y/pow(pow(y,2)+pow(x-m,2),1.5)-m*y/pow(pow(y,2)+pow(x-m+1,2),1.5)-2*v;
}
big rungekutta(vec r,vec v){
big f;
f.r.k1.x=v.x;
f.r.k1.y=v.y;
f.v.k1.x=accx(r.x,r.y,v.y);
f.v.k1.y=accy(r.x,r.y,v.x);
f.r.k2.x=v.x+f.v.k1.x*dt/2;
f.r.k2.y=v.y+f.v.k1.y*dt/2;
f.v.k2.x=accx(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.y+f.v.k1.y*dt/2);
f.v.k2.y=accy(r.x+f.r.k1.x*dt/2,r.y+f.r.k1.y*dt/2,v.x+f.v.k1.x*dt/2);
f.r.k3.x=v.x+f.v.k2.x*dt/2;
f.r.k3.y=v.y+f.v.k2.y*dt/2;
f.v.k3.x=accx(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.y+f.v.k2.y*dt/2);
f.v.k3.y=accy(r.x+f.r.k2.x*dt/2,r.y+f.r.k2.y*dt/2,v.x+f.v.k2.x*dt/2);
f.r.k4.x=v.x+f.v.k3.x*dt;
f.r.k4.y=v.y+f.v.k3.y*dt;
f.v.k4.x=accx(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.y+f.v.k3.y*dt);
f.v.k4.y=accy(r.x+f.r.k3.x*dt,r.y+f.r.k3.y*dt,v.x+f.v.k3.x*dt);
return f;
}
vec moto(vec r,runge rk){
r.x+=(rk.k1.x+2*rk.k2.x+2*rk.k3.x+rk.k4.x)*dt/6;
r.y+=(rk.k1.y+2*rk.k2.y+2*rk.k3.y+rk.k4.y)*dt/6;
return r;
}
Beim Zeichnen der Ergebnisse erhalte ich nur eine Spirale, während ich bei Verwendung der angegebenen Eingaben eine Hufeisenbahn erhalten sollte. Ich habe viele verschiedene Eingaben versucht (m = 0,0001 und m = 0,000003, letzteres entspricht den tatsächlichen Werten der Erd- und Sonnenmasse (Sonnenmasse ist 1 m)).
Ich kann einfach nicht sehen, was los ist (wahrscheinlich alles: D). Kann mir bitte jemand helfen?
quelle
Antworten:
Ich werde weder Ihre Gleichungen noch Ihren Code sorgfältig durchgehen, aber da Sie die genauen Startpositionen, die Sie verwendet haben, nicht erwähnt haben, ist meine Vermutung, dass Sie nicht zuerst nach einem stabilen Hufeisen-Orbit- Zustandsvektor gesucht haben.
Die Positionen für Sonne und Erde sind im synodischen (rotierenden) Rahmen festgelegt. Aber wo hast du dein Astroid angefangen? Für eine bestimmte Startposition müssen Sie dem Asteroiden auch die richtige Startgeschwindigkeit (Geschwindigkeit und Richtung) geben, da er sonst möglicherweise nicht stabil ist und einfach abwandert oder sich wie erwähnt dreht.
Schauen Sie sich diese Antwort an, in der ich die richtigen Gleichungen zeige und dann nach stabilen Hufeisenbahnen suche, bevor ich sie zeichne.
Ich beginne mit einer Position auf der x-Achse gegenüber der Erde und etwas weiter oder näher von der Sonne entfernt als die Erde. Dann starte ich es mit einer geringen Geschwindigkeit (im rotierenden Rahmen) genau in Richtung der y-Achse. Ich beobachte, wie es zur Erde driftet und dann zurück in die Region des Startgebiets bumerang. Wenn es die x-Achse wieder kreuzt, überprüfe ich, wie nahe die Geschwindigkeit an y liegt. Ich benutze einen Optimierungs-Nulllöser
brentq, um die Anfangsgeschwindigkeit anzupassen, bis die Rücklaufgeschwindigkeit genau in die entgegengesetzte Richtung ist.Wenn dies zutrifft, ist die Umlaufbahn periodisch und wiederholt sich, und wir können sie unter den Bedingungen des CR3BP als veraltete Umlaufbahn bezeichnen.
Aus dieser Antwort (es gibt dort viel mehr zu lesen, einschließlich mehrerer Referenzen!):
quelle