Simulace oběhu planety v Pascalu – Pascal – Fórum – Programujte.com

Zdravím!

Peru se se zápočtovým programem, který má zadání: "Simulace oběhu 1 planety a 1 jejího měsíce kolem Slunce." (začala jsem zatím pouze planetou). Samotné použití Pascalu v tomto případě mi přijde nanejvýš nešťastné, pravděpodobně už ho také nikdy jindy nepoužiji, ale tohle je podmínka...

Podařilo se mi snad přijít na to, jak by měl program fungovat (krokovat), nechala jsem si vypisovat polohy a měly by zhruba odpovídat. Ale jakmile použiju grafiku, tak je problém s velkými čísly - jistě při vykreslení jsem používala měřítko, ale narážím stále na chybu 217: "Unhandled Exception Occurred" během.

Tedy alespoň myslím, že je chyba ve velkých (a malých) číslech.

Nemá někdo zkušenosti s takovouto grafikou?

 
program simulace;

uses Graph, Crt, Math;

var Gd, Gm, m, n, Radius: Integer;
    kappa, m_Sl, vx0, vy0, ax, ay, perihel, x0, y0, x, y, vx, vy, F: Extended;
    dt, k, l, kk, ll: LongInt;

begin
        Gd:=Detect;

        m:=0;                  {posunutí středu soustavy souřadné "doprava"}
        n:=0;                  {posunutí středu soustavy souřadné "dolů"}

        perihel:=0.9141*(1.49598E+11);
        dt:= 60;
        m_Sl:= 2E+19;          {hmotnost Slunce}
        kappa:= 6.67;          {gravitacni konstanta}
        vx0:= 0;
        vy0:= 44000;
        x0:= -perihel;
        y0:= 0;

        InitGraph(Gd, Gm, '');

        if GraphResult <> grOK then
            Halt(1);

        Setcolor (14);                           {vykreslení SLUNCE}
              For Radius:=0 to 20 do
                 begin
                 Circle (0,0,Radius);
                 end;

        k:=round(x0);           {zaokrouhlení na celé pixely}
        l:=round(y0);
        Circle (k, l, 3);

        repeat
            begin
                F:=(-kappa*m_Sl)/power(x0*x0+y0*y0, 1.5);

                ax:=F*x0;
                ay:=F*y0;

                vx:=vx0+ax*dt;
                vy:=vy0+ay*dt;

                x:=x0+vx*dt;
                y:=y0+vy*dt;

                kk:=round(x);
                ll:=round(y);
                Circle (kk, ll, 3);

                x0:=x;
                y0:=y;
                vx0:=vx;
                vy0:=vy;

                Delay (1000);
            end;
         Until Keypressed;
         CloseGraph;
end.

Nejsem programátor, tohle je pro mě jen nutné zlo. Ráda si sestavím nějaký program, když se mi chce a mám nápad. Ale u tohohle mě tlačí čas a vyučující :/

Předem díky za rady!

#1 Anetace
Teď nemám čas se na to kouknout, ale jen tak zběžně bych si tipnul, že nepřepočítáváš reálné hodnoty pozice planety do rozsahu obrazovky...

#1 Anetace
Koukám, že tam máš úplně nesmyslná čísla. Slunce je mnohem těžší, gravitační konstanta menší. Zemi máš sice ve správné vzdálenosti, ale je tak rychlá, že by uletěla a stal by se z ní nomád. Zkus při ladění používat hodnoty z reálného světa a teprve potom laborovat.

Význam výrazu power(x0*x0+y0*y0, 1.5) mi zcela uniká. Zkus používat standardní fyzikální vzorečky.

#3 Kit
tipuju, že chtěla spočítat přeponu pravoúhlého trojúhelníka a nikdo jí neřekl, že pascal umí druhou mocninu i odmocninu, ale vzpomněla si na tu fintu, že odmocnina se dá spočítat jako mocnina na 1/2 akorát nějak popletla tu polovinu?

#2 JoDiK
Takže opavdu je tomu tak, pokoušíš se zaokrouhlovat čísla někde x10^-15, což se funkci Round už nelíbí, navíc je to nesmysl, protože by všechno vyšlo 0.

Doporučuju si nejdřív to vykreslování nazkoušet na něčem jednoduchém, třeba po kružnici a postupně zapojovat  složitější výpočty...

#2 JoDiK
Ano, to jsem psala. Zkoušela jsem měnit rozměry při vykreslení, aby se vešly na obrazovku, ale hlásí to pořád stejnou chybu. 

#6 Anetace
Ve kterém řádku ti to hlásí tu chybu?

#3 Kit
Ale to jsou standardní fyzikální vzorečky :) Je to Newtonův gravitační zákon...
A co se týče hmotnosti Slunce a gravitační konstanty, vzhledem k tomu, že se násobí, Slunce má hmotnost 2E+30 a gravitační konstanta je 6,67E-11, tak jsem to zkrátila, aby hmotnost Slunce nebyla tak velká a abych nemusela pracovat s tak malým číslem. 

V perihelu má Země maximální rychlost a ta je 31 km/s, tzn. 31000 m/s. Pravda, v této verzi mám rychlost moc velkou, protože jsem v průběhu zkoušela i jiné planety, tak mi tam zůstala špatná. 

#4 JoDiK
Inspirovala jsem se tímto textem na straně 12. Dává to smysl...

#7 JoDiK
Program se spustí, ale ihned se ukončí s chybou 217: Unhandled exception occurred

