Suomipelit – Fysiikkaa peliohjelmoijille - Osa 1: Sijainti ja liike

1. Johdanto #

Fysiikan mallinnus on tärkeässä osassa kaikissa uusissa peleissä. Realistinen fysiikka lisää peliin jopa enemmän realismia kuin realistinen grafiikka. Se tuo peliin vahvan paikalla olon tunteen. Tämä artikkeli opettaa sinulle fysiikkaa peliohjelmointia varten. Käsittelemme fysiikkaa sekä 3D-avaruudessa että 2D-avaruudessa, jossa asiat muuttuvat huomattavasti helpommiksi. Keskitymme erityisesti kiinteiden kappaleiden mekaniikkaan. Tarkemmin klassiseen eli Newtonin mekaniikkaan. Kiinteä kappale (englanniksi rigid body) on kappale, jonka muoto ei muutu. Se ei siis törmätessään hajoa tai painu kasaan. Mekaniikka taas tarkoittaa kappaleiden välistä vuorovaikusta voimien ja törmäysten välityksellä. Newtonin mekaniikka tarkoittaa mekaniikkaa, joka pätee arkipäisillä mittasuhteilla vaikka luodeille, jalkapalloille, ihmisille ja autoille. Se ei siis toimi atomien tai galaksien mittakaavassa.

Tämä artikkelisarja tulee alustavien suunnitelmien mukaan sisältämään seuraavat osat.

1. Sijainti ja liike (tämä artikkeli)
2. Kiihtyvyys ja voimat
3. Törmäykset
4. Lepokontakti ja kitka

Asioiden yksinkertaistamiseksi oletamme, että kaikki kappaleet ovat joko palloja (ympyröitä 2D-avaruudessa) tai suorakulmaisia särmiöitä eli laatikoita (suorakulmioita 2D-avaruudessa). Vaikka esittelemme vain nämä kaksi muotoa, voidaan tekniikat yleistää minkä muotoisiin kappaleisiin tahansa. Tämä artikkeli ei ota kantaa siihen kuinka palloja tai laatikoita piirretään.

Luomme tässä artikkelissa kaksi luokkaa tai pidemminkin rakennetta: RIGIDBODY2D ja RIGIDBODY3D. Nämä rakenteet pitävät sisällään kaikki tarvittavat tiedot kiinteästä kappaleesta sen simuloimiseen. Alussa, kun rakenteet sisältävät vain kappaleen muodon ja koon, ne näyttävät tältä:

typedef struct
{
  int shape;   // Muoto. 0=ympyrä ja 1=suorakulmio 
  float radius[2];  // ”Säde” metreinä.
} RIGIDBODY2D;

typedef struct
{
  int shape;  // Muoto. 0=pallo ja 1=laatikko
  float radius[3];  // ”Säde” metreinä.
} RIGIDBODY3D;

2. Sijainti ja asento #

Ensimmäinen käsiteltävä asia on, kuinka esitämme kappaleen sijaintia ja asentoa 2D- tai 3D-avaruudessa. Tämä saattaa tuntua itsestään selvältä, mutta jo kappaleen asennon esittäminen 3D-avaruudessa on ongelmallinen.

Kiinteän kappaleen sijaintia (englanniksi position) esitetään sen (massa)keskipisteen koordinaattien avulla. Ei siis voisi olla paljon helpompaa. 2D-avaruudessa tarvitaan kaksi lukua: x ja y-koordinaatti ja 3D-avaruudessa vielä kolmas z-koordinaatti. Sijaintikoordinaattien yksikkönä käytämme metriä (lyhenne m). Kun lisäämme nämä tiedot rakenteisiimme niistä tulee seuraavanlaiset:

typedef struct
{
  int shape;   // Muoto. 0=ympyrä ja 1=suorakulmio 
  float radius[2];  // "Säde" metreinä.
  float position[2];  // Sijainti metreinä
} RIGIDBODY2D;

typedef struct
{
  int shape;  // Muoto. 0=pallo ja 1=laatikko
  float radius[3];  // "Säde" metreinä.
  float position[3];  // Sijainti metreinä
} RIGIDBODY3D;

