program bounce; uses sysutils; {$R+} {$Assertions on} // ********************************************************** // *** Klicovy datovy typ pro ulozeni informaci o puntiku // ********************************************************** type tPuntik = record x,y : real; // 2D poloha vx,vy : real; // 2D rychlost r : real; // polomer kruznice m : real; // hmotnost (dulezita pri srazce) end; // pomocna funkce oznacujici mista, kde mate doplnit kod function zdeChybi(const s:string):real; begin zdeChybi := 0; assert(false,s); end; // ********************************************************** // *** modifikace rychlosti pri odrazu dotykajicich se puntiku // ********************************************************** procedure odraz(var a,b : tPuntik); {radek c. 27} // https://en.wikipedia.org/wiki/Elastic_collision#Two-dimensional_collision_with_two_moving_objects var skal_soucin,faktor : real; begin Assert(abs( sqr(a.x-b.x)+sqr(a.y-b.y)-sqr(a.r+b.r) ) < 1E-10,'Kruznice se nedotykaji!'); skal_soucin := zdeChybi('skal. soucin vystupujici v (3)'); faktor := zdeChybi('hodnota faktoru Q podle (3)'); a.vx := zdeChybi('x-slozka podle (1)'); a.vy := zdeChybi('y-slozka podle (1)'); b.vx := zdeChybi('x-slozka podle (2)'); b.vy := zdeChybi('y-slozka podle (2)'); end; // ********************************************************** // ** vypocet casu zbyvajiciho do srazky dvou puntiku // ********************************************************** function cas_do_srazky(const a,b : tPuntik):real; var skal_soucin, diskriminant, v2,d2,x2 :real; begin cas_do_srazky := -1; // pokud se v budoucnu nepotkaji vraci funkce zaporne cislo d2 := sqr(a.r+b.r); x2 := zdeChybi('|r1-r2|^2 vystupujici v (5)'); Assert( x2-d2 > -1E-10, 'Disky se prekryvaji!' ); v2 := zdeChybi('|v1-v2|^2 vystupujici v (5)'); if v2=0 then exit; // disky se vuci sobe nepohybuji skal_soucin := zdeChybi('(v1-v2).(r1-r2) vystupujici v (4-5)'); diskriminant := zdeChybi('hodnota D podle (5)'); if diskriminant<=0 then exit; cas_do_srazky := (-skal_soucin-sqrt(diskriminant))/v2; end; // ********************************************************** // *** rychly test spravnosti funkci cas_do_srazky a Odraz // ********************************************************** procedure testSrazeni; const a : tPuntik = ( x:10.0; y: 71.6; vx: 7 ; vy: 1 ; r: 4; m: 1); b : tPuntik = ( x:41.2; y: 20.0; vx: 1 ; vy: 9 ; r: 6; m: 2); tpq = 5.2; eps = 1E-11; var Ek,px,py,Fk,qx,qy,cas:real; begin // 1. kontrola, ze se spravne pocita cas do srazky cas := cas_do_srazky(a,b); a.x:=a.x + a.vx*cas; a.y:=a.y + a.vy*cas; b.x:=b.x + b.vx*cas; b.y:=b.y + b.vy*cas; if abs( tpq - cas ) > eps then begin writeln('Test funkce cas_do_srazky nedopadnul dobre.'); writeln('Misto ',tpq:8:5,' vratila ', cas:8:5 ); halt; end; // 2. kontrola, ze se pri srazce zachovava hybnost a energie px := a.m*a.vx+b.m*b.vx; py := a.m*a.vy+b.m*b.vy; Ek := 0.5*a.m*(sqr(a.vx)+sqr(a.vy))+0.5*b.m*(sqr(b.vx)+sqr(b.vy)); odraz(a,b); qx := a.m*a.vx+b.m*b.vx; qy := a.m*a.vy+b.m*b.vy; Fk := 0.5*a.m*(sqr(a.vx)+sqr(a.vy))+0.5*b.m*(sqr(b.vx)+sqr(b.vy)); if abs( px-qx ) + abs( py-qy ) > eps then begin writeln('Test funkce odraz nedopadnul dobre.'); writeln('Misto [',px:8:5,',',py:8:5,'] vratila po srazce hybnost [',qx:8:5,',',qy:8:5,']' ); halt; end; if abs( Ek-Fk ) > eps then begin writeln('Test funkce odraz nedopadnul dobre.'); writeln('Misto ',Ek:8:5,' vratila po srazce energii ',Fk:8:5 ); halt; end; end; // ********************************************************** // *** Seznam puntiku (glob. promenna) // *** -- blok velkych puntiku // *** -- pocatecni poloha a rychlost maleho puntiku // ********************************************************** const malyP : tPuntik = (x: -119; y: 0; vx: 239; vy: 2; r: 0.5; m: 1); velkyP : tPuntik = (x: 0; y: 0; vx: 0; vy: 0; r: 4; m: 100); const npx = 6; npy = 6; rozestup = 10; type tVsechnyPuntiky = array[0..npx*npy] of tPuntik; var puntiky : tVsechnyPuntiky; procedure provedRozmisteni; var a: tPuntik; i,j,n : integer; begin puntiky[0] := malyP; n := 1; for i := 0 to npx-1 do for j := 0 to npy-1 do begin a := velkyP; a.x := (i+0.5-npx/2)*rozestup; a.y := (j+0.5-npy/2)*rozestup; puntiky[n] := a; n := n+1; {radek c. 133} end; end; // ********************************************************** // *** Vystup do souboru // ********************************************************** var vystup : text; cas : real; cas_vypisu: real; const dt_vypis = 1/16; // Vypis polohy bodu posunute o cas dt za prepokladu pohybu bez srazek procedure vypisPolohyPuntiku(dt:real); var n:integer; begin for n:=0 to high(puntiky) do begin Writeln(vystup, puntiky[n].x+puntiky[n].vx*dt:7:3,' ', puntiky[n].y+puntiky[n].vy*dt:7:3,' ',puntiky[n].r:7:3); end; writeln(vystup); writeln(vystup); end; // ********************************************************** // *** Evoluce // ********************************************************** procedure posunVCase(dt:real); var n:integer; begin for n:=0 to high(puntiky) do begin puntiky[n].x:=puntiky[n].x + puntiky[n].vx*dt; puntiky[n].y:=puntiky[n].y + puntiky[n].vy*dt; end; cas := cas + dt; end; function dalsiSrazka:boolean; var m,n,m1,n1:integer; dt_min, dt: real; const nikdy = 1E20; begin // 1. najdi dvojici disku, ktere se srazi nejdrive dt_min := nikdy; for m:=0 to high(puntiky) do for n:=m+1 to high(puntiky) do begin {radek c. 182} dt := cas_do_srazky(puntiky[m],puntiky[n]); if (dt>0) and (dt cas_vypisu do begin vypisPolohyPuntiku(cas_vypisu-cas); cas_vypisu:= cas_vypisu + dt_vypis; end; posunVCase(dt_min); odraz(puntiky[m1],puntiky[n1]); end; end; // ********************************************************** // *** Hlavni program // ********************************************************** var k:integer; begin testSrazeni; provedRozmisteni; Assign(vystup,'puntiky.txt'); Rewrite(vystup); while dalsiSrazka do; // pokud neverite, ze se jiz nesrazi, lze v animaci pokracovat: //for k:=1 to 200 do vypisPolohyPuntiku(cas_vypisu-cas+k*dt_vypis); close(vystup); end. (* ************************************************************ * Kviz pro studenty * ************************************************************ 1. K cemu je na {radek c. 27} 2. K cemu slouzi n := n+1 na {radek c. 133} 3. Kde se berou meze cyklu na {radek c. 182} 4. V Kterém místě kódu je určena počáteční rychlost malého disku? *) (* ************************************************************ * instrukce pro pouzit GNUPLOT * ************************************************************ # Prikazy pro vykresleni souboru puntiky.txt # nejprve nastavime zpusob zobrazeni a barvy set size ratio -1 set palette model RGB defined ( 0 'red', 4 'blue' ) set cbrange [0:4] unset colorbox # opakovane muzeme pak prikazem (ano dlouhym) spustit animaci stats 'puntiky.txt' nooutput; do for [k=0:STATS_blocks-2] { plot [-80:80][-80:80] 'puntiky.txt' index k using 1:2:3:3 with circles palette fs solid; pause 0.1 } # az bude dobre vypadat, vyrobime animovany gif nasledujicimi prikazy (chvili to trva...) set term gif size 800,800 animate; set output 'puntiky.gif'; do for [k=0:STATS_blocks-2] { plot [-80:80][-80:80] 'puntiky.txt' index k using 1:2:3:3 with circles palette fs solid } ; unset term # komu se nelibi jeden dlouhy radek, lze jej rozdelit pouzitim \ set term gif size 800,800 animate; \ set output 'puntiky.gif'; \ do for [k=0:STATS_blocks-2] { \ plot [-80:80][-80:80] 'puntiky.txt' index k using 1:2:3:3 with circles palette fs solid; \ pause 0.1 \ } ;\ unset term *)