How to provide provide extra acceleration at certain time?

I am trying to get my satellite to Jupiter and I would like to provide an extra acceleration (boost) when a certain time is reached. I know I would need to use an if statement, but how would I go about providing this thrust for a certain time and then reverting back to my original equations of motion? I would like to play with certain values for the time that the thrust is applied to see when I could get my satellite closest to Jupiter.

If I used something like if t>some value, would I need to create another for-loop so that the values for the other accelerations/velocities/positions are still be calculated? So would I need the entire numerical calculation repeated for however long I plan on firing the engine? Then how would I go back to calculating the acceleration without the boost involved?

My idea:

Put this right after my a3 is calculated:

 if tc >= some value && tc <= some value
     a3(:,i) = a3(:,i) + boost;
 else

Now, would this apply the boost once, and then continue calculating the other values in the for-loop? Or would it get stuck inside of the if statement?

My code:

 clear all;
close all;
clc

 Ms = 1.98892E30; % [kg] Mass of the sun
mp1 = 6.42E24; % [kg] Mass of earth
mp2 = 2E28; % [kg] Mass of super jupiter
msat = 10000; %[kg] Mass of probe
dt = 2500; % [s] time step
n = 10080; % [=30 years] number of time steps
ti = 0; % [s] Initial time
tf = ti + (n-1) * dt; % [s] Final time
xu=2.992E8; % [km] upper x axis limit
xl=-xu; % [km] lower x axis limit
yu=2.992E8; % [km] upper y axis limit
yl=-yu; % [km] lower y axis limit
limits = [xl xu yu yl]; %predefines limits of axes of graph in array to be called later

 %initial conditions, positions in [km], velocities in [km/s]
pip1 = [1.5E8; 0]; %initial position vector for planet 1 (earth like)
vip1 = [0; 29.75]; %initial velocity vector for planet 1 (earth like)
pip2 = [0; 7.788E8]; %initial position vector for planet 2 (jupiter)
vip2 = [-13.07; 0]; %initial velocity vector for planet 2 (jupiter)
psun = [0;0]; %initial position of sun
psat1 = [1.5E8 + 384400; 0]; %initial orbit of probe
vsat1 = [0; 29.75+ 1]; %initial velocity of probe
boost = [-4;30] %boost vector 

 %predefine arrays
t = linspace(ti,tf,n); %time array
pp1 = zeros(2,n); %position of earth like planet array
vp1 = zeros (2,n); %velocity of earth like planet array
a1 = zeros (2,n); %acceleration of earth like planet array
pp2 = zeros(2,n); %position of  jupiter planet array
vp2 = zeros (2,n); %velocity of jupiter planet array
a2 = zeros (2,n);%acceleration of jupiter planet array
pp3 = zeros(2,n); %position of probe array
vp3 = zeros (2,n); %velocity of probe array
a3 = zeros (2,n); %acceleration of probe array

 % set intial conditions equal to first value in arrays
pp1(:,1) = pip1; %sets first row of array equal to initial position vector of planet 1
vp1 (:,1) = vip1; %sets first row of array equal to initial velocity vector of planet 1
pp2 (:,1) = pip2; %sets first row of array equal to initial position vector of planet 2
vp2 (:,1) = vip2; %sets first row of array equal to initial velocity vector of planet 2
pp3 (:,1) = psat1; %sets first row of array equal to inital positon vectory of probe
vp3 (:,1) = vsat1; %sets first row of array equal to initial velocity vector of probe

 for i = 1:n

   a1(:,i) = grava(pp1(:,i),psun(:,1),Ms) + grava(pp1(:,i),pp2(:,i),mp2); %calculates acceleration on earth
   a2(:,i) = grava(pp2(:,i),psun(:,1),Ms) + grava(pp2(:,i),pp1(:,i),mp1); %calculates acceleration on jupiter
   a3(:,i) = grava(pp3(:,i),pp1(:,i),mp1) + grava(pp3(:,i),pp2(:,i),mp2) + grava(pp3(:,i),psun(:,1),Ms); %caluclates acceleration of probe due to sun/earth/jupiter forces 

   vp1(:,i+1) = vp1(:,i) + a1(:,i) * dt; %calculates velocity of earth
   vp2(:,i+1) = vp2(:,i) + a2(:,i) * dt; %calculates velocity ofjupiter
   vp3(:,i+1) = vp3(:,i) + a3(:,i) * dt;

   pp1(:,i+1) = pp1(:,i) + vp1(:,i+1) * dt; %calculates position of earth using Euler-cromer
   pp2(:,i+1) = pp2(:,i) + vp2(:,i+1) * dt; %calculates position of jupiter using euler-cromer
   pp3(:,i+1) = pp3(:,i) + vp3(:,i+1) * dt;

   tc(i+1) = t(i) + dt;
end

 figure %creates figure
axis([xl xu yl yu]) %sets limits on axis (2 au)
plot(pp1(1,:),pp1(2,:),'r') %plots path of earth (x,y)
hold on; %holds onto above plot
plot(pp2(1,:),pp2(2,:),'b') %plots path of jupiter (x,y)
hold on; %holds onto above plot
plot(0,0,'*') %plots a single point to stand for the sun
hold on; %holds onto above plot
plot (pp3(1,:),pp3(2,:),'o') %probe plot
hold on;
xlabel('X (km)') %labels x axis
ylabel('Y (km)') %lables y axis
legend ('Earth','Jupiter','Sun','Probe'); %adds a legend
title('Orbits of Earth/Jupiter/Probe about Sun') %gives plot title

My function:

 function a = grava(p,p0,M)

 G = 6.673E-11/(1000^3); %divided by 1000^3 to make units consistant for using km
R = p - p0;
a = -G*M * R / (sum(R.*R)^(3/2));

 return

[Rik]

Code restored from Google cache.

Answer 1

[Walter Roberson]

Write your code to always have a thrust at each time step, so that you only have one version of the equations of motion. But the thrust value can be a function of the time, or it could be an array indexed at the time step number, and you would use [0,0,0] for the portion where there was no thrust. A vector of values because the thrust is directional in three-space. (If your thrust is to be taken relative to a fixed position on the spacecraft and that craft is spinning then you will need additional parameters to change the frame of reference)

If you are planning to code in atmospheric drag then you might as well include that as well, by way of functions of position and velocity that can return [0,0,0] for now.

[Ibro Tutic]

Thanks for the reply. I think that I am going to use an array indexed at a time value, however, would I need to use an if statement to begin inputting thrust values into the array? So something like:

 If tc >= some value && tc <= some value
    Boost(i)= thrust
 Else
    Boost(i) = 0
 End

Would this code stop within the if statement or would it see that time is Inbetween those values, set boost = thrust, and then proceed on with the numerical integration ?

[Walter Roberson]

time_array = linspace(0, last_time, number_of_points);

Boost_magnitude = zeros(size(time_array));
Boost_magnitude(time_array >= boost_start_time & time_array <= boost_end_time) = thrust;

and decide how you are going to move between magnitude vector and x/y/z components.