Aplikacja GRAWITACJA służy do przeprowadzania eksperymentów fizycznych na komputerze. Pozwala symulować ruch obiektów w polu grawitacyjnym oraz loty kosmiczne. Projekt powstał jako praca dyplomowa z fizyki, którą studiowałem na Wydziale Fizyki i Informatyki Stosowanej na krakowskiej AGH. 1500 wierszy kodu to niewiele ale stopień skomplikowania bardzo duży. Program kosztował dwa miesiące wakacji, po 10 godzin dziennie!
Inne uwagi:
- Program napisany przy użyciu VB i pracuje w środowisku MS Excela. Ułatwia to obliczenia i pozwala na szybkie modyfikowanie kodu lecz stwarza problemy w obszarze animacji.
- Program z założenia miał być swoistym rodzajem laboratorium elektronicznym do przeprowadzania eksperymentów i dlatego większą uwagę poświęcono liczbom.
- Możliwość symulowania czterech obiektów kosmicznych i rakiety – to raczej niewiele, ale celem miało być zasymulowanie lotu na Marsa.
- Tory ruchu planet wyliczane są z równań ruchu, które wywodzą się wprost z prawa powszechnego ciążenia Newtona.
- Do wyliczania współrzędnych na płaszczyźnie wykorzystano algorytmy: Eulera, Verleta i Runge-Kutta.
- Siła ciągu rakiety uwzględnia zmianę jej masy, co jest zgodne z wzorami wyprowadzonymi przez Konstantina Ciołkowskiego.
- W trakcie lotu możemy obracać rakietą wokół osi oraz sterować siłą ciągu.
- Wszelkie zmiany mogą być dokonywane w trakcie symulacji.
- 1500 wierszy kodu źródłowego, nie jest może liczbą imponującą lecz otwarta konstrukcja umożliwia różnorodne zastosowania.
Wszelkie szczegółowe informacje zostały opisane dokładnie w pracy dyplomowej, której fragmenty znajdują się w poniższych załącznikach:
mgr inż. Wacław Libront
Fragment kodu źródłowego:
'RAKIETA *********************************************************************
'wyliczanie chwilowego przyspieszenia pochodzącego od siły ciągu rakiety
'siła ciągu zawsze działa zgodnie z kierunkiem pochylenia rakiety
'który może być inny niż kierunek lotu wyliczony z grawitacji
Sub AA_rakieta()
Dim ax, ay, vx, vy, a, v, vx1, vy1, ac, acx, acy As Double
On Error GoTo błąd
'wartości chwilowe przyspieszeń i prędkości rakiety
ax = AA(0, 1)
ay = AA(0, 2)
vx = VV(0, 1)
vy = VV(0, 2)
'wypadkowe
a = (ax ^ 2 + ay ^ 2) ^ (1 / 2)
v = (vx ^ 2 + vy ^ 2) ^ (1 / 2)
'przyspieszenia po osiach rakiety wyliczamy z Talesa
'proporcjonalne do wektorów prędkości rakiety
'przyspieszenie zgodne ze zwrotem przepustnicy rakiety
'i uwzględniamy obrót automatyczny
'kąt = KĄT_GRANICA(KĄT_RAKIETY + KĄT_SKOKU_OBROTU)
'kąt = KĄT_GRANICA(KĄT_RAKIETY + KIERUNEK_OBROTU)
kąt = KĄT_RAKIETY
'wyliczone kierunki prędkości
vx1 = WERSOR_X(kąt)
vy1 = WERSOR_Y(kąt)
'nowe prędkości składowe po zmianie kąta
'po co to było liczone?
'vx = v * Abs(Cos(kąt * PI / 180)) * vx1
'vy = v * Abs(Sin(kąt * PI / 180)) * vy1
'przyspieszenie uwzględnia masę chwilową rakiety - Ciołkowski
'nie przeliczamy na inne jednostki, bo za szybko
'jest w N a my chcemy na kilogramy?
ac = SIŁA_CIĄGU / MASA_RAKIETY
'przyspieszenie ciągu działa w tą stronę, w którą jest zwrócona rakieta
acx = ac * Abs(Cos(kąt * PI / 180)) * vx1
acy = ac * Abs(Sin(kąt * PI / 180)) * vy1
'nowe składowe przyspieszeń rakiety
'po uwzględnieniu siły ciągu i kąta odchylenia
AA(0, 1) = AA(0, 1) + acx
AA(0, 2) = AA(0, 2) + acy
błąd:
End Sub
'obliczenie ilości paliwa i masy rakiety, gdy pracuje silnik
'i zaznaczono, że ma pobierać paliwo
Sub PALIWO_oblicz()
If PALIWO_POBIERAJ And SILNIK_PRACUJE Then
If MASA_PALIWA_AKT <= 0 Then
PALIWO_POBIERAJ = False
SILNIK_onoff 'wyłączenie silnika
Else
'powinniśmy uwzględniać skok czasu, ale bardzo szybko by spalało
'nie uwzględniamy czasu spalania dlatego bez * CZAS_SKOKU
MASA_PALIWA_AKT = MASA_PALIWA_AKT - SZYBKOŚĆ_SPALANIA * SKOK_CZASU
MASA_RAKIETY = MASA_RAKIETY - SZYBKOŚĆ_SPALANIA * SKOK_CZASU
MA(0) = MASA_RAKIETY 'zmniejszyć masę rakiety do obliczeń
End If
PALIWO_ustaw
End If
End Sub
'obliczanie kątów, położenia i kierunku rakiety
'gdy pracuje silnik
'i próby automatycznego wychodzenia z korkociągu,
'gdy prędkość spada do zera a przeciwny zwrot
Sub KĄTY_oblicz()
Dim vx, vy, vx1, vy1 As Variant
'położenie rakiety
vx = Sheets("GRAWITACJA").Range("E10")
vy = Sheets("GRAWITACJA").Range("F10")
'położenie planety środkowej
vx1 = Sheets("GRAWITACJA").Cells(10 + WYKRES_ŚRODEK_NR, 5)
vy1 = Sheets("GRAWITACJA").Cells(10 + WYKRES_ŚRODEK_NR, 6)
'nowe położenie rakiety
vx = vx - vx1
vy = vy - vy1
'v = (vx ^ 2 + vy ^ 2) ^ (1 / 2)
'wyliczenie kąta lotu
KĄT_LOTU = KĄT_oblicz(vx, vy)
'jeśli niema kompensacji to nie ma problemów z korkociągami
'bo sobie sam oblicza i reguluje
'tu wstawić czy kompensacja obrotów rakiety
If KOMPENSACJA Then
'prawidłowa reakcja na korkociąg gdy silnik pracuje
If Not (SILNIK_PRACUJE) Then
KĄT_RAKIETY = KĄT_GRANICA(KĄT_LOTU + ODCHYLENIE)
Else
ODCHYLENIE = KĄT_GRANICA(KĄT_RAKIETY - KĄT_LOTU)
End If
End If
'obraca się zgodnie z szybkością zmian kąta lotu
If OBRACAĆ = True Then
KĄT_RAKIETY = KĄT_GRANICA(KĄT_RAKIETY + KĄT_SKOKU_OBROTU * KIERUNEK_OBROTU)
ODCHYLENIE = KĄT_GRANICA(KĄT_RAKIETY - KĄT_LOTU)
Else
End If
End Sub
'ALGORYTMY *******************************************************************
'obliczanie odległości pomiędzy obiektami
Function R_oblicz(i, j As Variant) As Double
r2 = (XY(i, 1) - XY(j, 1)) ^ 2 + (XY(i, 2) - XY(j, 2)) ^ 2
R_oblicz = r2 ^ (1 / 2)
End Function
'obliczanie przyspieszeń z wzoru na grawitację
'a=-G*M/R^2
Sub A_oblicz()
On Error GoTo błąd
Dim R, R3 As Double
For i = 0 To WYM
AA(i, 1) = 0
AA(i, 2) = 0
For j = 0 To WYM
'sam ze sobą nie oddziaływuje i sprawdzamy tylko obiekty z masą >0
If i <> j And MA(i) > 0 And MA(j) > 0 Then
'odległości pomiędzy planetami
R = R_oblicz(i, j)
R3 = R ^ 3
'suma wszystkich składowych przyspieszeń daje przyspieszenie wypadkowe
For k = 1 To 2
'sumujemy przyspieszenie chwilowe
AA(i, k) = AA(i, k) + -GRAW * MA(j) * (XY(i, k) - XY(j, k)) / R3
Next k
End If
Next j
Next i
Exit Sub
błąd:
'gdy błąd, bo np promień równy 0 to omijamy i realizujemy dalej
Resume Next
End Sub
'RUNGE-KUTTY 4
'zapamiętaj początkowe
Sub XV0_oblicz()
For i = 0 To WYM
For k = 1 To 2
XY0(i, k) = XY(i, k)
VV0(i, k) = VV(i, k)
Next k
Next i
End Sub
'nowe prędkości z początkowym V i obliczonym K
Sub KXV1_oblicz()
For i = 0 To WYM
For k = 1 To 2
VV1(i, k) = AA(i, k) * SKOK_CZASU
XY(i, k) = XY0(i, k) + VV1(i, k) * SKOK_CZASU / 2
Next k
Next i
End Sub
Sub KXV2_oblicz()
For i = 0 To WYM
For k = 1 To 2
VV2(i, k) = AA(i, k) * SKOK_CZASU
XY(i, k) = XY0(i, k) + VV2(i, k) * SKOK_CZASU / 2
Next k
Next i
End Sub
Sub KXV3_oblicz()
For i = 0 To WYM
For k = 1 To 2
VV3(i, k) = AA(i, k) * SKOK_CZASU
XY(i, k) = XY0(i, k) + VV3(i, k) * SKOK_CZASU
Next k
Next i
End Sub
Sub K4_oblicz()
For i = 0 To WYM
For k = 1 To 2
VV4(i, k) = AA(i, k) * SKOK_CZASU
Next k
Next i
End Sub
'prędkość w metodzie RK
Sub VXRK_oblicz()
For i = 0 To WYM
For k = 1 To 2
VV(i, k) = VV0(i, k) + _
(VV1(i, k) + 2 * VV2(i, k) + 2 * VV3(i, k) + VV4(i, k)) / 6
XY(i, k) = XY0(i, k) + VV(i, k) * SKOK_CZASU
Next k
Next i
End Sub
'METODA Runge-Kutty4
Sub RK4()
'zapamiętaj poprzednie V i X
XV0_oblicz
'1 krok RK
A_oblicz
If (SILNIK_PRACUJE) And (SKOK_CZASU > 0) Then AA_rakieta
KXV1_oblicz
'2 krok RK
A_oblicz
If (SILNIK_PRACUJE) And (SKOK_CZASU > 0) Then AA_rakieta
KXV2_oblicz
'3 krok RK
A_oblicz
If (SILNIK_PRACUJE) And (SKOK_CZASU > 0) Then AA_rakieta
KXV3_oblicz
'4 krok RK
A_oblicz
'uwzględniamy silnik
If (SILNIK_PRACUJE) And (SKOK_CZASU > 0) Then AA_rakieta
K4_oblicz 'K4
'a teraz obliczenie V i X jako 1/6
VXRK_oblicz 'do tablic VV ze startowego V i 1/6(K1+2K2+2K3+k4)
'do tablic XY z przyrostem czasu dt
End Sub
