I need help with a MATLAB code. I am trying to figure out how to add perturbation to the following code. I want to add J2 Perturbation, Drag, SRP, and third body effects and use that to replace the forbit function. Does MATLAB have an in built function that can help me here. I heard about the numericalPropagator function of MATLAB. Can that be used here? % Keplerian Elements a = 29599.8; e = 0.0001; i = 0.9774; Omega = 1.3549; w = 0; M = 0.2645; [RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M); X = [RECI;VECI]*1e3; %========================================================================== % Propagation %========================================================================== h = 300; % Step Size (sec) steps = 300; % Number of Steps [X_RK] = RK_4(X,h,steps); % Orbit in ECI (PV) %% Funcitons function [X_RK] = RK_4(X,h,steps) X_init = X; for t=1:steps X=rk4orbit(t,X,h); XX(t)=X(1); YY(t)=X(2); ZZ(t)=X(3); VX(t)=X(4); VY(t)=X(5); VZ(t)=X(6); plot3(XX,YY,ZZ); % Plot numerical orbit pause(0.0000000005) end title('Simulated Orbit in ECI') % Create PV Matrix X_RK = vertcat(XX,YY,ZZ,VX,VY,VZ); X_RK = horzcat(X_init,X_RK); end function X = rk4orbit(t,X,h) % Numerical Solution, Runge-Kutta 4th Order k1 = orbitalDynamics(t,X); k2 = orbitalDynamics(t+h/2,X+k1*h/2); k3 = orbitalDynamics(t+h/2,X+k2*h/2); k4 = orbitalDynamics(t+h,X+k3*h); % Step forward in time X = X+(h/6)*(k1+2*k2+2*k3+k4); end function Xdot = forbit(t,X) mu = 3.986004418*10^14; % Gravitational constant (m^3/s^2) r = X(1:3); % Position (m) v = X(4:6); % Velocity (ms^2) dr = v; dv = (-mu/(norm(r))^3).*r; % Newton's law of gravity Xdot = [dr; dv]; end function [RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M) mu = 398600.441799999971; % Earth's Standard Gravitational Parameter (GM) Eps = 2.22044604925031e-8; % Machine epsilon E = SolveKepler(M,e,Eps); % Eccentric Anomaly from Mean Anomaly M = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); % True Anomaly obtained from Eccentric Anomaly p = a*(1 - e^2); % Semiparameter (km) % Position Coordinates in Perifocal Coordinate System x = (p*cos(M)) / (1 + e*cos(M)); % x-coordinate (km) y = (p*sin(M)) / (1 + e*cos(M)); % y-coordinate (km) z = 0; % z-coordinate (km) vx = -(mu/p)^(1/2) * sin(M); % velocity in x (km/s) vy = (mu/p)^(1/2) * (e + cos(M)); % velocity in y (km/s) vz = 0; % velocity in z (km/s) % Transformation Matrix (3 Rotations) t_rot = [cos(Omega)*cos(w)-sin(Omega)*sin(w)*cos(i) ... (-1)*cos(Omega)*sin(w)-sin(Omega)*cos(w)*cos(i) ... sin(Omega)*sin(i); ... sin(Omega)*cos(w)+cos(Omega)*sin(w)*cos(i) ... (-1)*sin(Omega)*sin(w)+cos(Omega)*cos(w)*cos(i) ... (-1)*cos(Omega)*sin(i); ... sin(w)*sin(i) cos(w)*sin(i) cos(i)]; % Transformation Perifocal xyz to ECI RECI = t_rot*[x y z]'; VECI = t_rot*[vx vy vz]'; end function E = SolveKepler(M,e,Eps) % Iterative Solution of Kepler's equation (M = E-e*sin(E)) %========================================================================== % Inputs: M = Mean Anomaly (rad) % e = Eccentricity % Eps = Machine epsilon % Output: E = Eccentric Anomaly (rad) %========================================================================== E0 = M; E1 = E0 - (E0 - e*sin(E0) - M)/(1 -e*cos(E0)); while (abs(E1-E0) > Eps) E0 = E1; E1 = E0 - (E0 - e*sin(E0) - M)/(1 - e*cos(E0)); end E = E1; end

