orbit Program to compute the orbit of a comet MATLAB ro i

%orbit - Program to compute the orbit of a comet ( MATLAB)

ro = input(\'Enter initial radial distance (AU):\');

vo = input(\'Enter initial tangential velocity (AU/yr):\');

r = [ro 0];

v = [0 vo];

GM = 4*pi^2;

m = 1;

t = 0;

nstep = input(\'Enter number of steps:\');

tau = input(\'Enter time step (yr):\');

NM = menu(\' Method:\',\'Euler\',\'Euler-Cromer\',\'verlet\',\'RK4\', \'RKA\');

%NM= input(\'Enter 1 for Euler, 2 for Euler-Comer,3 for verlet,4 for RK4,or 5 for RKa\')

r_old = r - v*tau;

state = [r(1) r(2) v(1) v(2)];

err = 1.e-3;

for istep = 1:nstep

rplot(istep)=norm(r);

thplot(istep) = atan2(r(2),r(1));

tplot(istep) = t;

K(istep) = 0.5*m*norm(v)^2;

P(istep) = -GM*m/norm(r);

E = K + P;

if( NM == 1) % Euler step

a = -GM*r/norm(r)^3;

r = r + tau*v;

v = v + tau*a;

t = t + tau;

elseif( NM == 2) % Euler-Cromer step

a = -GM*r/norm(r)^3;

v = v + tau*a;

r = r + tau*v;

t = t + tau;

elseif( NM == 3) %Verlet method

a = -GM*r/norm(r)^3;

r_new = 2*r - r_old + a*tau^2;

v = (r_new - r_old)/(2*tau);

r_old = r;

r = r_new;

t = t + tau;

elseif(NM == 4) %RK4

state = rk4(state, t, tau, \'gravrk\', GM);

r = [state(1) state(2)];

v = [state(3) state(4)];

t = t + tau;

   elseif(NM == 5) %RKa

[state, t, tau] = rka(state, t, tau, err, \'gravrk\', GM);

r = [state(1) state(2)];

v = [state(3) state(4)];

end

end

polar(thplot,rplot,\'+\')

grid

pause

plot(tplot,E,tplot,K,tplot,P)

legend(\'E\', \'K\', \'P\');

Modi tee orbit program So so Tummimy fxed mumber time steps, ttee Modity toe orbit to Compute He period eccentricity axis

Solution

% orbit - Program to compute the orbit of a comet. clear all; help orbit; % Clear memory and print header %* Set initial position and velocity of the comet. r0 = input(\'Enter initial radial distance (AU): \'); v0 = input(\'Enter initial tangential velocity (AU/yr): \'); r = [r0 0]; v = [0 v0]; state = [ r(1) r(2) v(1) v(2) ]; % Used by R-K routines %* Set physical parameters (mass, G*M) GM = 4*pi^2; % Grav. const. * Mass of Sun (au^3/yr^2) mass = 1.; % Mass of comet adaptErr = 1.e-3; % Error parameter used by adaptive Runge-Kutta time = 0; %* Loop over desired number of steps using specified % numerical method. nStep = input(\'Enter number of steps: \'); tau = input(\'Enter time step (yr): \'); NumericalMethod = menu(\'Choose a numerical method:\', ... \'Euler\',\'Euler-Cromer\',\'Runge-Kutta\',\'Adaptive R-K\'); for iStep=1:nStep %* Record position and energy for plotting. rplot(iStep) = norm(r); % Record position for polar plot thplot(iStep) = atan2(r(2),r(1)); tplot(iStep) = time; kinetic(iStep) = .5*mass*norm(v)^2; % Record energies potential(iStep) = - GM*mass/norm(r); %* Calculate new position and velocity using desired method. if( NumericalMethod == 1 ) accel = -GM*r/norm(r)^3; r = r + tau*v; % Euler step v = v + tau*accel; time = time + tau; elseif( NumericalMethod == 2 ) accel = -GM*r/norm(r)^3; v = v + tau*accel; r = r + tau*v; % Euler-Cromer step time = time + tau; elseif( NumericalMethod == 3 ) state = rk4(state,time,tau,\'gravrk\',GM); r = [state(1) state(2)]; % 4th order Runge-Kutta v = [state(3) state(4)]; time = time + tau; else [state time tau] = rka(state,time,tau,adaptErr,\'gravrk\',GM); r = [state(1) state(2)]; % Adaptive Runge-Kutta v = [state(3) state(4)]; end end %* Graph the trajectory of the comet. figure(1); clf; % Clear figure 1 window and bring forward polar(thplot,rplot,\'+\'); % Use polar plot for graphing orbit xlabel(\'Distance (AU)\'); grid; pause(1) % Pause for 1 second before drawing next plot %* Graph the energy of the comet versus time. figure(2); clf; % Clear figure 2 window and bring forward totalE = kinetic + potential; % Total energy plot(tplot,kinetic,\'-.\',tplot,potential,\'--\',tplot,totalE,\'-\') legend(\'Kinetic\',\'Potential\',\'Total\'); xlabel(\'Time (yr)\'); ylabel(\'Energy (M AU^2/yr^2)\');