#5 JoDiK
Jo jo, to jsem si myslela. Já jsem si s grafikou hrála už předtím, jednodušeji, ale vyučující mi to hodil na hlavu, že to se simulací nemá nic společného. Můžu ukázat... dokázala jsem, aby pomocí Keplerovy rovnice planety obíhala kolem Slunce, stejně tak zároveň Měsíc kolem planety, ale on vyloženě vyžaduje, aby byl použít Newtonův gravitační zákon :( Bohužel, tam takto vysoká a malá čísla.

program simulace;
uses Graph, Crt;

var Gd, Gm, m, n, a, b, c, d, f, g, i, j, Radius, x, y: Integer;
    E, Z, pomX, pomY, pomI, pomJ: real;

    
begin

       Gd:=Detect;
       a:=400;
       b:=200;
       c:=40;
       d:=50;
       m:=900;
       n:=350;
       f:=200;
       g:=0;
       InitGraph(Gd, Gm, '');

       if GraphResult <> grOK then
            Halt(1);

        Setcolor (14);                           {SLUNCE}
        For Radius:=0 to 20 do
            begin
            Circle (m-a,n,Radius);
            end;



        E:=0;                                   {počáteční "rychlost" planety}
        Z:=0;                                   {počáteční "rychlost" měsíce}

        Repeat
            pomX:= (a*-cos(E))-f;
            x:=Trunc(pomX)+m;             {x-ová souřadnice polohy planety}
            pomY:= b*-sin(E);
            y:=Trunc(pomY)+n;             {y-ová souřadnice polohy planety}
            
                PutPixel (x,y,Green);     {vykreslení dráhy planety}
                Setcolor(Green);
                Circle (x,y,10);          {vykreslení planety}

            Delay(20);

            E:=E+0.01 ;                   {"rychlost" planety}

            pomI:= (c*-cos(Z))-g;
            i:=Trunc(pomI);             {x-ová souřadnice polohy měsíce}
            pomJ:= d*-sin(Z);
            j:=Trunc(pomJ);             {y-ová souřadnice polohy měsíce}

                PutPixel (x+i,y+j,Lightblue);    {vykreslení dráhy měsíce}
                Setcolor(Lightblue);
                {Ellipse (x, y, 0, 360, c, d);}
                Circle (x+i,y+j,3);               {vykreslení měsíce}

            Delay(20);

            Setcolor(Black);
            Circle (x,y,10);             {přemaže stávající polohu planety}
            Circle (x+i,y+j,3);          {přemaže stávající polohu měsíce}
            {Ellipse (x, y, 0, 360, c, d);}

            Z:=Z+0.1336;                  {"rychlost" měsíce}
            
            
        until Keypressed;
            




      CloseGraph;

end.

#9 Anetace
Tak už je to jasné r/r**3 totiž není r**1.5, ale 1/r**2

#12 Kit
Ale to tučné r je vektor polohy planety, který je potom nahrazen souřadnicemi x a y, tedy se nemůže zkrátit.

#9 Anetace
Asi si nerozumíme, ty výsledky asi mají smysl, ale je nesmyslné je zaokrouhlovat a pokoušet se v těchto hodnotách něco vykreslovat.

Většinou se to dělá tak, že si zjistíš maximální rozsah těch výsledků a podle nich si vypočítáš měřítko. Tímto měřítkem budeš násobit všechny vypočítané výsledky a pak teprve budeš zaokrouhlovat. Tzn. že například vzdálenost planety od slunce ti vyjde mezi 1,2x10^8 a 1,6x10^8 a grafika má rozsah 800x600. Tak si určíš měřítko třeba 1x10^-6 a tímto vynásobíš každou vypočítanou vzdálenost a pak teprve zaokrouhlíš. Tim se dostaneš do rozsahu obrazovky.

#11 Anetace
Doporučuji nepoužívat { komentáře }, ale (* komentáře *), protože složené závorky dost pletou programátory, kteří píší v jiných jazycích.

Kromě toho se komentáře spíš používají, když programátor chce označit chybu ve svém programu.

#15 Kit
Dobře, díky. Já nikdy v ničem jiném neprogramovala a ty komentáře jsou víceméně jen proto, abych se v tom vyznala já. Příště  se poučím :)

#14 JoDiK
Dobře, ještě se o to pokusím. Ale bohužel až večer, teď musím pryč. Zatím děkuji za rady!

#16 Anetace
Učitelé většinou nutí psát komentáře, ale trendem v programování je psát tak čitelně, aby komentáře vůbec nebyly potřebné. Názvy proměnných x, x0, vx, ax jsou výstižné (ze vzorečků), ale například s názvy a, b, c, d je to podstatně horší.

Neboj se nazvat proměnnou (spíš konstantu) názvem GRAVITACNI_KONSTANTA a dát ji do bloku const. Stejně tak i HMOTNOST_SLUNCE. Tím úplně odpadne potřeba komentářů tohoto typu.

Tak jo, díky vaší pomoci (hlavně JoDik) se mi podařilo rozběhat planetu kolem Slunce :) Díky!
Teď ještě kolem té planety měsíc... s tím ti teď hraju a hrát ještě budu. Říkám si, že to už nemůže být složité, když už mám ten mechanismus oběhu, ale není to tak jednoduché, jak se zdá. Pohyb okolo pohybujícího se tělesa... nic moc. Pokud se mi nepodaří, nejspíš vás ještě poprosím, jestli byste se na to nekoukli. Nejdřív se ale budu snažit sama.

Tak mi funguje oběh planety kolem Slunce a ač se zdá oběh měsíce analogický, nedaří se mi :( Možná bude chyba v posunu souřadnic, možná v tom měřítku, možná v rychlosti vykreslení... Určitě vím, že kdybych zachovala správné měřítko u obou těles, tak by v poměru k oběžné dráze planety měla dráha měsíce průměr asi 2 pixely, takže tam to bude chtít měřítko trochu změnit. Je tam něco, co vás praští do očí? Něco co je hned špatně? Nepsala bych, ale už jsem celkem zoufalá a fakt mě to nebaví :/

program simulace;

uses Graph, Crt, Math;

var Gd, Gm, m, n, Radius, i: Integer;
    dt, dt_mesic, AU, kappa, HmotnostSlunce, HmotnostPLanety, vx0, vx0_mesic, vy0, vy0_mesic, ax, ax_mesic, ay, ay_mesic, perihel, perigeum, x0, x0_mesic, y0, y0_mesic, x, x_mesic, y, y_mesic, vx, vx_mesic, vy, vy_mesic, F, F_mesic: Extended;

begin
        Gd:=Detect;

        m:=650;                  (*posunutí středu soustavy souřadné "doprava"*)
        n:=350;                  (*posunutí středu soustavy souřadné "dolů"*)

        AU:=1.49598E+11;
        perihel:=0.9141*AU;
        perigeum:=0.0024*AU;
        dt:=40000;
        dt_mesic:=40000;
        HmotnostSlunce:= 2E+30;
        HmotnostPlanety:= 6E+24;
        kappa:= 6.67E-11;
        vx0:= 0;
        vy0:= -31000;
        x0:= perihel;
        y0:=0;
        vx0_mesic:=0;
        vy0_mesic:=-1082;
        x0_mesic:=x0+perigeum;
        y0_mesic:=y0;
        
        i:=1;

        InitGraph(Gd, Gm, '');

        if GraphResult <> grOK then
            Halt(1);

        Setcolor (14);                           (*vykreslení SLUNCE*)
              For Radius:=0 to 10 do
                 begin
                 Circle (m,n,Radius);
                 end;

        PutPixel (m-round(x0*2E-09), n-round(y0*2E-09), Green);
        (*PutPixel (m-round(x0_mesic*2.2E-09), n-round(y0_mesic*2.2E-09), Red);*)

        WriteLn('x_0 :',round(x0*3E-09));
        WriteLn('y_0 :',round(y0*3E-09));
        
        (*WriteLn('x_0_mesic: ',m-round(x0_mesic*2.2E-09));
        WriteLn('y_0_mesic: ',n-round(y0_mesic*2.2E-09));*)
        
        repeat
            begin
                F:=-kappa*HmotnostSlunce/power((m-x0)*(m-x0)+(n-y0)*(n-y0), 1.5);
                ax:=F*x0;
                ay:=F*y0;
                vx:=vx0+ax*dt;
                vy:=vy0+ay*dt;
                x:=x0+vx*dt;
                y:=y0+vy*dt;

                (*F_mesic:=-kappa*HmotnostPlanety/power((m-x-x0_mesic)*(m-x-x0_mesic)+(n-y-y0_mesic)*(n-y-y0_mesic), 1.5);
                ax_mesic:=F_mesic*x0_mesic;
                ay_mesic:=F_mesic*y0_mesic;
                vx_mesic:=vx0_mesic+ax_mesic*dt_mesic;
                vy_mesic:=vy0_mesic+ay_mesic*dt_mesic;
                x_mesic:=x0_mesic+vx_mesic*dt_mesic;
                y_mesic:=y0_mesic+vy_mesic*dt_mesic;
                
                PutPixel (round((x-x_mesic)*2E-09), round((y-y_mesic)*2E-09), Red);*)

                SetColor(Green);
                Circle (m-round(x*2E-09), n-round(y*2E-09), 3);
                SetColor(Black);
                Circle(m-round(x*2E-09), n-round(y*2E-09), 3);
                Putpixel(m-round(x*2E-09), n-round(y*2E-09), Green);

                WriteLn('x_',i,': ',round(x*2E-09));
                WriteLn('y_',i,': ',round(y*2E-09));

                (*WriteLn('x_mesic_',i,': ',round((x-x_mesic)*2E-09));
                WriteLn('y_mesic_',i,': ',round((y-y_mesic)*2E-09));*)

                x0:=x;
                y0:=y;
                vx0:=vx;
                vy0:=vy;
                (*x0_mesic:=x_mesic;
                y0_mesic:=y_mesic;
                vx0_mesic:=vx_mesic;
                vy0_mesic:=vy_mesic; *)
                i:=i+1;

                Delay (10);
            end;

         Until Keypressed;
         CloseGraph;
end.

Zakomentované jsou moje pokusy o měsíc...

#20 Anetace
Nezapomeň, že měsíc obíhá kolem planety. K jeho rychlosti musíš tedy rychlost planety přičíst. Na počátku tedy nepoletí 1082 m/s, ale o 31000 m/s víc. Je to o správné volbě vztažné soustavy.

No a to mi právě nejde :/ Rychlost jsem změnila, ale mám zřejmě špatně zvolené některé souřadnice...

Pouzil bych matematickou rotaci pres sin, cos. S vektory bych se nedrbal. Rychlost rotace pak urcuje narust uhlu. A udelal bych is na to objekty. Funkci pak predal jmeno objektu

obj.uhel += 0.2;
obj.x = floor (obj.x0 + obj.rx * cos(obj.uhel))
obj.y = floor (obj.y0 + obj.ry * sin(obj.uhel))

spocitej (slunce);
kresli (slunce);
spocitej (zeme);
kresli (zeme);
spocitej (mesic);
kresli (mesic);

#23 peter
Já jsem to původně dělala přes vykreslení elipsy pomocí Keplerovy rovnice, ale bylo mi to vráceno zpět, že mám použít gravitační zákon...

Mno, prepsal jsem to do js. Ale tez tam mam nekde chybu :) Pro x to pocita nejaky straslivy cisla nez pro y a ujizdi na stranu. (Kit - se nesmej!)
Ty tam mas nejake nepresnosti v pocitani sily, umocneni 1.5 se mi nezda, kdyz ma vzorecek bud r na druhou nebo r na druhou * r/r. A pak uz je treba pocitat aktualni polohu a ne nulovy stav (to ti psal Kit s tou vztaznou soustavou). 