typedef struct
{
  int shape;  // Muoto. 0=ympyrä ja 1=suorakulmio 
  float radius[2];  // Säde" metreinä.
  float position[2];  // Sijainti metreinä
  float orientation;  // Asento radiaaneina
} RIGIDBODY2D;

typedef struct
{
  int shape;  // Muoto. 0=pallo ja 1=laatikko
  float radius[3];  // "Säde" metreinä.
  float position[3];  // Sijainti metreinä
  float orientation[4];  // Asento kvaterniona
} RIGIDBODY3D;

void multQuaternion(float result[4],
                    const float q1[4], const float q2[4])
{
  result[0] = q1[0]*q2[0] - q1[1]*q2[1] – 
              q1[2]*q2[2] - q1[3]*q2[3];
  result[1] = q1[0]*q2[1] + q1[1]*q2[0] + 
              q1[2]*q2[3] - q1[3]*q2[2];
  result[2] = q1[0]*q2[2] + q1[2]*q2[0] + 
              q1[3]*q2[1] - q1[1]*q2[3];
  result[3] = q1[0]*q2[3] + q1[3]*q2[0] + 
              q1[1]*q2[2] - q1[2]*q2[1];
}

float normalizeQuaternion(float q[4])
{
  float m=sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]);
  if (m!=0)
  {
    q[0]/=m;
    q[1]/=m;
    q[2]/=m;
    q[3]/=m;
  }
  return m;
}

3. Johdatus differentiaaliyhtälöihin #

Differentiaaliyhtälöt ovat yhtä vaikeita kuin miltä kuulostavatkin. Onneksi meidän ei tarvitse yrittää ratkaista sellaista, mutta meidän täytyy tietää mitä ne ovat. Kun olen esitellyt differentiaaliyhtälöt, saatat miettiä mihin niitä oikein tarvitaan. Palaset loksahtavat kuitenkin lopulta paikalleen.

Ensimmäinen ja tärkein asia on käsite nimeltä derivaatta (englanniksi derivative). Derivaatta tarkoittaa muutosta. Se ilmaiseen miten jokin muuttuu suhteessa johonkin toiseen. Otetaan esimerkiksi suora y(x)=2x. Ks. kuva.

Funktio		Derivaatta
kx		k
x^k		kx^(k-1)
Sin x		Cos x

4. Numeerinen integrointi #

Alkuperäisen funktion päätteleminen sen derivaatasta, eli integrointi, on siis varsin hankala prosessi. Tarkoitusperiämme varten emme ole kuitenkaan kiinnostuneita koko alkuperäisestä funktiosta vaan ainoastaan sen arvosta tietyssä pisteessä.

Integrointi ja derivointi eivät ole aivan täydellisesti vastakkaisia toimintoja. Derivoinnissa nimittäin menetetään vakiotermi, jota integroinnissa ei voi palauttaa. Esim. funktioiden x^2+1, x^2+2, x^2+3 jne. kaikkien derivaatta on 2x. Emme täten voi derivaatasta 2x mitenkään päätellä täydellisesti mikä alkuperäinen funktio oli. Se on x^2+jotain. Tämän takia differentiaaliyhtälön ratkaisussa tarvitaan derivaatan lisäksi tietää yksi piste funktiolta f(x) ns. alkuarvo (englanniksi initial value). Täytyy siis tietää mitä on f(0). Jos tiedämme, että funktion f(x) derivaatta on 2x ja f(0) on 1, voimme yksikäsitteisesti päätellä, että f(x) on x^2+1, koska nyt f(0) on 0^2+1=1, eli alkuarvo toteutuu.

Numeerinen integraattori on tietokoneohjelma, joka saa syötteenään funktion f(x) derivaatan f’(x), alkuarvon f(0) ja pyrkii näiden perusteella arvioimaan paljonko f(x) on jollain tietyllä x:n arvolla. Esittelen kolme erilaista numeerista integraattoria helpoimmasta vaikeampaan. Ne ovat Eulerin menetelmä, keskipistemenetelmä ja Runge-Kutta-menetelmä. Katsotaanpa ensin miltä numeerisen integraattorin prototyyppi näyttäisi.

