Contents
Universal anomaly, Lagrange f and g coefficients
Given the initial position R0 and velocity V0 vector of a satelite at time t0 in earth centered earth fixed reference frame. In this example we show how to calculate position and velocity vectors of the satellite dt minutes later using Lagrange f and g coefficients in terms of the universal anomaly
clc; clear all; R0 = [7200, -13200 0]; %[km] Initial position V0 = [3.500, 2.500 1.2]; %[km/s] Initial velocity dt = 36000; %[s] Time interval mu = 398600; %[km^3/s^2] Earth’s gravitational parameter r0 = norm(R0); v0 = norm(V0); vr0 = R0*V0'/r0; alfa = 2/r0 - v0^2/mu; if (alfa > 0 ) fprintf('The trajectory is an ellipse\n'); elseif (alfa < 0) fprintf('The trajectory is an hyperbola\n'); elseif (alfa == 0) fprintf('The trajectory is an parabola\n'); end X0 = mu^0.5*abs(alfa)*dt; %[km^0.5]Initial estimate of X0 Xi = X0; tol = 1E-10; % Tolerance while(1) zi = alfa*Xi^2; [ Cz,Sz] = Stumpff( zi ); fX = r0*vr0/(mu)^0.5*Xi^2*Cz + (1 - alfa*r0)*Xi^3*Sz + r0*Xi -(mu)^0.5*dt; fdX = r0*vr0/(mu)^0.5*Xi*(1 - alfa*Xi^2*Sz) + (1 - alfa*r0)*Xi^2*Cz + r0; eps = fX/fdX; Xi = Xi - eps; if(abs(eps) < tol ) break end end fprintf('Universal anomaly X = %4.3f [km^0.5] \n\n',Xi) % Lagrange f and g coefficients in terms of the universal anomaly f = 1 - Xi^2/r0*Cz; g = dt - 1/mu^0.5*Xi^3*Sz; R = f*R0 + g*V0; r = norm(R); df = mu^0.5/(r*r0)*(alfa*Xi^3*Sz - Xi); dg = 1 - Xi^2/r*Cz; V = df*R0 + dg*V0; fprintf('Position after %4.2f [min] R = %4.2f*i + %4.2f*j + %4.2f*k[km] \n',dt/60,R); fprintf('Velocity after %4.2f [min] V = %4.3f*i + %4.3f*j + %4.2f*k[km/s] \n',dt/60,V);
The trajectory is an ellipse Universal anomaly X = 1922.210 [km^0.5] Position after 600.00 [min] R = -6781.27*i + -11870.72*j + -3270.69*k[km] Velocity after 600.00 [min] V = 3.488*i + -3.362*j + 0.41*k[km/s]
Ploting trajectory
t = 0; step = 100; %[s] Time step ind = 1; while (t < dt) [R V] = UniversalLagrange(R0, V0, t); Ri(ind,:) = R; [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R); ind = ind + 1; t = t + step; end % Earth 3D Plot Earth3DPlot(1); scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r'); hold off; % Earth topographic map EarthTopographicMap(2,820,420); scatter(Long ,Lat,'.r'); text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',... 'HorizontalAlignment','left','FontSize',14 ); title('Satellite ground track'); hold off; %
Hyperbolic trajectory
clear Ri Long Lat R0 = [7200, -6200 0]; %[km] Initial position V0 = [5.500, 7.500 1.2]; %[km/s] Initial velocity dt = 12000; %[s] Time interval % Ploting trajectory t = 0; step = 100; %[s] Time step ind = 1; while (t < dt) [R V] = UniversalLagrange(R0, V0, t); Ri(ind,:) = R; [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R); ind = ind + 1; t = t + step; end % Earth 3D Plot Earth3DPlot(3); scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r'); % Earth topographic map EarthTopographicMap(4,820,420); scatter(Long ,Lat,'.r'); text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',... 'HorizontalAlignment','left','FontSize',14 ); title('Satellite ground track');
Equations from the book Orbital Mechanics for Engineering Students, Second Edition Aerospace_Engineering
where these functions
[R V] = UniversalLagrange(R0, V0, t);
[Long(ind,:), Lat(ind,:)] = R2RA_Dec(R);