<html>
<style>
#graph {border:1px solid #f00; width:600px; height:500px; position:relative;}
#graph div {position:absolute;}
</style>
<div id="graph">
  <div id="o0"><b>o</b></div>
  <div id="o1">o</div>
  <div id="o2">o</div>
</div>
<div id="label">
  <div id="l0"></div>
  <div id="l1"></div>
  <div id="l2"></div>
  <div>
    <input type=button value=start onclick="if (timer!=null) {clearInterval(timer);} timer = setInterval(rotace,k.timer);">
    <input type=button value=stop onclick="clearInterval(timer);">
  </div>
</div>

<script>
function kresli(obj)
{
obj.label.innerHTML = obj.name + ": x = %1, y = %2".replace('%1',obj.x).replace('%2',obj.y);
obj.el.style.left   = k.x0 + Math.floor(obj.x * k.zoomx) + 'px';
obj.el.style.top    = k.y0 + Math.floor(obj.y * k.zoomy) + 'px';
}

function rotace()
{
var r2, F, ax, ay, vx, vy, x,y;
// zeme
r2 = Math.pow(slunce.x - zeme.x, 2) + Math.pow(slunce.y - zeme.y, 2);
F = - k.kappa * slunce.m * zeme.m / r2;
ax = F * zeme.x;
ay = F * zeme.y;
vx = zeme.vx + ax * k.dt;	// v = vo + a * t
vy = zeme.vy + ay * k.dt;
x  = zeme.x + vx * k.dt;
y  = zeme.y + vy * k.dt;
zeme.vx = vx;
zeme.vy = vy;
zeme.x = x;
zeme.y = y;

/*
// mesic
r = 
F = - k.kappa * slunce.m * zeme.m / (r * r);
ax = F * mesic.x0;
ay = F * mesic.y0;
vx = mesic.vx0 + ax * k.dt;
vy = mesic.vy0 + ay * k.dt;
x  = mesic.x0 + vx * k.dt;
y  = mesic.y0 + vy * k.dt;
mesic.x = x;
mesic.y = y;
*/

kresli(slunce);
kresli(zeme);
kresli(mesic);
}

var k = {
	// nasa units: http://nssdc.gsfc.nasa.gov/…t_notes.html
	// nasa earth: http://nssdc.gsfc.nasa.gov/…rthfact.html
	// nasa sun:   http://nssdc.gsfc.nasa.gov/…sunfact.html
	zoomx: 1e-41,	// x * zoom, y * zoom
	zoomy: 1e-8,	// x * zoom, y * zoom
	x0:   300,	//px
	y0:   200,	//px
	timer: 300,	//ms
	
        kappa: 6.67384E-11,	// Gravitational Constant [m3*kg-1*s-2] (N*m2*kg-2) - nasaU
        AU:    1.49597900E+11,	// AU [m] - nasaU
	s_m: 1.988500E+30,	// mass [kg] - nasaS
	z_m: 5.9726E+24,	// mass [kg] - nasaE
        z_rp: 1.4709E+11,	// peri [m] = 0.98329 * AU - nasaE
	z_vmax: -30.29E+3,	// rychlost v perihelu [m/s] - nasaE
        dt: 3600		// time step [s]
}

var slunce = {
        name:  'Slunce',
        el:    document.getElementById('o0'),
        label: document.getElementById('l0'),
	x: 0,
	y: 0,
	m: k.s_m
}
var zeme = {
        name:  'Zeme',
        el:    document.getElementById('o1'),
        label: document.getElementById('l1'),
	x:  k.z_rp,	//x0
	y:  0,		//y0
        vx: 0,
        vy: k.z_vmax,
	m:  k.z_m,
	x0: k.z_rp,
	y0: 0,
        vx0: 0,
        vy0: k.z_vmax
}

var mesic = {
        name:  'Mesic',
        el:    document.getElementById('o2'),
        label: document.getElementById('l2'),
	x:  k.z_rp,		//x0
	y:  0,		//y0
        vx: 0,
        vy: k.z_vmax,
	m:  k.m_m,
	x0: k.z_rp,
	y0: 0,
        vx0: 0,
        vy0: k.z_vmax
}

