Solar Navigation System

 

He Lin

Abstract
1 Introduction
1.1 Lagrange point
1.2 Global Positioning System
1.3 Solar sail propulsion
2 Simulations of a solar sail trajectory to Mars L4 point.
2.1 Simplified SRP force model
2.2 Mission design
2.3 Earth departure
2.4 Heliocentric transfer
Conclusion
References
Appendix Matlab Codes

Abstract

This project is focused on establishing a new satellite navigation system aimed at the positioning task of deep space missions. The proposed solar navigation system works with the same theory as GPS system and will have two navigation beacons at the Sun-Mars L4/L5 Lagrange point to decrease geometry dilution of precision with solar sail as propulsion. This project paper introduced basic theories of Lagrange point, the GPS dilution of precision and solar sail. And a minimal solar sail trajectory to place the satellite beacon at Mars Lagrange point is simulated with Matlab.

Introduction

As we launch more and more probes in to the solar system, interstellar flight mission controls requires more and more precision in orbit determination. Even if future technology developments enable autonomous flight, orbit determination is still one crucial factor in space flights.

Nowadays precise orbit determination largely relies on the Deep Space Network (DSN) which is a network of antennas that tracks interplanetary spacecraft missions. DSN angular accuracy is shown in Fig 1 below. We can see from the figure below that orbit determination accuracy at Mars orbit still has about 2.250 km (year 2000, estimated), which is not accurate enough to support precise maneuvers.



Fig 1. DSN Angular Accuracy vs. Time (courtesy Nerom)

On the other hand, the fast development of Global Positioning System (GPS) has well proven its capability in near-earth orbital determination missions. Its nature of passive positioning brings more benefit as the system can support a very large group of users.

1.1 Lagrange point

The Lagrange points are five stationary positions in a restricted three body system where the third body has negligible mass. A figure below demonstrates the layout of the problem.


Fig 2. Demonstration of restricted three body problem

In 1772, Lagrange discovered a set of equations that describes the dynamics of the third body.



The above equations have five solutions when  are all zero. The Lagrange points mark positions where the combined gravitational pull of the two large masses provides precisely the centripetal force required to rotate with them.

The positions of the Lagrange points L1 ~ L5 in Earth-Sun system are shown in figure below.



Fig 3. Earth Sun system Lagrange points (courtesy Wikipedia)

And also from the potential curves (white lines in Fig 3) we can see that L4/L5 points are in the low potential areas, which means L4 and L5 are stable solutions. Objects on Earth orbit in the vicinity of L4/L5 points will have a tendency to stay there. This minimizes the need for altitude maneuvers to be carried out on navigation beacons in that position. Which means the space vehicles need to carry less fuel.

1.2 Global Positioning System

The Global Positioning System (GPS) is a satellite navigation system which uses a constellation of 24 Medium Earth Orbit satellites that transmit precise microwave signals, enabling GPS users to determine their current location, the time, and velocity.

1.2.1 Global Positioning System Dilution of Position

Dilution of Position (DOP) is an important factor in calculating the position of the users. It denotes the impact of the geometry of positioning satellites on the position solution accuracy of the users.
To study DOP we need to take a look into the linear model for position solution for GPS systems. [2]

The initial guess of the position and clock bias of the user is
The approximation based on the initial guess is
The true position and true clock bias is where  are unknown, we now have a set of linear equations:
The position of





Where



Now we have




So far we have k equations (k>=4), which can be used to solve for the position of the user.



The solution above can be used to update the estimates of the user position and clock bias and henceforth iterate this will yield the position and clock bias of the user.
Now recall that






Denoting the i-th diagonal element of H as, it exists that





And,




Taken into consideration the clock bias, the Geometric dilution of precision (GDOP) is


GDOP =


Position DOP is


PDOP =


Horizontal DOP is


HDOP =

1.2.2 GPS DOP for far-earth orbits

