Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

Zespół Szkół Ogólnokształcących w Bobowej

graw01Aplikacja 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:

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
Zobacz tutaj