%% satellite.m %% Compute the period of a satellite using Gauss-Chebyshev quadrature%%epsilon = 0.5;    % ellipic eccentricity of the satellite trajectoryp2      = 2.09919245059286;k       = 1.4e-3; % constant describing the flattening of the earthomega2 = 0.25*(1-epsilon) + k/6;p1     = k / ( 6 * omega2 * p2); n = 10; % number of quadrature pointsint = 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;endint = 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' )