Ruch ciała w polu sił centralnych
pochodzących od dwóch centrów nieruchomych

Autor: Mateusz Wojtaszek

Student IV roku Fizyki Fazy Skondensowanej na Wydziale Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego

Data utworzenia: Marzec 2009

Króki opis programu:
Program przedstawia zagadnienie oddziaływania trzech ciał na płaszczyźnie. Rozpatrywany jest przypadek ruchu ciała w polu sił centralnych pochodzących od dwóch centrów nieruchomych.

Dzięki opcji manipulate użytkownik może:
- rysować trajektorię ruchu ciała w czterech zakresach czasu: 10, 25, 50 i 100;
- rysować wektor prędkości początkowej, oraz ustawiać jego wartość (można wybrać opcję bez rysowania wektora);
- rysować przybliżony rozkład potencjału pochodzący od dwóch nieruchomych centrów (można wybrać opcję bez rozkładu potencjału);
- zmieniać masy wszystkich trzech ciał;
- zmieniać położenia wszystkich trzech ciał;
- resetować wyniki i powracać do ustawień początkowych.

Ciekawe przypadki:
Ciekawe przypadki są możliwe do zobaczenia po kliknieciu na odpowiedni przycisk w menu manipulate. Poniżej znajduje się krótki opis niektórych przykładów:
1) Obszar I - w tym obszarze wszystkie trajektorie znajdują się pomiędzy dwiema elipsami, w których ogniskach znajdują się spoczywające ciała (trajektorie okrążają obie masy centrów sił);
2) Obszar II - w tym obszarze wszystkie trajektorie zawarte są w przestrzeni ograniczonej przez elipsę, w której ogniskach znajdują się spoczywające ciała (trajektorie okrążają oba centra sił tworząc krzywe o kształcie ósemki);
3) Obszar III - w tym przypadku trajektorie obiegają jedno z centrów sił. Są one zawarte pomiedzy elipsą a hiperbolą. W ogniskach elipsy znajdują się centra sił.
4) Obszar IV - ruch, w którym obieg następuje jedynie wokół cięższego z dwóch ciał.

Zachęcam do szukania własnych ciekawych przykładów!

Bibliografia:
"Mechanik Klasyczna" - Grzegorz Białkowski