The GPS system was designed for navigation on earth surface rather than interplanetary navigation so although GPS system has very high accuracy on earth, its position solution result accuracy for interplanetary missions are quite poor. Viewed from planets like Mars, all GPS beacons are clustered around the earth, which have almost the same position vector. The following figure shows the possible best DOP factor for space vehicles at different planet orbits. The result shows intolerable bad position solution which prevents the GPS system to be used in far-earth missions.


Fig 4. GPS best DOP of orbit near different planets


Fig 5. GPS optimal HDOP of different orbits in inner solar system

1.2.3 Improving precision with beacon satellites far from Earth

From the last section we can conclude that the layout of current GPS system greatly restricted the usage of GPS in interplanetary mission. Beacons satellites using the same technique of GPS but located at the Sun-Mars Lagrange points will greatly contribute to the precision of the position solution.



Fig 6. Layout of the projected positioning beacons

This project only involves with deploying navigation beacons to the Sun-Mars L4/L5 position, so the vertical DOP will not be improved. The following calculating only involves with the HDOP to demonstrate the improvement of horizontal DOP.

For simplicity, calculation of HDOP is performed in the situation shown in the figure. The Earth is on the opposite position of the Sun to the Mars. Two navigation beacons are deployed at the L4/L5 points of Mars and the existing GPS orbiting Earth are considered one more source for the navigation solution. Calculations are performed on the positions of planets other than Earth and they are on the same line of Earth-Sun-Mars, on the same side of the Sun with Mars.


Fig 7. Positions of improved HDOP calculation

Fig 8 below shows the HDOP result of the projected solar navigation system. As seen from the figure, HDOP are greatly improved thanks to the more spread out configuration of navigation beacons.  Considering the HDOP of solar navigation system at distance of Pluto is even less than the HDOP of the HDOP of existing GPS system at the distance of Moon (Fig 5), solar navigation system can well handle the navigation needs for outer space missions.

Although future improvements in VDOP and GDOP are out of the reach of this project, simply put two beacons at a reasonable distance (several AUs) above/below the Sun with respect to the Sun-Mars orbital plane could solve the problem of poor VDOP of this proposed solar navigation system.



Fig 8. Result of the improved HDOP

1.3 Solar sail propulsion

Solar sail propulsion uses a big thin film to collect the pressure from photons and henceforth generate the thrust. Although solar sail is often mentioned in science fictions in the past, the fundamental theory was by no means new. James Maxwell first demonstrated the pressure of light in 1873[3], and the solar sail concept was proposed by two Imperial Russian rocket scientists Konstantin Tsiolkovsky and Fridrich Sander in the 1920s. Solar sail was not tested until 1993, when a 20 meter wide solar sail was deployed from the Soviet space station Mir[4]. In 2004, Japan successfully deployed two solar sails from a sounding rocket. Although no successful thrust test had been performed up to date.



Fig 9. Solar sails deployed from Japan sounding rocket (courtesy nasa.org)

1.3.1 Solar Sail force model

The main thrust force for solar sail is Sun Radiation Pressure (SRP). At distance r from the sun, the SRP is described with the function[5]:



Where .

Wright developed a higher-order force model for solar sails that uses a set of optical coefficients  to construct the characteristics of the solar sail material. [6]
 reflection coefficient
s: the specular reflection factor
: emission coefficient of the front side
: emission coefficient of the back side
: non-Lambertian coefficient of the front side, the angular distribution of the emitted and the diffusely reflected photons.
 non-Lambertian coefficient of the back side

The SRP force based on the above parameters is:



Fig 10. SRP force of a non-prefect force model



Where A is the sail area and




And





Define , the Wright SRP force model may be written as

2 Simulations of a solar sail trajectory to Mars L4 point.

2.1 Simplified SRP force model

According to Wright, a highly reflective aluminum-coated front side and with a highly emissive chromium-coated back side solar sail has  And calculated  This denotes that the non-perfect terms of the solar radiation force is rather small. For simplicity we will use a force model of in the simulation.

2.2 Mission design