float integrate( float initialValue, float x,
                 float (*derive)(float value, float x) );

void integrate( float *result, const float *initialValue,                         
                DWORD size, float x, 
                void (*derive)(float *result, 
                               const float *value,
                               float x, void *userData), 
                float *temp, void *userData );

void integrate( float *result, const float *initialValue,                         
                DWORD size, float x, 
                void (*derive)(float *result, 
                               const float *value,
                               float x, void *userData), 
                float *temp, void *userData )
{
  int i;
  
  // Kysy paljonko f'(0) on
  derive( temp, initialValue, 0, userData );
  
  // Laske f(0)+x*f'(x)
  for ( i=0; i<size; i++ )    
    result[i]=initialValue[i]+x*temp[i];
}

void integrate( float *result, const float *initialValue,                         
                DWORD size, float x, 
                void (*derive)(float *result, 
                               const float *value,
                               float x, void *userData), 
                float *temp, void *userData )
{
  int i;
  float *temp2;
  temp2=temp+size;
 
  // Tutki Eulerin menetelmän avulla paljonko on f(x/2)
  // Kysy paljonko on f'(0)
  derive( temp, initialValue, 0, userData );
  // laske f(0)+(x/2)*f'(0)
  for ( i=0; i<size; i++ )   
    temp2[i]=initialValue[i]+(x/2.0)*temp[i];
 
  // Nyt voimme soveltaa keskipistemenetelmää
  // Kysy paljonko on f'(x/2), kun f(x/2) tiedetään
  derive( temp, temp2, x/2.0, userData );
  // Laske paljonko on f(0)+x*f'(x/2)
  for ( i=0; i<size; i++ ) 
    result[i]=initialValue[i]+x*temp[i];
}

void integrate( float *result, const float *initialValue,                         
                DWORD size, float x, 
                void (*derive)(float *result, 
                               const float *value,
                               float x, void *userData), 
                float *temp, void *userData )
{
  int i;
  float *temp2, *k1, *k2, *k3, *k4;

  temp2=temp+size;
  k1=temp+2*size;
  k2=temp+3*size;
  k3=temp+4*size;
  k4=temp+5*size;

  derive( temp, initialValue, 0, userData );
  for ( i=0; i<size; i++ )
    k1[i]=x*temp[i];

  for ( i=0; i<size; i++ ) 
    temp2[i]=initialValue[i]+0.5*k1[i];
  derive( temp, temp2, 0.5*x, userData );
  for ( i=0; i<size; i++ )
    k2[i]=x*temp[i];

  for ( i=0; i<size; i++ )   
    temp2[i]=initialValue[i]+0.5*k2[i];
  derive( temp, temp2, 0.5*x, userData );
  for ( i=0; i<size; i++ )
    k3[i]=x*temp[i];

  for ( i=0; i<size; i++ )
    temp2[i]=initialValue[i]+k3[i];
  derive( temp, temp2, x, userData );
  for ( i=0; i<size; i++ )
    k4[i]=x*temp[i];

  for (i=0; i<size; i++)
    result[i]=initialValue[i]+
              (1.0/6.0)*k1[i]+
              (1.0/3.0)*k2[i]+
              (1.0/3.0)*k3[i]+
              (1.0/6.0)*k4[i];
}

void derive( float *result,
             const float value, 
             float x, void *userData )
{
  result[0]=2*x;  // Funktion x^2+3 derivaatta on 2x
}

float initialValue=3;
float output;
float temp[6];
integrate(&output, &initialValue, 1, 0.5, derive, temp,       
          NULL);

Funktio			Derivaatta	Alkuarvo f(0)
2x+3			2		3
x^2+5x+2		2x+5		2
3*x^3+4*x^2+8*x+7	9x^2+8x+8	7

