How to add the additional force into the code (possibly state vector)

4 Ansichten (letzte 30 Tage)
Ash matlab
Ash matlab am 27 Jan. 2020
Hi Y'all, I wanted to add a force that pulled the object to the left ( -x axis)
How would I add that into the code. I thought the best way would be in the state vector
%% 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)');

Antworten (0)

Kategorien

Mehr zu Gravitation, Cosmology & Astrophysics finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by