kresli(slunce);
kresli(zeme);
kresli(mesic);


var timer = null;
</script>
</html>

Vcera jsem se uz k premysleni moc nedostal. V podstate bych si s tim tak nehral. Jen mne zajima, jestli fyzikalni nesmysly dokazi prepsat do rozumne programove podoby. Vzdycky mne na fyzikalnich a matem. knizkach stvou jejich neuplne postupy, ktere se bez znalosti strasne tezko prepisuji.
Mno, zatim jsem prisel na tohle:
- slunce jako stred, poloha x=0, y=0
- poloha zeme v perihelu je x=-1.47E+11 [m], y=0, rychlost vx=0, vy=30.29E+3 (zaporne x, protoze je vlevo od slunce, kladne vy, protoze se otaci zleva doprava... aspon podle toho, co jsem vystoural v google na obrazku, takze jestli je real jinak, tak to treba zmenit)
- Fg je teda kladna, jak to mam k.kappa * slunce.m * zeme.m / r2
- Fg je treba rozlozit na x,y podle uhlu v jakem pusobi. uhel je pres arkus tangent. pomer bude poloha zeme vuci slunci, rozdil x,y
alpha = Math.atan(zeme.x / zeme.y)
Fx = Fg * Math.cos(alpha);
Fy = Fg * Math.sin(alpha);
a z toho teprve muzes pocitat ax, ay
- mno, a dal jsem se zatim nedostal, takze je to na tobe. Ten vypocet ax jako F * s se mi nezda. Zatim se mi to googlem nepodarilo najit.
Fyzika je mi nekolik let vzdalena, takze uz o ni zas vim houbec a vsechno googluju :)

#26 peter
Ahoj, já fyziku a matiku studuju - teda učitelství, ale i tak fyziku. Tou fyzikální částí počítání gravitační síly a zrychlení jsem si jistá, dokonce jsem sem už výše dávala text, podle kterého jsem zpočátku pracovala. V tom chyba není.. chápu, že vám není jasná ta mocnina, ale je to prostě tak, protože síla je F=m*a. Hmotnost se zkrátí a a je potom F*složka x nebo y. Proto je v gravitačním zákoně vzdálenost r na 3... podle pythagorovy věty je počítaná jako odmocnina a tak to je ve výsledku druhá odmocnina na třetí - takže r na 3/2=1,5.

Posun os je to, co mi v tom dělá bordel. Střed soustavy můžu nechat [0,0], ale neuvidím nic, protože [0,0] je v levém horním rohu okna.

Co se týče perihelu, mám ho počítaný přes AU (astronomickou jednotku) a rychlost mám teda ne zcela přesně, uznávám.

Oběh planety funguje! Takže fyzikálně je to správně... ten Měsíc už ne

#27 Anetace
To se ale řeší velmí snadno a přehledně (posunutí středu).

Napiš si svoje podprogramy pro zobrazování, posílej do nich vypočítané údaje a teprve z nich vykresluj do skutečných souřadnic na obrazovce s posunutým středem...

Procedure Circle(x,y,Radius:integer)
begin
  Graph.Circle(x+stredX,y+stredY,Radius);
end

Aha. Fyzik, co se placa s programovani. To znam, od kolegyne :) Pro nas jasna vec, ze si na kazdou cas napiseme vlastni proceduru / funkci, ktera program krasne zjednodusi.

Ten prepocet se dela, jak pise Jodik. Jen se obvykle predavaji funkci cele objekty a ne jen data. V tom mem programu je to radek
k.x0 + Math.floor(obj.x * k.zoomx)
k.x0 - konstanta pro posun stredu v px (pixelech obrazovky)
floor - zaokrouhlovani dolu
obj.x - souradnice x, spocitana fyzikou
k.zoom - nasobitel pro prepocet fyzikalniho cisla na pixely
Cili, proti jodikovi zaokrouhluji na pixely az pri vykresleni.
x + stredX (JoDiK)
On by to musel delat jeste pred funkci. Ve tvem to prave mas bez toho posunuti
round(x * 2E-09) (Anetace)
200 + round(x * 2E-09) - treba takto

mno, a problem teda je, ze mi tve vypocty nechteli fungovat. Ale mozna jsem tam mel chybu. Neresim :)
 

Jo, proc ti vnucujeme ty funkce/procedury. To je kvuli chybe pri prepisovani. snadno se pri opakovani kodu muzes prehlednout. A pak tam mas treba jine prepocty pro zoom 2 - 2.2 - 3. Coz je teda fyzikalne-matematicka chyba, duvod pochybovat i o dalsich vypoctech.

m-round(x0*2E-09), n-round(y0*2E-09), Green);
        (*PutPixel (m-round(x0_mesic*2.2E-09), n-round(y0_mesic*2.2E-09), Red);*)

        WriteLn('x_0 :',round(x0*3E-09));

Tak mi to nedalo a taky jsem se o něco pokusil. Psal jsem to v Delphi 7, takže to sice není čistý Pascal, ale popravdě, kdo v něm dneska ještě seriózně dělá? Na inspiraci to bude snad dobrý.

Podstatná část:

unit OrbitalSystem;

interface

const
  G             = 6.673e-11;  (* gravitational constant [N(m/kg)^2] *) 
  ReducedG      = 6.673;      (* by 10^-11 *) 
  def_TimeStep  = 60;         (* [s] *)
  cTrailSteps   = 100;
  def_TrailTime = 1000000;    (* [s] *)

type
  TVector3e = Array[0..2] of Extended;

  TVector3er = packed record
    X:  Extended;
    Y:  Extended;
    Z:  Extended;
  end;

  TTrailArray = Array[0..Pred(cTrailSteps)] of TVector3e;

const
  ZeroVector3e: TVector3e = (0.0,0.0,0.0);
  
//==============================================================================

