3 views (last 30 days)

Show older comments

This question was flagged by Walter Roberson

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

Walter Roberson
on 13 Dec 2015

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.

Walter Roberson
on 13 Dec 2015

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.

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

Start Hunting!