I need help with a MATLAB code. I am trying to figure out how to add perturbation to the following code. I want to add J2 Perturbation, Drag, SRP, and third body effects and use that to replace the forbit function. Does MATLAB have an in built function that can help me here. I heard about the numericalPropagator function of MATLAB. Can that be used here? % Keplerian Elements a = 29599.8; e = 0.0001; i = 0.9774; Omega = 1.3549; w = 0; M = 0.2645; [RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M); X = [RECI;VECI]1e3; %========================================================================== % Propagation %========================================================================== h = 300; % Step Size (sec) steps = 300; % Number of Steps [X_RK] = RK_4(X,h,steps); % Orbit in ECI (PV) %% Funcitons function [X_RK] = RK_4(X,h,steps) X_init = X; for t=1:steps X=rk4orbit(t,X,h); XX(t)=X(1); YY(t)=X(2); ZZ(t)=X(3); VX(t)=X(4); VY(t)=X(5); VZ(t)=X(6); plot3(XX,YY,ZZ); % Plot numerical orbit pause(0.0000000005) end title('Simulated Orbit in ECI') % Create PV Matrix X_RK = vertcat(XX,YY,ZZ,VX,VY,VZ); X_RK = horzcat(X_init,X_RK); end function X = rk4orbit(t,X,h) % Numerical Solution, Runge-Kutta 4th Order k1 = orbitalDynamics(t,X); k2 = orbitalDynamics(t+h/2,X+k1h/2); k3 = orbitalDynamics(t+h/2,X+k2h/2); k4 = orbitalDynamics(t+h,X+k3h); % Step forward in time X = X+(h/6)(k1+2k2+2k3+k4); end function Xdot = forbit(t,X) mu = 3.98600441810^14; % Gravitational constant (m^3/s^2) r = X(1:3); % Position (m) v = X(4:6); % Velocity (ms^2) dr = v; dv = (-mu/(norm(r))^3).r; % Newton's law of gravity Xdot = [dr; dv]; end function [RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M) mu = 398600.441799999971; % Earth's Standard Gravitational Parameter (GM) Eps = 2.22044604925031e-8; % Machine epsilon E = SolveKepler(M,e,Eps); % Eccentric Anomaly from Mean Anomaly M = 2atan(sqrt((1+e)/(1-e))tan(E/2)); % True Anomaly obtained from Eccentric Anomaly p = a(1 - e^2); % Semiparameter (km) % Position Coordinates in Perifocal Coordinate System x = (pcos(M)) / (1 + ecos(M)); % x-coordinate (km) y = (psin(M)) / (1 + ecos(M)); % y-coordinate (km) z = 0; % z-coordinate (km) vx = -(mu/p)^(1/2) * sin(M); % velocity in x (km/s) vy = (mu/p)^(1/2) * (e + cos(M)); % velocity in y (km/s) vz = 0; % velocity in z (km/s) % Transformation Matrix (3 Rotations) t_rot = [cos(Omega)cos(w)-sin(Omega)sin(w)cos(i) ... (-1)cos(Omega)sin(w)-sin(Omega)cos(w)cos(i) ... sin(Omega)sin(i); ... sin(Omega)cos(w)+cos(Omega)sin(w)cos(i) ... (-1)sin(Omega)sin(w)+cos(Omega)cos(w)cos(i) ... (-1)cos(Omega)sin(i); ... sin(w)sin(i) cos(w)sin(i) cos(i)]; % Transformation Perifocal xyz to ECI RECI = t_rot[x y z]'; VECI = t_rot[vx vy vz]'; end function E = SolveKepler(M,e,Eps) % Iterative Solution of Kepler's equation (M = E-esin(E)) %========================================================================== % Inputs: M = Mean Anomaly (rad) % e = Eccentricity % Eps = Machine epsilon % Output: E = Eccentric Anomaly (rad) %========================================================================== E0 = M; E1 = E0 - (E0 - esin(E0) - M)/(1 -ecos(E0)); while (abs(E1-E0) > Eps) E0 = E1; E1 = E0 - (E0 - esin(E0) - M)/(1 - ecos(E0)); end E = E1; end

Elements Of Electromagnetics

7th Edition

ISBN:9780190698614

Author:Sadiku, Matthew N. O.

Publisher:Sadiku, Matthew N. O.

ChapterMA: Math Assessment

Section: Chapter Questions

Problem 1.1MA

See similar textbooks

Related questions

Question

% Keplerian Elements