type
  TBody = class(TObject)
  private
    fCenterBody:        TBody;
    fIdentifier:        TGUID;
    fName:              String;
    fReducedValues:     Boolean;
    fMass:              Extended;
    fPosition:          TVector3e;
    fVelocity:          TVector3e;
    fAcceleration:      TVector3e;
    fElapsedTime:       Int64;
    fData:              Pointer;
    fInteractingBodies: Array of TBody;
    fStoreTrail:        Boolean;
    fTrailTime:         Int64;
    fTrailStep:         Int64;
    fTrailCount:        Integer;
    fTrailLast:         Int64;
    fTrail:             TTrailArray;
    Function GetInteractingBody(Index: Integer): TBody;
    Function GetInteractingBodiesCount: Integer;
    procedure SetStoreTrail(Value: Boolean);
    procedure SetTrailTime(Value: Int64);
    Function GetTrail(Index: Integer): TVector3e;
  protected
    Function IndexOfBody(Identifier: TGUID): Integer; virtual;
  public
    constructor Create(CenterBody: TBody; Mass: Extended; InitialPosition, InitialVelocity: TVector3e);
    destructor Destroy; override;
    procedure AddInteractingBody(Body: TBody); virtual;
    procedure RemoveInteractingBody(Body: TBody); virtual;
    procedure ClearInteractingBodies; virtual;
    procedure CalculateVectors(TimeDelta: LongWord = def_TimeStep); virtual;
    procedure MoveByVector(TimeDelta: LongWord = def_TimeStep); virtual;
    procedure MakeTrail; virtual;
    property CenterBody: TBody read fCenterBody;
    property Identifier: TGUID read fIdentifier;
    property Name: String read fName write fName;
    property ReducedValues: Boolean read fReducedValues write fReducedValues;
    property Mass: Extended read fMass write fMass;
    property Position: TVector3e read fPosition write fPosition;
    property PositionX: Extended read fPosition[0] write fPosition[0];
    property PositionY: Extended read fPosition[1] write fPosition[1];
    property PositionZ: Extended read fPosition[2] write fPosition[2];
    property Velocity: TVector3e read fVelocity write fVelocity;
    property VelocityX: Extended read fVelocity[0] write fVelocity[0];
    property VelocityY: Extended read fVelocity[1] write fVelocity[1];
    property VelocityZ: Extended read fVelocity[2] write fVelocity[2];
    property Acceleration: TVector3e read fAcceleration write fAcceleration;
    property AccelerationX: Extended read fAcceleration[0] write fAcceleration[0];
    property AccelerationY: Extended read fAcceleration[1] write fAcceleration[1];
    property AccelerationZ: Extended read fAcceleration[2] write fAcceleration[2];
    property ElapsedTime: Int64 read fElapsedTime;
    property InteractingBodies[Index: Integer]: TBody read GetInteractingBody;
    property InteractingBodiesCount: Integer read GetInteractingBodiesCount;
    property Data: Pointer read fData write fData;
    property StoreTrail: Boolean read fStoreTrail write SetStoreTrail;
    property TrailTime: Int64 read fTrailTime write SetTrailTime;
    property TrailStep: Int64 read fTrailStep;
    property TrailCount: Integer read fTrailCount;
    property Trail[Index: Integer]: TVector3e read GetTrail;
  end;
  
//==============================================================================  

  TBodyEvent = procedure(Sender: TObject; Body: TBody) of object;

  TOrbitalSystem = class(TObject)
  private
    fElapsedTime:   Int64;
    fBodies:        Array of TBody;
    fOnBodyCreate:  TBodyEvent;
    fOnBodyDestroy: TBodyEvent;
    Function GetBody(Index: Integer): TBody;
    Function GetBodiesCount: Integer;
  protected
    Function CheckDistances(Body: TBody): Boolean; virtual;
    procedure AddInteractingBodyToSystem(Body: TBody); virtual;
  public
    constructor Create;
    destructor Destroy; override;
    Function IndexOfBody(Body: TBody): Integer; overload; virtual;
    Function IndexOfBody(BodyIdentifier: TGUID): Integer; overload; virtual;
    Function IndexOfBody(BodyName: String): Integer; overload; virtual;
    Function AddBody(CenterBody: TBody; Mass: Extended; InitialPosition, InitialVelocity: TVector3e): TBody; overload; virtual;
    Function AddBody(CenterBodyIdx: Integer; Mass: Extended; InitialPosition, InitialVelocity: TVector3e): TBody; overload; virtual;
    Function RemoveBody(Body: TBody): Integer; virtual;
    procedure DeleteBody(Index: Integer); virtual;
    procedure ClearBodies; virtual;
    procedure ProcessStep(TimeDelta: LongWord = def_TimeStep); virtual;
    procedure SetReducedMode(Active: Boolean);
    property ElapsedTime: Int64 read fElapsedTime;
    property Bodies[Index: Integer]: TBody read GetBody;
    property BodiesCount: Integer read GetBodiesCount;
    property OnBodyCreate: TBodyEvent read fOnBodyCreate write fOnBodyCreate;
    property OnBodyDestroy: TBodyEvent read fOnBodyDestroy write fOnBodyDestroy;
  end;

  Function Vector3e(X,Y,Z: Extended): TVector3e;
  Function Vector3eToStr(Vector: TVector3e; Rounding: Integer = MAXINT): String;
  Function AddVectors(Vec1,Vec2: TVector3e): TVector3e;
  Function SubstractVectors(Vec1,Vec2: TVector3e): TVector3e;
  Function VectorScalarMultiply(Vector: TVector3e; Scalar: Extended): TVector3e;
  Function PointsDistance(Pt1,Pt2: TVector3e): Extended;
  Function VectorMagnitude(Vector: TVector3e): Extended;
  Function OppositeVector(Vector: TVector3e): TVector3e;

implementation

uses
  SysUtils, Math;

Function Vector3e(X,Y,Z: Extended): TVector3e;
begin
Result[0] := X;
Result[1] := Y;
Result[2] := Z;
end;

//------------------------------------------------------------------------------

Function Vector3eToStr(Vector: TVector3e; Rounding: Integer = MAXINT): String;
begin
If Rounding < MAXINT then
  Result := '[X: ' + FloatToStr(RoundTo(Vector[0],Rounding)) +
           '; Y: ' + FloatToStr(RoundTo(Vector[1],Rounding)) +
           '; Z: ' + FloatToStr(RoundTo(Vector[2],Rounding)) + ']'
else
  Result := '[X: ' + FloatToStr(Vector[0]) +
           '; Y: ' + FloatToStr(Vector[1]) +
           '; Z: ' + FloatToStr(Vector[2]) + ']';
end;

//------------------------------------------------------------------------------

Function AddVectors(Vec1,Vec2: TVector3e): TVector3e;
begin
Result[0] := Vec1[0] + Vec2[0];
Result[1] := Vec1[1] + Vec2[1];
Result[2] := Vec1[2] + Vec2[2];
end;

Function SubstractVectors(Vec1,Vec2: TVector3e): TVector3e;
begin
Result[0] := Vec1[0] - Vec2[0];
Result[1] := Vec1[1] - Vec2[1];
Result[2] := Vec1[2] - Vec2[2];
end;

//------------------------------------------------------------------------------

Function VectorScalarMultiply(Vector: TVector3e; Scalar: Extended): TVector3e;
begin
Result[0] := Vector[0] * Scalar;
Result[1] := Vector[1] * Scalar;
Result[2] := Vector[2] * Scalar;
end;

//------------------------------------------------------------------------------

Function PointsDistance(Pt1,Pt2: TVector3e): Extended;
begin
Result := Sqrt(Sqr(Pt1[0] - Pt2[0]) + Sqr(Pt1[1] - Pt2[1]) + Sqr(Pt1[2] - Pt2[2]));
end;

//------------------------------------------------------------------------------

Function VectorMagnitude(Vector: TVector3e): Extended;
begin
Result := Sqrt(Sqr(Vector[0]) + Sqr(Vector[1]) + Sqr(Vector[2]));
end;

//------------------------------------------------------------------------------

Function OppositeVector(Vector: TVector3e): TVector3e;
begin
Result[0] := -Vector[0];
Result[1] := -Vector[1];
Result[2] := -Vector[2];
end;

//==============================================================================
//******************************************************************************
//==============================================================================

Function TBody.GetInteractingBody(Index: Integer): TBody;
begin
If (Index >= Low(fInteractingBodies)) and (Index <= High(fInteractingBodies)) then
  Result := fInteractingBodies[Index]
