Code covered by the BSD License  

Highlights from
A MATLAB Script for Propagating Interplanetary Trajectories from Earth to Mars

A MATLAB Script for Propagating Interplanetary Trajectories from Earth to Mars

by

 

Numerically integrate the orbital equations of motion of an Earth to Mars interplanetary trajectory.

geo_eqm (t, y)
function ydot = geo_eqm (t, y)

% geocentric equations of motion

% version for tcm_e2m.m

% input

%  t = simulation time (seconds)
%  y = state vector (kilometers and kilometers/second)

% output

%  ydot = integration vector

% Orbital Mechanics with MATLAB

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

global mmu smu xmu aunit req

global imoon isun iplanet

global csrp_geo jdtdb_wrk

agrav = zeros(3, 1);

amoon = zeros(3, 1);

asun = zeros(3, 1);

asrp = zeros(3, 1);

accp = zeros(3, 1);

% current tdb julian date

jdate = jdtdb_wrk + t / 86400.0;

% acceleration due to the earth

agrav = gravity (t, y);

if (isun == 1)
    
    %%%%%%%%%%%%%%%%%%%%%
    % solar perturbations
    %%%%%%%%%%%%%%%%%%%%%
    
    svsun = jplephem(jdate, 11, 3);
    
    rsun = aunit * svsun(1:3);
    
    % compute heliocentric position vector of the spacecraft
    
    rs2sc = y(1:3) - rsun(1:3);
    
    % f(q) formulation
    
    vtmp = y(1:3) - 2.0d0 * rsun(1:3);
    
    dot1 = dot(y(1:3), vtmp);
    
    dot2 = dot(rsun, rsun);
    
    qsun = dot1 / dot2;
    
    fsun = qsun * ((3.0d0 + 3.0d0 * qsun + qsun * qsun) ...
        / (1.0d0 + (1.0d0 + qsun)^1.5d0));
    
    d3sun = norm(rs2sc)^3;
    
    % point-mass gravity of the sun
    
    asun = -smu * (y(1:3) + fsun * rsun(1:3)) / d3sun;
    
end

if (imoon == 1)
    
    %%%%%%%%%%%%%%%%%%%%%
    % lunar perturbations
    %%%%%%%%%%%%%%%%%%%%%
    
    svmoon = jplephem(jdate, 10, 3);
    
    rmoon = aunit * svmoon(1:3);
    
    % compute selenocentric position vector of the spacecraft
    
    rm2sc = y(1:3) - rmoon;
    
    % f(q) formulation
    
    vtmp = y(1:3) - 2.0d0 * rmoon;
    
    dot1 = dot(y(1:3), vtmp);
    
    dot2 = dot(rmoon, rmoon);
    
    qmoon = dot1 / dot2;
    
    fmoon = qmoon * ((3.0d0 + 3.0d0 * qmoon + qmoon * qmoon) ...
        / (1.0d0 + (1.0d0 + qmoon)^1.5d0));
    
    d3moon = norm(rm2sc)^3;
    
    % point-mass gravity of the moon
    
    amoon = -mmu * (y(1:3) + fmoon * rmoon(1:3)) / d3moon;
    
end

if (csrp_geo > 0.0d0)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % solar radiation pressure perturbation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    svsun = jplephem(jdate, 11, 3);
    
    rsun = aunit * svsun(1:3);
    
    % geocentric distance of the sun (kilometers)
    
    rmsun = norm(rsun);
    
    % geocentric unit vector of the sun
    
    usun = svsun(1:3) / rmsun;
    
    % geocentric distance of the spacecraft
    
    rscm = norm(y(1:3));
    
    % compute unit position vector of spacecraft
    
    usat = y(1:3) / rscm;
    
    % determine shadow conditions
    
    a = usat(2) * usun(3) - usat(3) * usun(2);
    b = usat(3) * usun(1) - usat(1) * usun(3);
    c = usat(1) * usun(2) - usat(2) * usun(1);
    
    d = sqrt(a * a + b * b + c * c);
    
    e = dot(usat, usun);
    
    u = asin(0.00460983743 / (rmsun / aunit));
    
    p = asin(0.0046951089 / (rmsun / aunit));
    
    if (e > 0.0)
        
        q = -d;
        
    else
        
        q = d;
        
    end
    
    % increase the Earth's radius by 90 km
    % to account for the atmosphere
    
    ratm = req + 90.0;
    
    b = asin(ratm / rscm);
    
    v = b - u;
    w = b + p;
    
    xx = sin(v);
    yy = sin(w);
    
    % determine shadow conditions
    
    if (q <= yy && q > xx)
        
        % penumbra
        
        csrp = 0.0;
        
    elseif (q <= xx && q >= 0)
        
        % umbra
        
        csrp = 0.0;
        
    else
        
        % sunlight
        
        csrp = csrp_geo;
        
    end
    
    % compute heliocentric position vector and magnitude of the spacecraft
    
    rs2sc = y(1:3) - rsun;
    
    rs2scm = norm(rs2sc);
    
    rs2scm3 = rs2scm * rs2scm * rs2scm;
    
    % compute acceleration vector (km / sec^2)
    
    asrp = 0.001d0 * csrp * aunit^2 * rs2sc / rs2scm3;
    
end

if (iplanet == 1)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % planetary point-mass perturbations
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for i = 1:1:9
        
        svplanet = jplephem(jdate, i, 11);
        
        rplanet = svplanet(1:3);
        
        for j = 1:1:3
            
            rp(i, j) = aunit * rplanet(j);
            
        end
        
    end
    
    % compute planet-centered position vectors of spacecraft (au)
    
    for i = 1:1:9
        
        rp2sc(i, 1) = y(1) - rp(i, 1);
        
        rp2sc(i, 2) = y(2) - rp(i, 2);
        
        rp2sc(i, 3) = y(3) - rp(i, 3);
        
    end
    
    % compute f(q) functions for each planet
    
    for k = 1:1:9
        
        qp(k) = dot(y(1:3), y(1:3) - 2.0 * rp(k, :)') / dot(rp(k, :), rp(k, :));
        
        fp(k) = qp(k) * ((3.0 + 3.0 * qp(k) + qp(k) * qp(k)) / (1.0 + (1.0 + qp(k))^1.5));
        
        d3p(k) = norm(rp2sc(k, :)) * norm(rp2sc(k, :)) * norm(rp2sc(k, :));
        
    end
    
    % compute planetary perturbations
    
    for j = 1:1:3
        
        accp(j) = 0.0;
        
        for k = 1:1:9
            
            xmuwrk = aunit^3 * xmu(k + 1) / 86400.0d0 ^2;
            
            % note: exclude Earth perturbation
            
            if (k ~= 3)
                
                accp(j) = accp(j) - xmuwrk * (y(j) + fp(k) * rp(k, j)) / d3p(k);
                
            end
            
        end
        
    end
    
end

% compute total integration vector

ydot = [ y(4)
    y(5)
    y(6)
    agrav(1) + amoon(1) + asun(1) + accp(1) + asrp(1)
    agrav(2) + amoon(2) + asun(2) + accp(2) + asrp(2)
    agrav(3) + amoon(3) + asun(3) + accp(3) + asrp(3)];


Contact us