a = 29599.8;

e = 0.0001;

i = 0.9774;

Omega = 1.3549;

w = 0;

M = 0.2645;

[RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M);

X = [RECI;VECI]*1e3;

%==========================================================================

% Propagation

%==========================================================================

h = 300; % Step Size (sec)

steps = 300; % Number of Steps

[X_RK] = RK_4(X,h,steps); % Orbit in ECI (PV)

%% Funcitons

function [X_RK] = RK_4(X,h,steps)

X_init = X;

for t=1:steps

X=rk4orbit(t,X,h);

XX(t)=X(1);

YY(t)=X(2);

ZZ(t)=X(3);

VX(t)=X(4);

VY(t)=X(5);

VZ(t)=X(6);

plot3(XX,YY,ZZ); % Plot numerical orbit

pause(0.0000000005)

end

title('Simulated Orbit in ECI')

% Create PV Matrix

X_RK = vertcat(XX,YY,ZZ,VX,VY,VZ);

X_RK = horzcat(X_init,X_RK);

end

function X = rk4orbit(t,X,h)

% Numerical Solution, Runge-Kutta 4th Order

k1 = orbitalDynamics(t,X);

k2 = orbitalDynamics(t+h/2,X+k1*h/2);

k3 = orbitalDynamics(t+h/2,X+k2*h/2);

k4 = orbitalDynamics(t+h,X+k3*h);

% Step forward in time

X = X+(h/6)*(k1+2*k2+2*k3+k4);

end

function Xdot = forbit(t,X)

mu = 3.986004418*10^14; % Gravitational constant (m^3/s^2)

r = X(1:3); % Position (m)

v = X(4:6); % Velocity (ms^2)

dr = v;

dv = (-mu/(norm(r))^3).*r; % Newton's law of gravity

Xdot = [dr; dv];

end

function [RECI, VECI] = Kepler2RV(a, e, i, Omega, w, M)

mu = 398600.441799999971; % Earth's Standard Gravitational Parameter (GM)

Eps = 2.22044604925031e-8; % Machine epsilon

E = SolveKepler(M,e,Eps); % Eccentric Anomaly from Mean Anomaly

M = 2*atan(sqrt((1+e)/(1-e))*tan(E/2)); % True Anomaly obtained from Eccentric Anomaly

p = a*(1 - e^2); % Semiparameter (km)

% Position Coordinates in Perifocal Coordinate System

x = (p*cos(M)) / (1 + e*cos(M)); % x-coordinate (km)

y = (p*sin(M)) / (1 + e*cos(M)); % y-coordinate (km)

z = 0; % z-coordinate (km)

vx = -(mu/p)^(1/2) * sin(M); % velocity in x (km/s)

vy = (mu/p)^(1/2) * (e + cos(M)); % velocity in y (km/s)

vz = 0; % velocity in z (km/s)

% Transformation Matrix (3 Rotations)

t_rot = [cos(Omega)*cos(w)-sin(Omega)*sin(w)*cos(i) ...

(-1)*cos(Omega)*sin(w)-sin(Omega)*cos(w)*cos(i) ...

sin(Omega)*sin(i); ...

sin(Omega)*cos(w)+cos(Omega)*sin(w)*cos(i) ...

(-1)*sin(Omega)*sin(w)+cos(Omega)*cos(w)*cos(i) ...

(-1)*cos(Omega)*sin(i); ...

sin(w)*sin(i) cos(w)*sin(i) cos(i)];

% Transformation Perifocal xyz to ECI

RECI = t_rot*[x y z]';

VECI = t_rot*[vx vy vz]';

end

function E = SolveKepler(M,e,Eps)

% Iterative Solution of Kepler's equation (M = E-e*sin(E))

%==========================================================================

% Inputs: M = Mean Anomaly (rad)

% e = Eccentricity

% Eps = Machine epsilon

% Output: E = Eccentric Anomaly (rad)

%==========================================================================

E0 = M;

E1 = E0 - (E0 - e*sin(E0) - M)/(1 -e*cos(E0));

while (abs(E1-E0) > Eps)

E0 = E1;

E1 = E0 - (E0 - e*sin(E0) - M)/(1 - e*cos(E0));

end

E = E1;

end

Expert Solution

Step by step

Solved in 2 steps

SEE SOLUTION Check out a sample Q&A here