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

deltav_guess(oev1, alttar, fpatar)
function [dv, dtof] = deltav_guess(oev1, alttar, fpatar)

% initial guess for de-orbit delta-v vector and flight time

% input

%  oev1   = classical orbital elements of initial orbit
%  alttar = "target" altitude at EI (kilometers)
%  fpatar = "target" flight path angle at EI (radians)

% output

%  dv   = de-orbit delta-v vector (km/sec)
%  dtof = flight time from maneuver to EI (seconds)

% Orbital Mechanics with MATLAB

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

global mu req

cfpae = cos(fpatar);

if (oev1(2) < 1.0d-10)

    % initial circular orbit

    [r1, v1] = orb2eci(mu, oev1);

    rmag1 = norm(r1);

    rentry = req + alttar;

    % determine deorbit delta-v

    sigma = rmag1 / rentry;

    % local circular velocity at entry altitude

    vlce = sqrt(mu / rentry);

    % required de-orbit delta-v (kilometers/second)

    deltav = vlce * sqrt(1.0 / sigma) ...
        * (1.0 - sqrt(2.0 * (sigma - 1.0) / ((sigma / cfpae)^2 - 1.0)));

    vmag1 = norm(v1);

    vunit1 = v1 / vmag1;

    % compute velocity vector after delta-v

    vmag2 = vmag1 - deltav;

    v2 = vmag2 * vunit1;

    % orbital elements after maneuver

    r2 = r1;

    oev2 = eci2orb1(mu, r2, v2);

    sma2 = oev2(1);

    ecc2 = oev2(2);

    % true anomaly on deorbit trajectory at entry interface (radians)

    tmp1 = mu * (2.0 * sma2 * rentry - rentry^2 - sma2^2 * (1.0 - ecc2^2));

    tmp2 = sma2 * rentry^2;

    rdot = -sqrt(tmp1 / tmp2);

    sta = (rdot / ecc2) * sqrt(sma2 * (1.0 - ecc2^2) / mu);

    cta = sma2 * (1.0 - ecc2^2) / (ecc2 * rentry) - (1.0 / ecc2);

    ta2 = atan3(sta, cta);

    % flight time from impulse to entry (seconds)

    dtof = tof1(mu, sma2, ecc2, pi, ta2);

else

    % initial elliptical orbit

    [r1, v1] = orb2eci(mu, oev1);

    rp = oev1(1) * (1.0 - oev1(2));

    ra = oev1(1) * (1.0 + oev1(2));

    rentry = req + alttar;

    % local circular velocity at entry altitude (kilometers/second)

    vlce = sqrt(mu / rentry);

    % required deorbit delta-v (kilometers/second)

    alpha1 = ra / rentry;

    beta1 = rp / rentry;

    tmp1 = sqrt(2.0 * beta1 / (alpha1 * (alpha1 + beta1)));

    tmp2 = sqrt(2.0 * (alpha1 - 1.0) / (alpha1 * (alpha1^2 - cfpae^2)));

    deltav = vlce * (tmp1 - tmp2 * cfpae);

    % compute unit velocity vector prior to impulse

    vmag1 = norm(v1);

    vunit1 = v1 / vmag1;

    % compute velocity vector after delta-v

    vmag2 = vmag1 - deltav;

    v2 = vmag2 * vunit1;

    % orbital elements after maneuver

    r2 = r1;

    oev2 = eci2orb1(mu, r2, v2);

    sma2 = oev2(1);

    ecc2 = oev2(2);

    % true anomaly on deorbit trajectory at entry interface (radians)

    tmp1 = mu * (2.0 * sma2 * rentry - rentry^2 - sma2^2 * (1.0 - ecc2^2));

    tmp2 = sma2 * rentry^2;

    rdot = -sqrt(tmp1 / tmp2);

    sta = (rdot / ecc2) * sqrt(sma2 * (1 - ecc2^2) / mu);

    cta = sma2 * (1.0 - ecc2^2) / (ecc2 * rentry) - (1.0 / ecc2);

    ta2 = atan3(sta, cta);

    % flight time from maneuver to entry interface (seconds)

    dtof = tof1(mu, sma2, ecc2, pi, ta2);
    
end

% compute unit velocity vector prior to maneuver

vmag1 = norm(v1);

vunit1 = v1 / vmag1;

% de-orbit delta-v vector (kilometers/second)

dv = -deltav * vunit1;

Contact us