x:			0	0.25	0.5	0.75	1
Oikea tulos		3	3.5	4	4.5	5
Eulerin menetelmä	3	3.5	4	4.5	5
Keskipistemenetelmä	3	3.5	4	4.5	5
Runge-Kutta		3	3.5	4	4.5	5

x:			0	0.25	0.5	0.75	1
Oikea tulos		2	3.31	4.75	6.31	8
Eulerin menetelmä	2	3.25	4.5	5.75	7
Keskipistemenetelmä	2	3.31	4.75	6.31	8
Runge-Kutta		2	3.31	4.75	6.31	8

x:			0	0.25	0.5	0.75	1
Oikea tulos		7	9.3	12.38	16.52	22
Eulerin menetelmä	7	9	11	13	15
Keskipistemenetelmä	7	9.29	12.28	16.20	21.25
Runge-Kutta		7	9.3	12.38	16.52	22

5. Tasainen liike #

Kuten olet ehkä huomannutkin, tämän artikkelin ja tämän luvun otsikossa esiintyy sana liike. Liike tarkoittaa sitä, että kappaleet eivät välttämättä pysy paikallaan vaan niiden sijainti muuttuu ajan suhteen. Siis niiden sijainti MUUTTUU ajan SUHTEEN. Tässä vaiheessa minä vähän elättelin varovaista toivetta, että lukija olisi kokenut seuraavanlaisen ahaa-elämyksen: "Hetkinen! Tässähän jokin (sijainti) muuttuu suhteessa johonkin (aika). Sehän on sama asia kuin derivaatta!" Sijainnilla on siis derivaatta suhteessa aikaan. Sijainnin aikaderivaatta on sen verran yleinen ilmiö, että sillä on arkikielessä ihan oma sana: nopeus (englanniksi velocity). Autossa ei siis sanota olevan sijainnin aikaderivaattamittaria vaan nopeusmittari. Nopeutta merkitään yleensä symbolilla v ja sen yksikkö on metriä sekunnissa (lyhenne m/s). Jos esitämme sijaintia symbolilla s ja aikaa symbolilla t voimme esittää lausetta "sijainnin aikaderivaatta on nopeus" kaavalla s'(t)=v. Artikkelisarjan tässä osassa oletamme liikkeen olevan tasaista eli nopeus pysyy kokoajan vakiona. 2D-avaruudssa nopeus on kaksi komponenttinen vektori, jonka suunta esittää liikkeen suunnan ja pituus vauhtia. Vastaavasti 3D-avaruudessa käytämme kolme komponenttista vektoria. Seuraavassa nopeudella täydennetyt rakenteet.

typedef struct
{
  int shape;  // Muoto. 0=ympyrä ja 1=suorakulmio 
  float radius[2];  // ”Säde” metreinä.
  float position[2];  // Sijainti metreinä
  float orientation;  // Asento radiaaneina
  float velocity[2];  // Nopeus metreinä sekunnissa
} RIGIDBODY2D;

typedef struct
{
  int shape;  // Muoto. 0=pallo ja 1=laatikko
  float radius[3];  // ”Säde” metreinä.
  float position[3];  // Sijainti metreinä
  float orientation[4];  // Asento kvaterniona
  float velocity[3];  // Nopeus metreinä sekunnissa
} RIGIDBODY3D;

void simulate2D(RIGIDBODY2D *b, float t)
{
  static float initialValue[2];
  static float output[2];
  static float temp[6*2];
  // Kopioi nykyinen sijainti alkuarvoksi.
  initialValue[0]=b->position[0];
  initialValue[1]=b->position[1];
  // Interoi aika-askeleen verran. Koska ”derive” tarvitsee
  // Kappaleen nopeuden, Annamme osoittimen kappaleeseen
  // ”userData”-parametrina.
  integrate(output, initialValue, 2, t, derive, 
             temp, (void *)b);
  // Kopioi tulos uudeksi sijainniksi
  b->position[0]=output[0];
  b->position[1]=output[1];
}