else raise Exception.Create('TBody.GetInteractingBody: Index (' + IntToStr(Index) + ') out of bounds.');
end;

//------------------------------------------------------------------------------

Function TBody.GetInteractingBodiesCount: Integer;
begin
Result := Length(fInteractingBodies);
end;

//------------------------------------------------------------------------------

procedure TBody.SetStoreTrail(Value: Boolean);
begin
If Value <> fStoreTrail then
  begin
    If Value then
      begin
        fTrailCount := 0;
        fTrailLast := -fTrailStep + 1;
      end
    else fTrailCount := 0;
    fStoreTrail := Value;
  end;
end;

//------------------------------------------------------------------------------

procedure TBody.SetTrailTime(Value: Int64);
begin
fTrailCount := 0;
fTrailTime := Value;
fTrailStep := Round(Value / cTrailSteps);
fTrailLast := -fTrailStep + 1;
end;

//------------------------------------------------------------------------------

Function TBody.GetTrail(Index: Integer): TVector3e;
begin
If (Index >= 0) and (Index < fTrailCount) then
  Result := fTrail[Index]
else raise Exception.Create('TBody.GetTrail: Index (' + IntToStr(Index) + ') out of bounds.');
end;

//==============================================================================

Function TBody.IndexOfBody(Identifier: TGUID): Integer;
begin
If Length(fInteractingBodies) > 0 then
  For Result := Low(fInteractingBodies) to High(fInteractingBodies) do
    If IsEqualGUID(Identifier,fInteractingBodies[Result].Identifier) then Exit;
Result := -1;
end;

//==============================================================================

constructor TBody.Create(CenterBody: TBody; Mass: Extended; InitialPosition, InitialVelocity: TVector3e);
begin
inherited Create;
fCenterBody := CenterBody;
CreateGUID(fIdentifier);
fName := GUIDToString(fIdentifier);
fReducedValues := True;
fMass := Mass;
If Assigned(CenterBody) then
  begin
    fPosition := AddVectors(InitialPosition,CenterBody.Position);
    fVelocity := AddVectors(InitialVelocity,CenterBody.Velocity);
  end
else
  begin
    fPosition := InitialPosition;
    fVelocity := InitialVelocity;
  end;
fAcceleration := ZeroVector3e;
fElapsedTime := 0;
SetLength(fInteractingBodies,0);
fStoreTrail := True;
SetTrailTime(def_TrailTime);
fTrailCount := 0;
fTrailLast := -fTrailStep + 1;
end;

//------------------------------------------------------------------------------

destructor TBody.Destroy;
begin
ClearInteractingBodies;
inherited;
end;

//------------------------------------------------------------------------------

procedure TBody.AddInteractingBody(Body: TBody);
begin
If not IsEqualGUID(Body.Identifier,Self.Identifier) and (IndexOfBody(Body.Identifier) < 0) then
  begin
    SetLength(fInteractingBodies,Length(fInteractingBodies) + 1);
    fInteractingBodies[High(fInteractingBodies)] := Body;
  end;
end;

//------------------------------------------------------------------------------

procedure TBody.RemoveInteractingBody(Body: TBody);
var
  Index,i:  Integer;
begin
If Body = fCenterBody then fCenterBody := nil;
Index := IndexOfBody(Body.Identifier);
If Index >= 0 then
  begin
    For i := Index to Pred(High(fInteractingBodies)) do
      fInteractingBodies[i] := fInteractingBodies[i + 1];
    SetLength(fInteractingBodies,Length(fInteractingBodies) - 1)
  end;
end;

//------------------------------------------------------------------------------

procedure TBody.ClearInteractingBodies;
begin
SetLength(fInteractingBodies,0);
end;

//------------------------------------------------------------------------------

procedure TBody.CalculateVectors(TimeDelta: LongWord = def_TimeStep);
(*
  Ai = (G * Mj * (Pj - Pi)) / (Rij ^ 3)

  Ai  - resulting acceleration vector
  G   - gravitational constant
  Mj  - mass of interacting body
  Pj  - position vector of interacting body
  Pi  - position vector of this body
  Rij - distance between both bodies
*)
var
  TempAccVector:  TVector3e;
  i,j:            Integer;
  Distance:       Extended;
  BigG:           Extended;
begin
If ReducedValues then BigG := ReducedG
  else BigG := G;
fAcceleration := ZeroVector3e;
For i := Low(fInteractingBodies) to High(fInteractingBodies) do
  begin
    Distance := PointsDistance(fInteractingBodies[i].Position,Self.Position);
    TempAccVector := ZeroVector3e;    
    For j := Low(TVector3e) to High(TVector3e) do
      TempAccVector[j] := (BigG * fInteractingBodies[i].Mass * (fInteractingBodies[i].Position[j] - Self.Position[j])) / Power(Distance,3);
    fAcceleration := AddVectors(fAcceleration,TempAccVector);
  end;
fVelocity := AddVectors(fVelocity,VectorScalarMultiply(fAcceleration,TimeDelta));
end;

//------------------------------------------------------------------------------

procedure TBody.MoveByVector(TimeDelta: LongWord = def_TimeStep);
begin
fPosition := AddVectors(fPosition,VectorScalarMultiply(fVelocity,TimeDelta));
fElapsedTime := fElapsedTime + TimeDelta;
end;

//------------------------------------------------------------------------------

procedure TBody.MakeTrail;
begin
If fStoreTrail then
  If (fElapsedTime - fTrailLast) > fTrailStep then
    begin
      If fTrailCount >= cTrailSteps then
        begin
          fTrailCount := cTrailSteps;
          Move(fTrail[1],fTrail[0],(cTrailSteps - 1) * SizeOf(TVector3e));
        end
      else Inc(fTrailCount);
      If Assigned(fCenterBody) then
        fTrail[Pred(fTrailCount)] := SubstractVectors(fPosition,fCenterBody.Position)
      else
        fTrail[Pred(fTrailCount)] := fPosition;
      fTrailLast := fElapsedTime;
    end;
end;

//==============================================================================
//******************************************************************************
//==============================================================================

Function TOrbitalSystem.GetBody(Index: Integer): TBody;
begin
If Length(fBodies) > 0 then
  begin
    If (Index >= Low(fBodies)) and (Index <= High(fBodies)) then
      Result := fBodies[Index]
    else
      raise Exception.Create('TOrbitalSystem.GetBody: Index (' + IntToStr(Index) + ') out of bounds.');
  end
else Result := nil;
end;

//------------------------------------------------------------------------------

Function TOrbitalSystem.GetBodiesCount: Integer;
begin
Result := Length(fBodies);
end;

//==============================================================================

Function TOrbitalSystem.CheckDistances(Body: TBody): Boolean;
var
  i:  Integer;
begin
Result := True;
For i := Low(fBodies) to High(fBodies) do
  If PointsDistance(fBodies[i].Position,Body.Position) = 0 then
    begin
      Result := False;
      Break;
    end;
end;

//------------------------------------------------------------------------------

procedure TOrbitalSystem.AddInteractingBodyToSystem(Body: TBody);
var
  i:  Integer;
