satellite.m

Integraton routines in matlab

%
% satellite.m 
%
% Compute the period of a satellite using Gauss-Chebyshev quadrature
%
%

epsilon = 0.5;    % ellipic eccentricity of the satellite trajectory
p2      = 2.09919245059286;
k       = 1.4e-3; % constant describing the flattening of the earth

omega2 = 0.25*(1-epsilon) + k/6;
p1     = k / ( 6 * omega2 * p2); 


n = 10; % number of quadrature points

int = 0;

for i = 1 : n
    x = cos( pi*(2*i-1)/(2*n) );
    r = 1 + 0.5*( 1 + x ) * ( p2 - 1 );
    g = 2 * omega2*( 1 - p1 / r);
    f = 1 / sqrt( 2 * g ) ;

    int = int + f;
end

int = pi/n * int;




 
fprintf(' Approximation of half the period, computed with \n' )
fprintf(' %d  point Gauss Chebyshev quadrature = %14.13f \n\n', n, int )
fprintf(' Value of half the period (quoted in Deuflhard/Hohmann) \n' )
fprintf(' is given by                           4.4395413186376 \n' )