void derive( float *result,
             const float *value, 
             float x, void *userData )
{
  // Osoitin kappaleeseen saatiin “userData”:n kautta
  RIGIDBODY2D *b=(RIGIDBODY2D *)userData;
  // Sijainnin derivaatta on nopeus.
  result[0]=b->velocity[0];
  result[1]=b->velocity[1];
}

typedef struct
{
  int shape;  // Muoto. 0=ympyrä ja 1=suorakulmio 
  float radius[2];  // ”Säde” metreinä.
  float position[2];  // Sijainti metreinä
  float orientation;  // Asento radiaaneina
  float velocity[2];  // Nopeus metreinä sekunnissa
  float angularVelocity;  // Kulmanopeus radiaaneina sekunnissa
} RIGIDBODY2D;

void simulate2D(RIGIDBODY2D *b, float t)
{
  static float initialValue[3];
  static float output[3];
  static float temp[6*3];
  // Kopioi nykyinen sijainti ja asento alkuarvoksi.
  initialValue[0]=b->position[0];
  initialValue[1]=b->position[1];
  initialValue[2]=b->orientation;
  // Interoi aika-askeleen verran. Koska ”derive” tarvitsee
  // kappaleen nopeuden, Annamme osoittimen kappaleeseen
  // ”userData”-parametrina.
  integrate(output, initialValue, 3, t, derive, 
             temp, (void *)b);
  // Kopioi tulos uudeksi sijainniksi ja asennoksi
  b->position[0]=output[0];
  b->position[1]=output[1];
  b->orientation=output[2];
}

void derive( float *result,
             const float *value, 
             float x, void *userData )
{
  // Osoitin kappaleeseen saatiin “userData”:n kautta
  RIGIDBODY2D *b=(RIGIDBODY2D *)userData;
  // Sijainnin derivaatta on nopeus.
  result[0]=b->velocity[0];
  result[1]=b->velocity[1];
  // Asennon derivaatta on kulmanopeus
  result[2]=b->angularVelocity;
}

typedef struct
{
  int shape;  // Muoto. 0=pallo ja 1=laatikko
  float radius[3];  // ”Säde” metreinä.
  float position[3];  // Sijainti metreinä
  float orientation[4];  // Asento kvaterniona
  float velocity[3];  // Nopeus metreinä sekunnissa
  float angularVelocity[3]; // Kulmanopeus radiaaneina sekunnissa
} RIGIDBODY3D;

void simulate3D(RIGIDBODY3D *b, float t)
{
  static float initialValue[7];
  static float output[7];
  static float temp[6*7];
  // Kopioi nykyinen sijainti ja asento alkuarvoksi.
  initialValue[0]=b->position[0];
  initialValue[1]=b->position[1];
  initialValue[2]=b->position[2];
  initialValue[3]=b->orientation[0];
  initialValue[4]=b->orientation[1];
  initialValue[5]=b->orientation[2];
  initialValue[6]=b->orientation[3];
  // Integroi aika-askeleen verran. Koska”derive” tarvitsee
  // Kappaleen nopeuden, Annamme osoittimen kappaleeseen
  // ”userData”-parametrina.
  integrate(output, initialValue, 7, t, derive, 
             temp, (void *)b);
  // Kopioi tulos uudeksi sijainniksi ja asennoksi
  b->position[0]=output[0];
  b->position[1]=output[1];
  b->position[2]=output[2];
  b->orientation[0]=output[3];
  b->orientation[1]=output[4];
  b->orientation[2]=output[5];
  b->orientation[3]=output[6];
  normalizeQuaternion(b->orientation);
}

void derive( float *result,
             const float *value, 
             float x, void *userData )
{
  static float tempQ[4];
  // Osoitin kappaleeseen saatiin “userData”:n kautta
  RIGIDBODY3D *b=(RIGIDBODY3D *)userData;
  // Sijainnin derivaatta on nopeus.
  result[0]=b->velocity[0];
  result[1]=b->velocity[1];
  result[2]=b->velocity[2];
  // Asennon derivaatta on 0.5*kulmanopeus*asento
  tempQ[0]=0;
  tempQ[1]=0.5*b->angularVelocity[0];
  tempQ[2]=0.5*b->angularVelocity[1];
  tempQ[3]=0.5*b->angularVelocity[2];
  multQuaternion(&result[3], tempQ, &value[3]);
}