Clear["Global`*"];(*Zapobiega powtórzeniu sie takich samych znaków wcześniej wczytanych*)
trajektoria[
{{x10_,y10_},{x20_,y20_},{x30_,y30_}},
{{vx10_,vy10_},{vx20_,vy20_},{vx30_,vy30_}},
{m1_, m2_, m3_},
T_,Vec_,Pot_,Podp_,Planet_] := 
Module[{rownania,Tmax,funkcja},
rownania=NDSolve[
{x1'[t]==vx1[t],y1'[t]==vy1[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x20))/((((((x1[t]-x20))^2)+(((y1[t]-y20))^2)))^(3/2)))-((m1 m3(x1[t]-x30))/((((((x1[t]-x30))^2)+(((y1[t]-y30))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y20))/((((((x1[t]-x20))^2)+(((y1[t]-y20))^2)))^(3/2)))-((m1 m3(y1[t]-y30))/((((((x1[t]-x30))^2)+(((y1[t]-y30))^2)))^(3/2))),
x1[0]==x10,
y1[0]==y10,
vx1[0]==vx10,
vy1[0]==vy10},
{x1,y1,vx1,vy1},{t,0,T}];(*Rozwiązanie numeryczne*)
Tmax = rownania[[1, 1, 2, 1, 1, 2]];(*Wyłuskanie czasu maksymalnego*)
funkcja={{x1[t],y1[t]}/.rownania[[1]],{x20,y20},{x30,y30}};(*Zapisanie danych do wykresu z wyniku obliczeń numerycznych*)

(*Wykres trajektorii*)
ruch=ParametricPlot[Evaluate[funkcja](*Rysowanie wykresu*),
{t,0,Tmax},(*W takim zakresie*)
PlotStyle->{Red}(*KOlor trajektorii ruchu*)];

(*Strzałka prędkości początkowej*)
If[Vec,{strzalka=Graphics[{Black,Arrowheads[Medium],Arrow[{{x10,y10},{x10+vx10,y10+vy10}}]}]},{strzalka=Graphics[]}];

(*Wykres potencjału*)
If[Pot,{potencjal=ContourPlot[m2/Norm[{x,y}-{x20,y20}]+m3/Norm[{x,y}-{x30,y30}],{x,Min[x20,x30]-Max[2,0.5*(Abs[x20-x30])],Max[x20,x30]+Max[2,0.5*(Abs[x20-x30])]},{y,Min[y20,y30]-Max[2,Abs[x20-x30]],Max[y20,y30]+Max[2,Abs[x20-x30]]},Contours->10]},{potencjal=Graphics[]}];

If[Podp,{podpis=Graphics[{Text[Style["(Subscript[M, 2])",Black,Medium,Bold],{x20,y20+0.2}],Text[Style["(Subscript[M, 3])",Black,Medium,Bold],{x30,y30+0.2}],Text[Style["(Subscript[m, 1])",Medium,Bold],{x10,y10+0.2}]}]},{podpis=Graphics[]}];

If[Planet,{planety=Graphics[{{Blue,Disk[{x20,y20},2*m2/1000]},{Blue,Disk[{x30,y30},2*m3/1000]},{Red,Disk[{x10,y10},2*m1/1000]}}]},{planety=Graphics[]}];

(*Pokazanie wszystkich wykresów*)
Show[potencjal,planety,ruch,strzalka,podpis,PlotRange->All,ImageSize->700,AspectRatio -> Automatic,Frame->False,Axes -> False](*Pokazanie wykresów*)
]//Quiet
Manipulate[
trajektoria[{xy10,xy20,xy30},
{{pred1,pred2},{0,0},{0,0}},
{masa1,masa2,masa3}, T,Vec,Pot,Podp,Planet] ,
{{T,0.001,"Czas"},0.001, CzasMax},(*Zakres czasu*)

{CzasMax,{10,25,50,100}},(*Zwiększenie zakresu czasu*)
Delimiter,
"Opcje dotyczące animacji: ",
"Wektor prędkości początkowej:",
{Vec,{True,False}},

"Rozkład potencjału:",
{Pot,{True,False}},

"Podpisy mas:",
{Podp,{True,False}},

"Ukazanie mas:",
{Planet,{True,False}},
Delimiter,
(*Zdefiniowanie mas*)
{{masa1,0.01,"(Subscript[m, 1])"},0.001,100},
{{masa2,0.01,"(Subscript[M, 2])"},0.001,100},
{{masa3,0.01,"(Subscript[M, 3])"},0.001,100},
Delimiter,
(*Zdefiniowanie pedkości ciała pierwszego*)
{{pred1,0,"Vx"},-2,2},
{{pred2,0,"Vy"},-2,2},
Delimiter,
(*Położenia cząstek*)
{{xy10,{ 0,0}},Locator},
{{xy20,{1,1}},Locator},
{{xy30,{-1,1}},Locator},
Delimiter,
(*Przyciski pomagające w obsłudze*)
Button["Resetuj do początkowych wartości",{pred1=0,pred2=0,masa1=0.001,masa2=0.001,masa3=0.001,T=0.001}],
Button["Resetruj położenia",{xy10={0,0},xy20={1,1},xy30={-1,1}}],
Delimiter,
"",
Button["Obszar I",{pred1=2,pred2=0,T=0.1,masa1=0.001,masa2=0.8,masa3=4.5,xy10={0,0},xy20={1,1},xy30={-1,1}}],
Button["Obszar II",{pred1=1.9,pred2=-2,T=0.1,masa1=9.2,masa2=40,masa3=40,xy10={0.1,1.2},xy20={1,1},xy30={-1,1}}],
Button["Obszar III",{pred1=0,pred2=0,T=0.1,masa1=1,masa2=1,masa3=1,xy10={1.2,0.3},xy20={1,1},xy30={-1,1}}],
Button["Obszar IV",{pred1=0,pred2=0.3,T=0.1,masa1=1,masa2=5,masa3=0.001,xy10={-2,0},xy20={1,1},xy30={-1,1}}],

ControlPlacement->{Left}
]