Code covered by the BSD License  

Highlights from
Cowell's Method - MICE Version

Cowell's Method - MICE Version

by

 

PDF document and MATLAB script that demonstrates using Cowell’s method to predict orbital motion.

cowell_mice.m
% cowell_mice.m    June 18, 2012

% this script demonstrates geocentric orbital motion
% using Cowell's method, JPL Mice routines and DE421

% simulation includes perturbations due to

%  non-spherical earth gravity
%  point-mass sun and moon gravity

% Orbital Mechanics with Matlab

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

clear all;

global dtr atr emu smu mmu

global j2 j4 lgrav mgrav etime0

global rkcoef ccoef scoef

% astrodynamic and utility constants

om_constants;

% angular conversion factor

atr = dtr / 3600;

% load the SPK, PCK and LSK kernels

cspice_furnsh('naif0010.tls');

cspice_furnsh('de421.bsp');

cspice_furnsh('de-403-masses.tpc');

cspice_furnsh('pck00008.tpc');

cspice_furnsh('earth_tod.tf');

% initialize rkf78 integrator

rkcoef = 1;

% define number of differential equations

neq = 6;

% numerical integration tolerance

tetol = 1.0d-12;

% read gravity model coefficients

gmfile = 'egm96.dat';

[ccoef, scoef] = readegm(gmfile);

% extract j2 coefficient

j2 = -ccoef(3, 1);

j4 = -ccoef(5, 1);

% Earth gravity model order and degree

lgrav = 8;

mgrav = 8;

% initial elapsed ephemeris time (seconds)

cdstr0 = '28 Oct 2008 12:20:07.945';

etime0 = cspice_str2et(cdstr0);

% final elapsed ephemeris time (seconds)

cdstrf = '02 Nov 2008 11:50:38.120';

etimef = cspice_str2et(cdstrf);

% initial position vector (EME2000; kilometers)

ri(1)= -460.02497433214972;

ri(2) = 5800.0988304484881;

ri(3) = 3031.577870969686;

% initial velocity vector (EME2000; kilometers/second)

vi(1) = -10.75021388334612;

vi(2) = -1.5266358237317599;

vi(3) = 1.3061896349718098;

% begin simulation

clc; home;

% total simulation time (seconds)

tsim = etimef - etime0;

% load initial position and velocity vectors

yi(1:3) = ri';

yi(4:6) = vi';

clc; home;

fprintf('\n\n  working ...\n');

while(1)
    
    % step size guess

    h = 10.0;

    % initial time
    
    ti = 0.0;
    
    % final time
    
    tf = tsim;

    % integrate from ti to tf

    yfinal = rkf78('ceqm_mice', neq, ti, tf, h, tetol, yi);

    % reset ininitial conditions
    
    yi = yfinal;

    % check for end of simulation

    if (tf >= tsim)
        
        break;
        
    end
    
end

% compute final state vector and orbital elements

rf = yfinal(1:3);

vf = yfinal(4:6);

oev2 = eci2orb1(emu, rf, vf);

ptarg = mice_spkezr('Moon', etimef, 'J2000', 'NONE', 'Earth');

rmoon = ptarg.state(1:3);

vmoon = ptarg.state(4:6);

% selenocentric EME2000 state vector of the spacecraft

rtmp = rf - rmoon';

vtmp = vf - vmoon';

% print results

fprintf('\nprogram cowell_mice\n');

fprintf('\n< geocentric orbital motion - Cowell''s method >\n');

fprintf('\n\ninitial UTC epoch       ');

disp(cdstr0);

fprintf('------------------------------------------------\n');

fprintf('\ngeocentric state vector and orbital elements\n');

fprintf('(Earth mean equator and equinox of J2000)\n');

svprint(ri, vi);

oev = eci2orb1(emu, ri, vi);

oeprint1(emu, oev);

fprintf('\n\nfinal UTC epoch         ');

disp(cdstrf);

fprintf('------------------------------------------------\n');

fprintf('\nselenocentric state vector and orbital elements\n');

fprintf('(Earth mean equator and equinox of J2000)\n');

svprint(rtmp, vtmp);

oev = eci2orb1(mmu, rtmp, vtmp);

oeprint1(mmu, oev);

fprintf('\n\ndegree of gravity model    %12.2f', lgrav);

fprintf('\norder of gravity model     %12.2f \n\n', mgrav);

% unload ephemeris file

cspice_unload('de421.bsp');

Contact us