begin
For i := Low(fBodies) to High(fBodies) do
  begin
    fBodies[i].AddInteractingBody(Body);
    Body.AddInteractingBody(fBodies[i]);
  end;
end;

//==============================================================================

constructor TOrbitalSystem.Create;
begin
inherited Create;
fElapsedTime := 0;
SetLength(fBodies,0);
end;

//------------------------------------------------------------------------------

destructor TOrbitalSystem.Destroy;
begin
ClearBodies;
inherited;
end;

//------------------------------------------------------------------------------

Function TOrbitalSystem.IndexOfBody(Body: TBody): Integer;
begin
If Length(fBodies) > 0 then
  For Result := Low(fBodies) to High(fBodies) do
    If fBodies[Result] = Body then Exit;
Result := -1;
end;

//---                                                                        ---

Function TOrbitalSystem.IndexOfBody(BodyIdentifier: TGUID): Integer;
begin
If Length(fBodies) > 0 then
  For Result := Low(fBodies) to High(fBodies) do
    If IsEqualGUID(fBodies[Result].Identifier, BodyIdentifier) then Exit;
Result := -1;
end;

//---                                                                        ---

Function TOrbitalSystem.IndexOfBody(BodyName: String): Integer;
begin
If Length(fBodies) > 0 then
  For Result := Low(fBodies) to High(fBodies) do
    If AnsiSameText(fBodies[Result].Name, BodyName) then Exit;
Result := -1;
end;

//------------------------------------------------------------------------------

Function TOrbitalSystem.AddBody(CenterBody: TBody; Mass: Extended; InitialPosition, InitialVelocity: TVector3e): TBody;
begin
Result := nil;
If Assigned(CenterBody) and (IndexOfBody(CenterBody) < 0) then
  raise Exception.Create('TOrbitalSystem.AddBody: Passed CenterBody is not in list of system bodies.');
Result := TBody.Create(CenterBody,Mass,InitialPosition,InitialVelocity);
If not CheckDistances(Result) then
  begin
    FreeAndNil(Result);
    raise Exception.Create('TOrbitalSystem.AddBody: Detected zero distance to other body.');
  end;
SetLength(fBodies,Length(fBodies) + 1);
fBodies[High(fBodies)] := Result;
AddInteractingBodyToSystem(Result);
If Assigned(fOnBodyCreate) then fOnBodyCreate(Self,Result); 
end;

//---                                                                        ---

Function TOrbitalSystem.AddBody(CenterBodyIdx: Integer; Mass: Extended; InitialPosition, InitialVelocity: TVector3e): TBody;
begin
If (Length(fBodies) > 0) and (CenterBodyIdx >= Low(fBodies)) and (CenterBodyIdx <= High(fBodies)) then
  Result := AddBody(fBodies[CenterBodyIdx],Mass,InitialPosition,InitialVelocity)
else
  Result := AddBody(nil,Mass,InitialPosition,InitialVelocity);
end;

//------------------------------------------------------------------------------

Function TOrbitalSystem.RemoveBody(Body: TBody): Integer;
begin
Result := IndexOfBody(Body);
If Result >= 0 then DeleteBody(Result);
end;

//------------------------------------------------------------------------------

procedure TOrbitalSystem.DeleteBody(Index: Integer);
var
  i:  Integer;
begin
If Length(fBodies) > 0 then
  begin
    If (Index >= Low(fBodies)) and (Index <= High(fBodies)) then
      begin
        For i := Low(fBodies) to High(fBodies) do
          fBodies[i].RemoveInteractingBody(fBodies[Index]);
        If Assigned(fOnBodyDestroy) then fOnBodyDestroy(Self,fBodies[Index]);
        fBodies[Index].Free;
        For i := Index to Pred(High(fBodies)) do
          fBodies[i] := fBodies[i + 1];
        SetLength(fBodies,Length(fBodies) - 1)
      end
    else
      raise Exception.Create('TOrbitalSystem.DeleteBody: Index (' + IntToStr(Index) + ') out of bounds.');
  end;
end;

//------------------------------------------------------------------------------

procedure TOrbitalSystem.ClearBodies;
var
  i:  Integer;
begin
For i := Low(fBodies) to High(fBodies) do
  begin
    If Assigned(fOnBodyDestroy) then fOnBodyDestroy(Self,fBodies[i]);
    fBodies[i].Free;
  end;
fElapsedTime := 0;
SetLength(fBodies,0);
end;

//------------------------------------------------------------------------------

procedure TOrbitalSystem.ProcessStep(TimeDelta: LongWord = def_TimeStep);
var
  i:  Integer;
begin
For i := Low(fBodies) to High(fBodies) do
  fBodies[i].CalculateVectors(TimeDelta);
For i := Low(fBodies) to High(fBodies) do
  fBodies[i].MoveByVector(TimeDelta);
For i := Low(fBodies) to High(fBodies) do
  fBodies[i].MakeTrail;  
fElapsedTime := fElapsedTime + TimeDelta;
end;

//------------------------------------------------------------------------------

procedure TOrbitalSystem.SetReducedMode(Active: Boolean);
var
  i:  Integer;
begin
For i := Low(fBodies) to High(fBodies) do
  fBodies[i].ReducedValues := Active;
end;   

end.

Funguje to bez problémů na jakýkoliv počet těles. Ale moc jsem neřešil přesnost, takže pokud se to bude překládat v něčem, co 80bit Extended nahradí za 64bit Double, tak by mohl být problém. A taky pozor na výpočetní náročnost, roste exponenciálně s každým dalším tělesem (výkonovou optimalizaci jsem nedělal a dělat nebudu).

A tady je komplet program kde je vše včetně kreslení s proměnlivým zoomem atd.(není moc otestovaný ale snad poběží): http://uloz.to/…imulator-zip
 

Problem delphi a pascalu je v tom, ze na to potrebujes extra program. Ale webovy prohlizec ma dneska kazdy :) Takze fakt nechapu, jak nekdo jeste ve skole muze do studentu hustit pascal/freepascal.

Díky všem, kteří se zapojili! S pomocí spousta lidí, kteří mě každý posunul trochu jiným směrem jsem program dodělala, včera odevzdala a udělala zkoušku. S Pascalem už nechci mít nic společnýho a mám ho z krku!

Sem to prepsal konecne asi spravne podle puvodnich rovnic. Konecne to rotuje. Ale mesic nejak podivne, tak se moc nesmejte, precijen nejsem fyzik :) (navic, jestli to spravne chapu, tak je to jen rovnice pro rotaci kolem orbity bez uvazovani okolnich teles; To uz by byl lepsi klasicky sin, cos, ne?) Co tam mam spatne? 

<html>
<style>
#graph {border:1px solid #f00; width:600px; height:500px; position:relative;}
#graph div {position:absolute;}
</style>
<div id="graph">
  <div id="o0"><b>o</b></div>
  <div id="o1">o</div>
  <div id="o2">o</div>
</div>
<div id="label">
  <div id="l0"></div>
  <div id="l1"></div>
  <div id="l2"></div>
  <div>
    <input type=button value=start onclick="if (timer!=null) {clearInterval(timer);} timer = setInterval(rotace,k.timer);">
    <input type=button value=stop onclick="clearInterval(timer);">
  </div>
