Hi there,

I have come up with both scripts but something is causing it not to work. I have been trying for a while but to no avail. Was wandering if someone could tell me what is wrong with my code. Thanks!

My scripts are as follows: satellite.m:

function [T,X,Y,Z,U,V,W] = satellite(Xo,Yo,Zo,Uo,Vo,Wo,tstart,tend,maxthrust)

%% SATELLITE has inputs which are components of the initial position

%%(Xo, Yo, Zo) and initial velocity (Uo, Vo, Wo). The engine is turned

%%on at time (tstart) and off at time (tend) with a maximum value

%%of thrust (maxthrust).

G = 6.67408*10^-11;

Me = 5.97*10^24;

Re = 6.37*10^6;

m = 1500;

dt =1;

U=Uo;V=Vo;W=Wo;

X=Xo;Y=Yo;Z=Zo;

T = 1;

total = 0;

for n = 1:10^10

U(n+1) = U(n)*(Xthrust/m) -G*Me*(X(n)/((X(n).^2+Y(n).^2+Z(n).^2).^3/2))*dt;

V(n+1) = V(n)*(Ythrust/m) -G*Me*(Y(n)/((X(n).^2+Y(n).^2+Z(n).^2).^3/2))*dt;

W(n+1) = W(n)*(Zthrust/m) -G*Me*(Z(n)/((X(n).^2+Y(n).^2+Z(n).^2).^3/2))*dt;

X(n+1) = X(n) + U(n+1)*dt;

Y(n+1) = Y(n) + V(n+1)*dt;

Z(n+1) = Z(n) + W(n+1)*dt;

T(n+1) = T(n) + 1;

x=X(n+1)-X(n);

yt=Y(n+1)-Y(n);

z=Z(n+1)-Z(n);

h = sqrt.(X^2 + Y^2 + Z^2) - Re;

Vmag = sqrt(U^2 + V^2 + W^2);

Acc = Vmag/dt;

n = n+1;

total = total + sqrt(x.^2 + y.^2 + z.^2);

if total >= 4.2*10^8

break

end

end

end

% end function satellite

read_input.m

function [Xo,Yo,Zo,Uo,Vo,Wo,tstart,tend,maxthrust] = read_input(inputfile,sat_id)

%READ_INPUT takes the inputs inputfile (string denoting name of the file to

%be read) and sat_id (integer indicating the satellite number).

%It returns outputs [Xo,Yo,Zo,Uo,Vo,Wo,tstart,tend,maxthrust] where

%(Xo,Yo,Zo) is the initial position, (Uo,Vo,Wo) is the initial velocity,

%(tstart,tend) are the start and end times of the engine and (maxthrust)

%is the maximum thrust produced.

%Call format: [Xo,Yo,Zo,Uo,Vo,Wo,tstart,tend,maxthrust]=read_input(inputfile,sat_id)

read = importdata(inputfile);

for n=1:length(read.data)

if rem(read.data(n,1),1)==0

initial=read.data(sat_id,2:end);

sat_id=inital(1);Xo=initial(2);Yo=initial(3);Zo=initial(4);

Uo=initial(5);Vo=initial(6);Wo=initial(7);

tstart=initial(8);tend=initial(9);maxthrust=initial(10);

else

Xo=NaN;

Yo=NaN;

Zo=Nan;

Uo=NaN;

Vo=NaN;

Wo=Nan;

tstart=NaN;

tend=NaN;

maxthrust=NaN;

disp('Error:Invalid input');

end

end

end

% end function read_input

engine.m :

function [Xthrust, Ythrust, Zthrust] = engine(tstart, tend, maxthrust, t, u, v, w)

% ENGINE produces the thrust needed to change satellite orbits.

% Inputs are engine start and end times (tstart, tend), the maximum thrust

% produced (maxthrust), current time (t), and current velocity components

% (u,v,w).

% Call format: [Xthrust, Ythrust, Zthrust] = engine(tstart, tend, maxthrust, t, u, v, w)

period = tend - tstart;

tpeak = tstart + 0.5*period;

if t > tstart && t

Xthrust = maxthrust*exp(-((t-tpeak)/(0.25*period))^4) ...

*u/sqrt(u^2+v^2+w^2);

Ythrust = maxthrust*exp(-((t-tpeak)/(0.25*period))^4) ...

*v/sqrt(u^2+v^2+w^2);

Zthrust = maxthrust*exp(-((t-tpeak)/(0.25*period))^4) ...

*w/sqrt(u^2+v^2+w^2);

else

Xthrust = 0;

Ythrust = 0;

Zthrust = 0;

end

end % function engine

Contents of satellite_data are as follows.

# Satellite parameters in the following order from left to right:

# Satellite ID, components of initial position (X, Y, Z) (m), components of initial velocity (U, V, W) (m/s)

# Engine start time (s), end time (s), maximum thrust (N)

1 0.00000000e+00 -5.9495945e+06 -3.4350000e+06 7.61560659e+03 0.00000000e+00 0.00000000e+00 20500 24100 350

2 -5.9495945e+06 0.00000000e+00 -3.4350000e+06 0.00000000e+00 7.61560659e+03 0.00000000e+00 20500 24100 300

3 0.00000000e+00 3.43500000e+06 -5.9495945e+06 7.61560659e+03 0.00000000e+00 0.00000000e+00 20500 24100 250

4 -3.4350000e+06 -5.9495945e+06 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.61560659e+03 20500 24100 200

5 3.43500000e+06 0.00000000e+00 -5.9495945e+06 0.00000000e+00 7.61560659e+03 0.00000000e+00 20500 24100 150

6 0.00000000e+00 -5.9495945e+06 3.43500000e+06 7.61560659e+03 0.00000000e+00 0.00000000e+00 20500 24100 100