The mission consists of three main parts: after initial launch for earth, the space vehicle deploys its solar sail and start accelerating, gaining to escape from earth gravity field. After the space vehicle gains enough  and reaches the Sphere of Influence (SOI) boundary of Earth, the space vehicle enters a heliocentric orbit and continues to gain  In the last section of the orbit transfer the SV drops into the vicinity of the Mars orbit L4/L5 point and a minimal thrust is used to final push the SV into the L4/L5 Lagrange point. For simplicity, in the simulation Sun, Earth, Mars and the SV are all in a same orbital plane.

2.3 Earth departure

For this phase of the transfer, we use an ECI frame to calculate the accumulation of. There are two forces in the system which are gravitational force on the direction of and solar radiation force on the direction of. With this steering law, the solar sail is directly perpendicular to the sun-SV vector and hence collects the most energy which shortens the time needed to reach Earth’s SOI.



Fig 11. Earth departure

A Matlab script was written to simulate this phase of orbit. In the simulation, satellite mass is 200 kg and solar sail area is 100000. For simplicity, Earth center mass and Sun center mass models are used. Using the Ode45 function in Matlab to integrate accelerations caused by gravity and solar radiation pressure we get the following diagram. From fig 12 we can see that at approximately 116.3 days from deployment of solar sail the space vehicle escapes the Earth’s SOI. Fig 13 shows the space vehicle orbit trajectory for 173 days.


Fig 12. SV Distance from Earth vs epoch time


Fig 13. SV orbit during near-Earth acceleration phase

2.4 Heliocentric transfer

After the SV leaves the SOI of Earth, the coordinate system is changed from ECI to Sun Inertial frame. The position vector and velocity vector are changed (See Fig 11).


For a simple steering law model, the solar sail uses a fixed angle towards incoming sunlight (Fig 14). [7]



Fig 14. Angles in the solar sail
We have the following relationship:

The angleis the only one parameter in the heliocentric transfer stage. In order to minimize the  needed reentering Mars orbit, several  values are tried.  seems to have the best result. The trajectory in Fig 15 is for a. In the figure we can see that the space vehicle slowly accelerates because of the solar pressure and reaches the Mars orbit.

 

 
Fig 15. Heliocentric transfer orbit and details of the transfer orbit intersecting Mars orbit

 


Fig 16. Distance of SV from the Sun vs. epoch time

When the SV has reached the Mars orbit, a  will push the space vehicle into the Mars orbit. varies with . The table below shows several tested.

Table 1. needed to decelerate the SV vs. solar sail angle


 (deg)

 (m/s)

4.

Transfer orbit does not intersect Mars orbit

4.05

712

4.1

11052

4.2

15959

4.5

24504

5

32439

5.1

33591

As can be seen that  has the smallest  requirement.

Comparing to Hoffman Transfer from Earth orbit to Mars orbit, the solar sail has significantly smaller which is 712m/s comparing to 5594 m/s for Hoffman Transfer. And also, the low  didn’t come along with very long transfer time. Transfer time for solar sail is 314.5 days comparing to 258.9 days for Hoffman transfer.

Conclusion

The proposed solar navigation system will greatly enhance the position solution accuracy for space vehicles far from the Earth. And the special orbit of the beacons will reduce the need for altitude control. Simple simulation revealed the advantages of solar sail propulsion by means of both total needed and tolerable travel time. More precise simulations can be carried out by including 3-body effects and perturbations introduced by Sun, Earth and Mars, and improved mission design can include the gravity assist of the Moon, but those are out of the scope of this class project.

References

[1] Martin Lo, Shane Ross, The Lunar LI Gateway: Portal to the Stars and Beyond
[2] Pratap Misra, Per Enge, Global Positioning System, Signals Measurements and Performance, Second Edition
[3] Colin R. McInnes Solar Sailing Technology, Dynamics and Mission Applications
[4] Solar sail from Wikipedia. http://en.wikipedia.org/wiki/Solar_sail
[5] Bernd Dachwald, Bong Wie, Solar Sail Trajectory Optimization for Intercepting, Impacting, and Deflecting Near-Earth Asteroids
[6] Wright, J., Space Sailing, Gordon and Breach Science Publishers, Philadelphia, 1992.
[7] Scott Piggott, Solar Sail Trajectory Design and Control in Unrestricted Frames

