Code covered by the BSD License  

Highlights from
A MATLAB Script for Optimal Single Impulse De-orbit from Earth Orbits

from A MATLAB Script for Optimal Single Impulse De-orbit from Earth Orbits by David Eagle
optimal impulsive maneuver required to de-orbit a spacecraft in a circular or elliptical Earth orbit

deorbit_shoot(x)
function [f, g] = deorbit_shoot(x)

% objective function and EI equality constraints

% input

%  x(1) = current value of maneuver true anomaly (radians)
%  x(2) = x-component of deorbit delta-v vector (kilometers/second)
%  x(3) = y-component of deorbit delta-v vector (kilometers/second)
%  x(4) = z-component of deorbit delta-v vector (kilometers/second)
%  x(5) = flight time from maneuver to EI (seconds)

% output

%  f(1) = objective function (deorbit delta-v magnitude)
%  f(2) = EI geodetic altitude constraint error (kilometers)
%  f(3) = EI relative flight path angle constraint error (radians)

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global mu omega jdate0 gst0 oevpo

global alttar fpatar rei vei jdatei

% compute current state vector prior to deorbit delta-v

oewrk = oevpo;

oewrk(6) = x(1);

[rpo, vpo] = orb2eci(mu, oewrk);

% current state vector after deorbit delta-v

ri = rpo;

vi = vpo + x(2:4);

% current guess for flight time (seconds)

tof = x(5);

% propagate from maneuver to current guess for time at EI

options = odeset('RelTol', 1.0e-8, 'AbsTol', 1.0e-10);

[twrk, ysol] = ode45(@j2eqm, [0 tof], [ri vi], options);

% extract state vector of the spacecraft at entry interface

rei = ysol(length(twrk), 1:3);

vei = ysol(length(twrk), 4:6);

% scalar objective function (deorbit delta-v magnitude)

f(1) = norm(x(2:4));

% julian date at entry interface

jdatei = jdate0 + tof / 86400.0;

% greenwich apparent sidereal time at entry interface (radians)

gast = mod(gst0 + omega * tof, 2.0 * pi);

% compute relative flight path coordinates at entry interface

fpc = eci2fpc1(gast, rei, vei);

% compute geodetic coordinates at entry interface

rmag = fpc(5);

dec = fpc(2);

[alt, lat] = geodet1 (rmag, dec);

% geodetic altitude error at entry interface (kilometers)
    
f(2) = alt - alttar;

% relative flight path angle error at entry interface (radians)

f(3) = fpc(3) - fpatar;

% transpose objective function

f = f';

% no derivatives

g = [];

Contact us