sciencetalk

Code: Selecteer alles

//Te kiezen;

p = 6.957*10^8; // Het perihelium (of periastron);
MG = 1.327*10^(20); //Het product van de massa van het hemellichaam en de gravitatieconstante;
d = 0.5; // De genomen tijdstappen bij de baanberekening;

//Natuurconstante;

c = 2.998*10^8; // De lichtsnelheid;

//De situatie wordt weergegeven door de baanvector (x,y,v_x,v_y);

//Hulpfuncties;

function ox=otx(z) //Onttrekt x aan de baanvector;
    ox=z([1])
endfunction

function oy=oty(z) //Onttrekt y aan de baanvector;
    oy=z([2])
endfunction

function ovx=otvx(z) //Onttrekt v_x aan de baanvector;
    ovx=z([3])
endfunction

function ovy=otvy(z) //Onttrekt v_y aan de baanvector;
    ovy=z([4])
endfunction

function tang=tangb(z) //Berekent de tangens van de bewegingshoek;
    tang=(otvy(z))/(otvx(z))
endfunction

function nvxx=nuvxx(z)  //Zet v_x op de plaats van x terug in nulvector;
    nvxx = otvx(z)*[1,0,0,0]
endfunction

function nvyy=nuvyy(z)  //Zet v_y op de plaats van y terug in nulvector;
    nvyy = otvy(z)*[0,1,0,0]
endfunction

function r=lengr(z) //Berekent de lengte van de plaatsvector;
    r = sqrt((otx(z))^2 + (oty(z))^2)
endfunction

function er=ver(z) //Berekent de plaatseenheidsrichtingsvector;
    er = [otx(z)/lengr(z) , oty(z)/lengr(z)]
endfunction

function et=vet(z) //Berekent de tangentiële eenheidsrichtingsvector;
    et = [(otvx(z))/c , (otvy(z))/c]
endfunction

function va=veca(z) //Berekent het vectordeel van de versnelling;
    va = ver(z) - (ver(z)*((vet(z))'))*vet(z)
endfunction

function ac=acc(z)
    ac = ((-MG)/((lengr(z))^2))*veca(z) //Berekent de versnelling;
endfunction

function oax=otax(z) //Onttrekt a_x aan de versnellingsvector;
    oax=z([1])
endfunction

function oay=otay(z) //Onttrekt a_y aan de versnellingsvector;
    oay=z([2])
endfunction

function naxvx=nuaxvx(z)  //Zet a_x op de plaats van v_x terug in nulvector;
    naxvx = otax(z)*[0,0,1,0]
endfunction

function nayvy=nuayvy(z)  //Zet a_y op de plaats van v_y terug in nulvector;
    nayvy = otay(z)*[0,0,0,1]
endfunction


//De baanfunctie berekent de baanvector voor  t_n uit de baanvector voor  t_n-1;

function baan = baanf(n)
    if n == 0 baan = [0,p,c,0]
    else
      baan = baanf(n-1) + nuvxx(baanf(n-1))*d + nuvyy(baanf(n-1))*d + nuaxvx(baanf(n-1))*d + nuayvy(baanf(n-1))*d
        end
endfunction



// Punten van de baan vanaf het perihelium (of periastron);


tangb(baanf(0))
otx(baanf(0))
oty(baanf(0))
otvx(baanf(0))
otvy(baanf(0))


tangb(baanf(1))
otx(baanf(1))
oty(baanf(1))
otvx(baanf(1))
otvy(baanf(1))


tangb(baanf(2))
otx(baanf(2))
oty(baanf(2))
otvx(baanf(2))
otvy(baanf(2))


tangb(baanf(3))
otx(baanf(3))
oty(baanf(3))
otvx(baanf(3))
otvy(baanf(3))


tangb(baanf(4))
otx(baanf(4))
oty(baanf(4))
otvx(baanf(4))
otvy(baanf(4))


tangb(baanf(4))
otx(baanf(4))
oty(baanf(4))
otvx(baanf(4))
otvy(baanf(4))


tangb(baanf(5))
otx(baanf(5))
oty(baanf(5))
otvx(baanf(5))
otvy(baanf(5))