Home » Posts tagged 'Universal anomaly'

# Tag Archives: Universal anomaly

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

## Universal Kepler’s equation

In this example we solve the universal Kepler’s equation to determine universal anomaly after dt time with a given initial radius r0, velocity v0 and semimajor axis of a spacecraft.

```clc;
clear all;
mu       = 398600;        % Earth’s gravitational parameter [km^3/s^2]
% Initial conditions
r0       = 12000;  %[km] Initial radius
vr0      = 3;      %[km/s] Initial radial speed
dt       = 7200;   %[s] Time interval
a        = -20000; %[km] Semimajor axis
alfa = 1/a;        %[km^-1] Reciprocal of the semimajor axis;

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',Xi)

% Equations from the book
% Orbital Mechanics for Engineering Students,2nd Edition,Aerospace Engineering```
`Universal anomaly X = 173.324 [km^0.5]`

## Stumpff functions

The Stumpff functions C(z),S(z) developed by Karl Stumpff, are used for analyzing orbits using the universal variable formulation.

``` clc;
clear all;
z = -25:0.1:400;
n = size(z);
for i = 1:n(2)
if     (z(i) > 0)
S(i) = (z(i)^0.5 -sin(z(i)^0.5))/z(i)^1.5;
C(i) = (1 - cos(z(i)^0.5))/z(i);
elseif (z(i) < 0)
S(i) = (sinh((-z(i))^0.5) - (-z(i))^0.5)/((-z(i))^1.5);
C(i) = (cosh((-z(i))^0.5) - 1)/(-z(i));
else
S(i) = 1/6;
C(i) = 0.5;
end
end
% Plot
figure(1);
scatter(z(1:500),C(1:500),'.b');
hold on;grid on;
scatter(z(1:500),S(1:500),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(10,0.4,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

figure(2);
scatter(z(501:n(2)),C(501:n(2)),'.b');
hold on;grid on;
scatter(z(501:n(2)),S(501:n(2)),'.r');
xlabel('z');
ylabel('C(z),S(z)');
legend('C(z)','S(z)');
title('Stumpff functions C(z),S(z)');
text(250,0.0010,'smallsats.org','Color',[0 0 0], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
hold off;

```

Equations from the book Orbital Mechanics for Engineering Students, Second Edition, Aerospace Engineering