</div>

<script>
function vector(x,y)
{
this.x = x;
this.y = y;
return this;
}

function kresli(orbit,center)
{
orbit.ppx.x = center.ppx.x + Math.round(orbit.p.x * orbit.zoom.x);
orbit.ppx.y = center.ppx.y + Math.round(orbit.p.y * orbit.zoom.y);
orbit.label.innerHTML = orbit.name + ": x = %1, y = %2".replace('%1',orbit.p.x).replace('%2',orbit.p.y);
orbit.el.style.left   = orbit.ppx.x + 'px';
orbit.el.style.top    = orbit.ppx.y + 'px';
}

function spocitej(orbit,center)
{
var F,a,v,p,rr;
a = new vector();
v = new vector();
p = new vector();
rr  = Math.pow(Math.pow(center.p.x - orbit.p.x, 2) + Math.pow(center.p.y - orbit.p.y, 2), 1.5);
F   = - k.kappa * center.m / rr;
a.x = orbit.p.x * F;
a.y = orbit.p.y * F;
v.x = a.x * k.dt + orbit.v.x;
v.y = a.y * k.dt + orbit.v.y;
p.x = v.x * k.dt + orbit.p.x;
p.y = v.y * k.dt + orbit.p.y;
orbit.v = v;
orbit.p = p;
}

function rotace()
{
spocitej(zeme , slunce);
spocitej(mesic, zeme);
kresli(zeme  , slunce);
kresli(mesic , zeme);
}

var k = {
	// nasa units: http://nssdc.gsfc.nasa.gov/…t_notes.html
	// nasa earth: http://nssdc.gsfc.nasa.gov/…rthfact.html
	// nasa sun:   http://nssdc.gsfc.nasa.gov/…sunfact.html
	// nasa moon:  http://nssdc.gsfc.nasa.gov/…oonfact.html
	zoom : new vector(1e-9,1e-9),	// x * zoom, y * zoom
	timer: 100,			//program timer ms
	
        kappa : 6.67384E-11,	// Gravitational Constant [m3*kg-1*s-2] (N*m2*kg-2) - nasaU
        AU    : 1.49597900E+11,	// AU [m] - nasaU
	s_m   : 1.988500E+30,	// mass [kg] - nasaS
	z_m   : 5.9726E+24,	// mass [kg] - nasaE
        z_rp  : 1.4709E+11,	// peri [m] = 0.98329 * AU - nasaE
	z_vmax: 30.29E+3,	// rychlost v perihelu [m/s] - nasaE
	m_m   : 0.07342E+24,	// mass [kg] - nasaE
        m_rp  : 0.3633E+9,	// peri [m] = 0.98329 * AU - nasaE
	m_vmax: 1.076E+3,	// rychlost v perihelu [m/s] - nasaE
        dt    : 3600*24		// 24h time step [s]
}

var slunce = {
        name:  'Slunce',
        el:    document.getElementById('o0'),
        label: document.getElementById('l0'),
	m  : k.s_m,
	p  : new vector(0, 0),
	ppx: new vector(300, 200),
	zoom: k.zoom
	};
var zeme = {
        name : 'Zeme',
        el   : document.getElementById('o1'),
        label: document.getElementById('l1'),
	m  : k.z_m,
	p  : new vector(-k.z_rp, 0),
        v  : new vector(0, k.z_vmax),
	ppx: new vector(0, 0),
	zoom: k.zoom
	};
var mesic = {
        name:  'Mesic',
        el:    document.getElementById('o2'),
        label: document.getElementById('l2'),
	m  : k.m_m,
	p  : new vector(-k.m_rp, 0),
        v  : new vector(0, k.m_vmax),
	ppx: new vector(0, 0),
	zoom: new vector(k.zoom.x*10, k.zoom.y*10)
	}

kresli(slunce, slunce);
kresli(zeme  , slunce);
kresli(mesic , zeme);


var timer = null;
</script>
</html>

Kdyby ses podíval na to, co jsem postnul, a hned to nezavrhnul protože "pascal sucks, javascript ftw!" (mimochodem chtěl bych vidět jak uděláš v prohlížeči skutečnou simulaci), tak uvidíš, že žádná goniometrie není potřeba a ten můj postup počítá vliv všech simulovaných těles (viz http://en.wikipedia.org/…bit_modeling#…).

Co tam je za chybu nevim, ale koukám že neuchováváš vektor rychlosti ani zrychlení, což je jaksi potřeba, tudíž ty výpočty nemůžou dávat smysl.

Rychlost uchovavam, zrychleni nee, to se pocita pokazde z polohy.
Ja to prepsal do js, protoze je to pro mne jednodussi.
Nic proti tvemu reseni, ale stovky radku se mi nechtelo zkoumat. Ale tak treba to udelam a odhalim nejvetsi zahadu vesmiru, proc mi to nejde :)

Jo, spoustu radku bys usetril, kdybys pouzil canvas a nesnazil se to prekreslovat rucne..

Zrychlení není potřeba, to je pravda.
Nemusíš zkoumat celej kód, výpočty jsou kolem řádku 350 v tom uvedeným úryvku.
Překreslování - každej frame je jinej (nemění se toliko centrální kříž a osy) tak nevim kde se co dá podle tebe ušetřit. Jestli myslíš abych to kreslil přímo do canvasu formu, tak to nelze, neb to problikává (i při doublebuffered), a použít nějakou jinou komponentu - no to už to můžu nechat jak to je. Navíc takhle si můžu nechat "vyrenderovat" obrázek ve vyšším rozlišení a uložit ho.

Tak jsem se opět dlooouze podíval na ten tvůj JS kód a našel jsem ti tu chybu. Řádek 44 (rr = ...), používáš relativní souřadnice ale výpočet předpokládá absolutní. Abych to ujasnil, u soustavy Slunce - Země to nedělá problém protože Slunce je v souřadnicích [0,0], ale když počítáš Země - Měsíc tak najednou tam cpeš jako centrum Zemi na souřadnicích [nějaký miliardy, nějaký miliardy] ale měsíc máš na [nějaký statisíce, nějaký statisíce]. Proto ti z toho lezou bláboly. Nejjednodušší řešení je místo "center.p.x - orbit.p.x" dát prostě "- orbit.p.x" (máš relativní souřadnice takže můžeš vždy předpokládat centrum v bodě [0,0]) a bude to šlapat (totéž udělej pochopitelně u Y).

Jeee, parada, toci se, toci a ma modre oci. Udelal jsem jinou zmenu. Pridal jsem tam promenou 3.  A upravil F a funkce. A pak taky chybku na konci programu, kde mela byt zeme
1) function spocitej(orbit,orbit2,center)
2) F   = - k.kappa * orbit2.m / rr;
3) spocitej(zeme , slunce, slunce);
spocitej(mesic, zeme, slunce);
4) kresli(zeme  , slunce);
kresli(mesic , zeme);

Ted by to mohlo fungovat i se zmenou polohy slunce, snad :)

Dobre, tak to neni uplne ono, raci to opravim, jak pises.