6. Miksi se piti tehdä niin hankalasti? #

Tässä vaiheessa mieleen saattaa tulla, että miksi hitossa se piti tehdä noin hankalasti? Jos olet koskaan ohjelmoinut peliin jotain liikkuvaa, olet todennäköisesti käyttänyt kaavaa uusi_sijainti=vanha_sijainti+nopeus*aika. Olet ehkä jopa keksinyt tämän kaavan itse. Jos nyt vertaat tätä kaavaa Eulerin menetelmään f(x)=f(0)+f’(0)*x ja toteat sen tosiasian, että tässä tapauksessa f(x) on uusi_sijainti, f(0) on vanha_sijainti, f’(0) on nopeus ja x on aika, niin huomaat, että kyseessä ovat täsmälleen samat kaavat! Olet siis itse asiassa integroinut käyttäen Eulerin menetelmää, ilmeisesti jopa tietämättäsi. Tasaisen liikkeen tapaus on itse asiassa niin yksinkertainen, että Eulerin menetelmä paitsi antaa aivan oikean lopputuloksen, myös voidaan kirjoittaa helposti kaavaksi ilman varsinaista integraattoria. Artikkelisarjan seuraavissa osissa siirrymme monimutkaisempiin tilanteisiin, joissa Eulerin menetelmä alkaa antaa virheellisiä arvoja ja asiat on helpompi laskea käyttäen yleistä integraattoria.

7. Esimerkkiohjelma #

Esimerkkiohjelmamme ei ole päällepäin sangen kiintoisa: vain muutamia liikkuvia kappaleita. Ohjelman kiintoisin asia on pellin alla hyrräävä integraattorisysteemi.

Kuten jo johdannossa ilmoitin, en ota kantaa siihen kuinka kappaleita piirretään. Fysiikallinen simulaatio ilman visualisaatiota ei kuitenkaan olisi kovin kummoinen, joten tein avuksi tiedostot 2dpiirto.h ja 3dpiirto.h. Ne sisältävät muutamia OpenGL-rajapintaa käyttäviä funktioita grafiikan piirtämiseen, joten kääntääksesi ohjelmat sinun pitää linkittää myös kirjastojen opengl32.lib ja glu32.lib kanssa. Lisää OpenGL:stä löydät halutessasi artikkelisarjasta "OpenGL:n perusteet". Yhtä hyvin grafiikan piirtämiseen voitaisiin käyttää mitä tahansa muuta rajapintaa/grafiikkakirjastoa kuten SDL tai Direct3D.

2D-esimerkki screenshot

3D-esimerkki screenshot

8. Loppusanat #

Tässä osassa opit kuinka kappaleen asentoa voidaan esittää 2D- ja 3D-avaruudessa. Opit lisäksi numeerisesta integroinnista ja kuinka sitä sovelletaan tasaisen liikkeen laskemiseen. Seuraavassa osassa puhumme tilanteesta, jossa nopeus ei enää ole vakio eli puhumme kiihtyvyydestä ja voimista.

Raportoithan kaikki tästä artikkelista löytämäsi virheet (niin kirjoitus-, kuin asiavirheetkin) osoitteeseen markus.ilmola@pp.inet.fi , niin korjaan ne mahdollisimman nopeasti. Myös kaikki kommentit ja kysymykset ovat tervetulleita.

Lopuksi vielä linkkejä, joista löydät lisälukemista tässä artikkelissa opituista asioista. Varoitettakoon kuitenkin, että ne saattavat sisältää vaikeasti luettavia matemaattisia merkintöjä.

Katkelma kirjasta "Numerical Recipes in C: The Art of Scientific Computing":
http://www.library.cornell.edu/nr/bookcp...

David Barffin fysiikka artikkelisarjan liite:
http://www-2.cs.cmu.edu/~baraff/sigcours...