Appendix Matlab Codes

 

Code download:

dop.m

simulator.m

near_earth.m

solarforce.m

heliocentric_ss.m

elorb.m

test.m

 

function dd = dop
R =[42164,384400,91688940,41389422,78341163,628700338,1279796110,2725440592,4354851746,5766200977];
% R = [100:100:200000];
r = 26600;
a = ones(4,4);

%% GDOP
for index = 1:length(R)

a(1,2) = r;a(1,3) = 0;
a(2,2) = -r/2;a(2,3) = -r*sqrt(3)/2;
a(3,2) = -r/2;a(3,3) = r*sqrt(3)/2;
a(4,2) = 0;a(4,3) =0;
a(:,1) = R(index);
a = a/sqrt(R(index)^2 + r^2);
a(:,4) = 1;
a(4,1) = 1;
a(4,3) = 0;
%     disp(a);
D = inv(a'*a);
dopg(index) = sqrt(D(1,1) + D(2,2) + D(3,3)+D(4,4));
end
% disp(dopg);

figure(3);
loglog(R(1),dopg(1),'o',R(2),dopg(2),'o',R(3),dopg(3),'o',R(4),dopg(4),'o',R(5),dopg(5),'o',R(7),dopg(7),'o',R(9),dopg(9),'o',R(10),dopg(10),'o');
% grid on;
title('GPS best DOP of orbit near different planets');
xlabel('Distance from earth (km)');
legend('Geosync','Moon','Mercury','Venus','Mars','Saturn','Neptune','Pluto');
ylabel('GPS DOP');

figure(4);
loglog(R(1),dopg(1),'o',R(2),dopg(2),'o',R(4),dopg(4),'o',R(5),dopg(5),'o');
% grid on;
title('GPS best DOP of orbit near different planets');
xlabel('Distance from earth (km)');
legend('Geosync','Moon','Venus','Mars');
ylabel('GPS DOP');

%% HDOP
for index = 1:length(R)

a(1,2) = r;a(1,3) = 0;
a(2,2) = -r/2;a(2,3) = -r*sqrt(3)/2;
a(3,2) = -r/2;a(3,3) = r*sqrt(3)/2;
a(4,2) = 0;a(4,3) =0;
a(:,1) = R(index);
a = a/sqrt(R(index)^2 + r^2);
a(:,4) = 1;
a(4,1) = 1;
a(4,3) = 0;
%     disp(a);
D = inv(a'*a);
doph(index) = sqrt(D(1,1) + D(2,2));
end
% disp(doph);

figure(5);
loglog(R(1),doph(1),'o',R(2),doph(2),'o',R(3),doph(3),'o',R(4),doph(4),'o',R(5),doph(5),'o',R(6),doph(6),'o',R(7),doph(7),'o',R(9),doph(9),'o',R(10),doph(10),'o');
% grid on;
title('GPS best HDOP of different orbits in inner solar system');
xlabel('Distance from earth (km)');
legend('Geosync','Moon','Mercury','Venus','Mars','Jupitor','Saturn','Neptune','Pluto');
ylabel('GPS HDOP');

figure(6);
loglog(R(1),doph(1),'o',R(2),doph(2),'o',R(3),doph(3),'o',R(4),doph(4),'o',R(5),doph(5),'o');
grid on; grid minor;
title('GPS optimal HDOP of different orbits in inner solar system');
xlabel('Distance from earth (km)');
legend('Geosync','Moon','Mercury','Venus','Mars');
ylabel('GPS HDOP');



%% improved DOP

a = zeros(3,2);
r1 = [-149598023,0];
r2 = [227939186/2, 227939186*sqrt(3)/2];
r3 = [227939186/2, -227939186*sqrt(3)/2];


R = [57909083   108208601   149598023   227939186   778298361   1429394133  2875038615  4504449769  5915799000];

for i = 1:length(R)

r=[R(i),0];
a(1,1:2) = (r-r1)/norm(r-r1);
a(2,1:2) = (r-r2)/norm(r-r2);
a(3,1:2) = (r-r3)/norm(r-r3);
D = inv(a'*a);
doph_improved(i) = sqrt(D(1,1) + D(2,2));
end

figure(7)
loglog(R(1),doph_improved(1),'o',R(2),doph_improved(2),'o',R(4),doph_improved(4),'o',R(5),doph_improved(5),

'o',R(6),doph_improved(6),'o',R(7),doph_improved(7),'o',R(8),doph_improved(8),'o',R(9),doph_improved(9),'o');
grid on; grid minor;
title('improved optimal HDOP');
xlabel('Distance from earth (km)');
legend('Mercury','Venus','Mars','Jupitor','Saturn','Uranus','Neptune','Pluto');
ylabel('GPS HDOP');

 

function xv2 = simulator

%% initial parameters setup
global Er
Er = 149598023e3; % earth sun distance
global Mr
Mr = 227939186e3; % mars sun distance
global miue
miue = 3.986004415e14; %earth miu
global mius
mius = 132712440018e9; %sun miu
ifplot = 1; % whether generate plots or not
%% satellite parameters
global sv_m
sv_m = 200; % sv mass
global A
A = 100000; % solar sail area

%% integrate orbit within earth SOI
time = 0:10:10051520;%original guess is 1.5e7, at 10051520 sec the sv is out of earth SOI

r0 = [0,-8000e3]; %initial position vector wrt earth
v0 = [sqrt(miue/norm(r0)),0]; %initial speed vector wrt earth

tol=0.00001;
options=odeset('RelTol',tol,'AbsTol',[tol tol tol tol]);
[t,xv]=ode45('near_earth',time,[r0, v0],options);

%% plot these stuff
if(ifplot)
dist = zeros(1,length(xv)); %calc distance from earth
for ii = 1:length(xv)
dist(ii) = norm(xv(ii,1:2));
end
figure(1);plot(t/86400,dist,'b',[t(1),t(length(t))]./86400,[9.2465e8,9.2465e8],'r')
xlabel('time from start of solar sail deployment (days)');ylabel('distance (m)');title('Distance from Earth vs epoch time');
legend('SV distance from earth center','Earth SOI boundary');
Esoi =9.2465e8;
theta = 0.01:0.01:2*pi;
circle.x = Esoi*sin(theta);
circle.y = Esoi*cos(theta);
figure(2);plot(xv(:,1),xv(:,2),'b',circle.x,circle.y,'r');axis equal;
legend('SV orbit','Earth SOI boundary');
xlabel('x position');ylabel('y position');title('SV orbit during near-Earth acceleration phase');
end

%% heilocentric orbits integration

% position and velocity vector wrt earth at beginning of heliocentric
% transfer
r0 = [-9.224201687029069e8   0.641963434426086e8];%r0 = xv(length(t),1:2);
v0 = [-3.894661347310160e2  -0.540691492389036e2];%v0 = xv(length(t),3:4);
% position and velocity vector wrt sun
r0 = r0 + [0, Er];
v0 = v0 - [sqrt(mius/Er),0];% can also be minus

time = 0:100:24000000;

tol=0.00001;
options=odeset('RelTol',tol,'AbsTol',[tol tol tol tol]);
[t,xv2]=ode45('heliocentric_ss',time,[r0, v0],options);


theta = 0:0.01:2*pi;
circle.x = Mr*sin(theta);circle.y = Mr*cos(theta);
circle2.x = Er*sin(theta);circle2.y = Er*cos(theta);
figure(3);plot(xv2(:,1),xv2(:,2),'b',circle.x,circle.y,'r',circle2.x,circle2.y,'k');axis equal;
xlabel('x');ylabel('y');title('heliocentric acceleration');
legend('SV trajectory','Mars orbit','Earth orbit');

end

%% integragters
function xdot =near_earth(t,x) %#ok<INUSL>
global miue;
global Er;
global sv_m;
global A;
r0 = x(1:2);

ag = [-miue*r0(1)/(norm(r0))^3, -miue*r0(2)/(norm(r0))^3]; % gravational acceleration

r = [r0(1),r0(2)+Er]; % position vector of sv wrt sun
atemp = solarforce(A,r)/sv_m;
if r0(1)>0
as = [atemp/norm(r)*r(1),atemp/norm(r)*r(2)];% SRP accel, if x>0 it is in a leaving orbit. as != 0
elseif r0(1)<=0
as = [0,0];% SRP accel, if x < 0 it is in a return orbit, as = 0;
end
a = ag + as; % total accel is sum of two forces on the sv
xdot = [x(3);x(4);a(1);a(2) ];

end

function f = solarforce(A,r)
% calculates the SRP force on a solar sail
% hl
% 12.09.08
p = SRP(norm(r));
alpha = 0;% for max thrust solar sail is always perpendicular to position vector
% add complex steering laws here
f = 2*p*A*cos(alpha);
end
function p =SRP(r)
% calculates the spolar radiation pressure
% at distance r from the sun
% hl
% 12.09.08
s0 = 1368; %flux at earth orbit
r0 = 149598023e3;% earth sun distance
c = 3e8;
p = 2*s0/c*r0^2/r^2;
end

 

%% integragters
function xdot =heliocentric_ss(t,x) %#ok<INUSL>
% heliocentric integrator with solar sail
% hl
% 12.09.08
global mius;
global sv_m;
global A;
r = x(1:2); % position vector wrt sun

ag = [-mius*r(1)/(norm(r))^3, -mius*r(2)/(norm(r))^3]; % gravational acceleration

atemp = solarforce(A,r)/sv_m;

alpha = atan2(r(2),r(1));
beta = 4.05/180*pi; % 5 deg
theta =pi/2+alpha -beta;
as = [atemp*sin(beta)*sin(theta),atemp*sin(beta)*cos(theta)];
% as = [0,0]; % no accel
% as = [atemp/norm(r)*r(1),atemp/norm(r)*r(2)];% SRP accel ====fastest accel====
a = ag + as; % total accel is sum of two forces on the sv
xdot = [x(3);x(4);a(1);a(2) ];

end

 

function [a, e] = elorb(R0, V0, miu)
% in m, m/s, deg.
% specially modified for project. for normal elorb see assignment folder

global Mr;
R(1) = R0(1);R(2) = R0(2); R(3) = 0;
V(1) = V0(1);V(2) = V0(2); V(3) = 0;
H = cross(R, V);
h = norm(H);

N = cross([0;0;1],H);

v = norm(V);
r = norm(R);

E = ((v^2 - miu/r)*R - dot(R,V)*V)/miu;
e = norm(E);


xi = v^2/2 - miu/r;


if( norm(E) ~= 1)
a = -miu/2/xi;
p = a*(1 - norm(e)^2);
else
p = h^2/miu;
a = inf;
end

 

function deltav = test(xv)

Mr = 227939186e3;
mius = 132712440018e9; %sun miu
Vmars = sqrt(mius/Mr);

dist = zeros(1,length(xv)); %calc distance from earth
for ii = 1:length(xv)
dist(ii) = norm(xv(ii,1:2));
end
index = find(dist>=Mr);
if(isempty(index))
disp('not reach mars orbit');
end
index = index(1);

alpha = atan2(xv(index,4),xv(index,3));
beta = alpha + pi/2;
deltav = norm(xv(index,3:4)-[Vmars*sin(beta),Vmars*cos(beta)]);