program

%orbitARK - program untuk menggambarkan orbit planet
% dengan metoda Runge-Kutta adaptive
clear; help orbitARK
r0 = input('Masukkan posisi awal dalam AU - ');
v0 = input('Masukkan kecepatan awal dalam AU/thn - ');
h = input('Masukkan nilai time-step dalam thn - ');
r = [r0 0]; %vektor posisi
v = [0 v0]; %vektor kecepatan
GM = 4*pi^2;
m = 1;
time = 0;
stepmax = 200;
state = [r(1) r(2) v(1) v(2)];
err = 1.e-3;
for i=1:stepmax
rplot(i) = norm(r);
wplot(i) = atan2(r(2), r(1));
tplot(i) = time;
kinetik(i) = 0.5*m*norm(v)^2;
potensial(i) = -GM*m/norm(r);
[state, time, h] = ARK(state, time, h, err, 'gerak', GM);
r = [state(1) state(2)];
v = [state(3) state(4)];
time = time + h;
end
subplot(121)
polar(wplot, rplot, '+')
grid
ylabel('Posisi (AU)')
title('Gerak Orbit')
subplot(122)
Etotal = kinetik + potensial;
plot(tplot,kinetik,'-.', tplot,potensial,'--', tplot,Etotal,'-')
xlabel('Waktu (thn)')
ylabel('Energi')
title('EK(dot) EP(dash) ETot(solid)')
subplot(111)



%orbitE - program untuk menggambarkan orbit planet
% dengan metoda Euler
clear; help orbitE
r0 = input('Masukkan posisi awal dalam AU - ');
v0 = input('Masukkan kecepatan awal dalam AU/thn - ');
h = input('Masukkan nilai time-step dalam thn - ');
r = [r0 0]; %vektor posisi
v = [0 v0]; %vektor kecepatan
GM = 4*pi^2;
m = 1;
time = 0;
stepmax = 200;
for i=1:stepmax
rplot(i) = norm(r);
wplot(i) = atan2(r(2), r(1));
tplot(i) = time;
kinetik(i) = 0.5*m*norm(v)^2;
potensial(i) = -GM*m/norm(r);
a = -GM*r/norm(r)^3;
r = r + h*v;
v = v + h*a;
time = time + h;
end
subplot(121)
polar(wplot, rplot, '+')
grid
ylabel('Posisi (AU)')
title('Gerak Orbit')
subplot(122)
Etotal = kinetik + potensial;
plot(tplot,kinetik,'-.', tplot,potensial,'--', tplot,Etotal,'-')
xlabel('Waktu (thn)')
ylabel('Energi')
title('EK(dot) EP(dash) ETot(solid)')
subplot(111)

%orbitEC - program untuk menggambarkan orbit planet
% dengan metoda Euler-Cromer
clear; help orbitEC
r0 = input('Masukkan posisi awal dalam AU - ');
v0 = input('Masukkan kecepatan awal dalam AU/thn - ');
h = input('Masukkan nilai time-step dalam thn - ');
r = [r0 0]; %vektor posisi
v = [0 v0]; %vektor kecepatan
GM = 4*pi^2;
m = 1;
time = 0;
stepmax = 200;
for i=1:stepmax
rplot(i) = norm(r);
wplot(i) = atan2(r(2), r(1));
tplot(i) = time;
kinetik(i) = 0.5*m*norm(v)^2;
potensial(i) = -GM*m/norm(r);
a = -GM*r/norm(r)^3;
v = v + h*a;
r = r + h*v;
time = time + h;
end
subplot(121)
polar(wplot, rplot, '+')
grid
ylabel('Posisi (AU)')
title('Gerak Orbit')
subplot(122)
Etotal = kinetik + potensial;
plot(tplot,kinetik,'-.', tplot,potensial,'--', tplot,Etotal,'-')
xlabel('Waktu (thn)')
ylabel('Energi')
title('EK(dot) EP(dash) ETot(solid)')
subplot(111)

%orbitRK2 - program untuk menggambarkan orbit planet
% dengan metoda Runge-Kutta orde-2
clear; help orbitRK2
r0 = input('Masukkan posisi awal dalam AU - ');
v0 = input('Masukkan kecepatan awal dalam AU/thn - ');
h = input('Masukkan nilai time-step dalam thn - ');
r = [r0 0]; %vektor posisi
v = [0 v0]; %vektor kecepatan
GM = 4*pi^2;
m = 1;
time = 0;
stepmax = 200;
state = [r(1) r(2) v(1) v(2)];
for i=1:stepmax
rplot(i) = norm(r);
wplot(i) = atan2(r(2), r(1));
tplot(i) = time;
kinetik(i) = 0.5*m*norm(v)^2;
potensial(i) = -GM*m/norm(r);
state = RK2(state, time, h, 'gerak', GM);
r = [state(1) state(2)];
v = [state(3) state(4)];
time = time + h;
end
subplot(121)
polar(wplot, rplot, '+')
grid
ylabel('Posisi (AU)')
title('Gerak Orbit')
subplot(122)
Etotal = kinetik + potensial;
plot(tplot,kinetik,'-.', tplot,potensial,'--', tplot,Etotal,'-')
xlabel('Waktu (thn)')
ylabel('Energi')
title('EK(dot) EP(dash) ETot(solid)')
subplot(111)

%orbitRK4 - program untuk menggambarkan orbit planet
% dengan metoda Runge-Kutta orde-2
clear; help orbitRK4
r0 = input('Masukkan posisi awal dalam AU - ');
v0 = input('Masukkan kecepatan awal dalam AU/thn - ');
h = input('Masukkan nilai time-step dalam thn - ');
r = [r0 0]; %vektor posisi
v = [0 v0]; %vektor kecepatan
GM = 4*pi^2;
m = 1;
time = 0;
stepmax = 200;
state = [r(1) r(2) v(1) v(2)];
for i=1:stepmax
rplot(i) = norm(r);
wplot(i) = atan2(r(2), r(1));
tplot(i) = time;
kinetik(i) = 0.5*m*norm(v)^2;
potensial(i) = -GM*m/norm(r);
state = RK4(state, time, h, 'gerak', GM);
r = [state(1) state(2)];
v = [state(3) state(4)];
time = time + h;
end
subplot(121)
polar(wplot, rplot, '+')
grid
ylabel('Posisi (AU)')
title('Gerak Orbit')
subplot(122)
Etotal = kinetik + potensial;
plot(tplot,kinetik,'-.', tplot,potensial,'--', tplot,Etotal,'-')
xlabel('Waktu (thn)')
ylabel('Energi')
title('EK(dot) EP(dash) ETot(solid)')
subplot(111)