Heittoliike_drag_verlet.py

# Heittoliikkeen simulointia,  mukana ilmanvastus kv²
# k riippuu ilman tiheydestä, joka riippuvuu korkeudesta ja lämpötilasta
# lämpötila riippuvuu korkeudesta
# joten k:n arvo muuttuu korkeuden mukaan
# a = dv/dt = -g - kv²
# verlet second order method
# Juhani Kaukoranta 23.8.2021
import pygame, math
# pygame:ssa x kasvaa vasemmalta oikealle, y ylhäältä alas
pygame.font.init()
font = pygame.font.Font('freesansbold.ttf', 24)
# pallon ominaisuudet: (otettu baseballista, Cd=0.3)
# ei ilmanvastusta Cd = 0, paraabelirata
Cd = 0.3 # ilmanvastuskerroin, oletusarvo 0.3
roo = 1.225 # ilman tiheys kg/m³ (ICAO, 0m +15C)
A = 0.0043 # pallon poikkipinta-ala m²
m = 0.145 # pallon massa kg
r = 5 # piirretyn pallon säde pikseleinä
# ilman tiheys roo = 1.225*(288.15/(288.15-0.0065*h))**(-4.25578597)
k1 = (Cd*A)/(2*m) # pallon vastuskertoimen osa
# k = k1*roo # pallon vastuskerroin, mukana ilman tiheys roo
v0 = 300 # lähtönopeus
g = 9.80665 # putouskiihtyvyys
alfadeg = 30 # heittokulma (suhteessa vaakatasoon)
alfa = alfadeg*math.pi/180 # heittosuunta radiaaneina
# pystykiihtyvyys y''(t) = -g - k*vy² = -(g+k*vy²)
# vaakakiihtyvyys x''(t) = -k *vx²
t0 = 0.0 # aloitusaika
Dt = 1/100.0 # aika-askel
T = 7000 # animaation kesto 7000*1/100 = 70 s
# piirtoalueen asetus
(width,height) = (1280,720) # Pygamen piirtoala pikseleinä
screen = pygame.display.set_mode((width, height))
black = (0,0,0)
yellow =(255,255,0)
red = (255,0,0)
green = (0,255,0)
blue =(0,0,255)
magenta = (255,0,255)
screen.fill(black)
pygame.display.set_caption("pallon lento, ilmanvastus kv² ja tiheyden korkeusriippuvuus mukana: ")
# pallon (0,0) vastaa  Pygamen näytöllä (100,500)
(x0,y0) =(0,0) # pallon koordinaattien aloitusarvot, (0 metriä, 0 metriä)
(pygamex0,pygamey0) = (100,500) # Pygamen koordinaatiston origo
h = 0 # korkeus lähtötasosta tasosta,lähtötaso on merenpinnan taso, Pygamen y=500=gamey0
roo = 1.225*(288.15/(288.15-0.0065*h))**(-4.25578597) # ilman tiheyden alkuarvo
k = k1*roo # pallon ilmanvastuskertoimen aloitusarvo
vx = v0*math.cos(alfa) # vaakasuoran nopeuskomponentin alkuarvo
vy = v0*math.sin(alfa) # pystysuoran nopeuskomponentin alkuarvo
(x1,y1) =(x0 + vx*Dt,y0 + vy*Dt) # eka askel Eulerin avulla
hmax = y1 # pitää kirjaa pallon maksimikorkeudesta
h = y1 # pallon korkeus
timer = 0 # ajastin, 0.5 s välein tulostetaan näytölle pallon dataa
alku =90 # punainen lähtötasoviivan alku pikseleinä vähän ennen lähtöpaikkaa Pygamessa
loppu = 1000 # punaisen lahtötasoviivan loppupistepikseleinä  Pygamessa
pygame.draw.line(screen,red,(pygamex0-20,pygamey0+r),(pygamex0+1000,pygamey0+r),2) # vaakaviiva
# piirretty niin, että pallo on vaakaviivalla, keskipiste korkeudella r
text = font.render("Lähtökulma = "+str(round(alfadeg))+"  Lähtönopeus = "+str(round(v0)), True, yellow, black)
screen.blit(text,(100,40))
pygame.draw.circle(screen, yellow,(100+x0,500-y0),r) # piirretään pallo
pygame.display.flip() # päivitetään näyttö
pygame.time.delay(2)
pygame.draw.circle(screen,black,(100+x0,500-y0),r) # pallo  pois
pygame.draw.circle(screen,yellow,(100+x1,500-y1),r) # piirretään pallo
pygame.display.flip()
pygame.time.delay(2)
pygame.draw.circle(screen,black,(100+x1,500-y1),r) # pallo pois
# simulointi alkaa aikaan T asti
# eka arvojen (x0,y0) ja (x1,y1) jälkeen tulostuu päivitettyjä (x2,y2)
for j in range (0,T) :
    timer += 1 # puolen sekunnin paikkatietojen tulostusta varten
    x2 = 2*x1 - x0 - (k*vx**2)*Dt**2 # pallon uusi x-koordinaatti
    y2 = 2*y1 - y0 - (g + k*vy**2)*Dt**2 # pallon  uusi y -kordinaatti
    h = y2  # pallon korkeus
    if h > hmax:
        hmax = h
    vx = (x2 - x1)/Dt # päivitetty vaakanopeus
    vy = (y2 - y1)/Dt # päivitetty pystynopeus
    pygame.draw.circle(screen,yellow,(100+x2,500-y2),r) # piirretään pallo
    pygame.display.flip()
    pygame.time.delay(2)
    (xpaikka,ypaikka) = (round(abs(x2)),round(y2)) # paikka muistiin timer-tulostusta varten
    vaakamatka = round(abs(x2)) # muistiin timer-tulostusta varten
    korkeus = round(h) # muistiin timer-tulostusta varten
    x0 = x1 # askeleen jälkeen päivitetään koordinaatit
    x1 = x2
    y0 = y1
    y1 = y2
    roo = 1.225*(288.15/(288.15-0.0065*h))**(-4.25578597) # päivitetty ilman tiheys
    k = k1*roo # pallon päivitetty ilmanvastuskerroin
    if y2 >= hmax : # pallo lakipisteessä
        korkeus = round(h)
        korkeustxt0 = font.render("xxxxxxxxxxxxxxxxxxx",True,black,black)
        korkeustxt = font.render("maxkorkeus = "+str(korkeus)+" m", True, yellow, black)
        screen.blit(korkeustxt0,(100,70)) # puhdistetaan korkeustietoja
        screen.blit(korkeustxt,(100,70))
        hmax = y2
    if   y2  < 0  : # pallo putosi maahan
        break
    pygame.draw.circle(screen, yellow,(100+x2,500-y2),r) # piirretään pallo
    pygame.display.flip()
    pygame.time.delay(2)
    pygame.draw.circle(screen,black,(100+x2,500-y2),r) # edellinen pallon paikka pois
    if timer == 50: #puolisekunti kulunut,näytetään paikkatietoja
        text0 = font.render("xxxxxxxxxxxxxxxx",True,black,black)
        screen.blit(text0,(100,100)) # puhdistetaan vaakamatkatiedot
        text = font.render("Vaakamatka = "+str(vaakamatka),True,yellow,black)
        screen.blit(text,(100,100))
        screen.blit(text0,(100,120)) # puhdistetaan korkeustiedot
        txt = font.render("Lentoaika = "+str(round(100*j*Dt)/100),True,yellow,black)
        screen.blit(txt,(100,120))
        timer = 0