Walker-Delta Satellite Constellation for Earth

Observation:
SPE 510 Midterm Report
Kelvin Loh (20114534)
Room 3333, Department Of Aerospace Engineering
KAIST
[email protected]
April 9, 2012
Abstract
A suite of MATLAB functions are developed to analyze the Walker-Delta Con-
stellation. Since this is one of the more popular constellations, two cases are
also being thought out to determine if the constellation pattern can be a so-
lution for the cases that require extensive Earth observation. It is shown that
the constellation is sensitive to altitude changes.
1 Introduction
A constellation is a set of satellites distributed over space intended to work together
to achieve common objectives. The Walker-Delta Constellation is one such popular
constellation design since it is the most symmetric. In order to study the performance
of this type of constellation, a set of computational tools have to be developed. The
theoretical foundations are important as it directly aects the accuracy of the anal-
ysis. The functions are tested and the results of the tests are also reported in this
report. Additional models are also included in the code such as low-thrust maneu-
vers tangential to the spacecraft ight path, but they are only documented in the
comments of the MATLAB functions.
2 K. Loh
2 Theoretical Models
This section covers the theoretical basis for the codes used to evaluate the perfor-
mances of the Walker-Delta constellation to cover a target on Earth for observation.
Proofs of the models used are not provided in this report.
2.1 Equations of relative motion
The two-body equation of motion, namely, Newtons second law and his universal
law of gravitation form the starting point for any study of astrodynamics. Equation
1 describes the motion of a body (satellite) of negligible mass with respect to another
body (Earth) under the inuence of a gravitational force. This then assumes that an
inertial frame of reference can be centered at the barycenter of the massive body.

r =

r
2
r
r
(1)
In Newtons law of gravitation, the gravitational force can be expressed in terms
of the gradient of a scalar gravitational potential function of a body. For a
spherical body

F = =
_

r
_
(2)
For a non-spherical body, can contain additional terms which can be used to
determine

F, and subsequently

r. As the inertial reference frame is centered at the
center of the massive body, Equation 1 is also known as the fundamental equation
of relative two-body motion. The MATLAB function orbit.m uses this equation to
determine the satellite motion.
2.2 Keplerian Orbit Propagation
Keplers laws are derived from Equation 1. The orbital position as a function of
time can be decribed by the mean anomaly, M
e
. The propagation routines used
in the codes ground track2.m and satellite constellation.m are based on solving the
Keplers equation (Equation 3) for each time step.
M
e
= E e sin (E) (3)
where E is the eccentric anomaly, and e is the eccentricity.
Walker-Delta Satellite Constellation for Earth Observation 3
For elliptical orbits, it can be shown [1] that
M
e
=

2
h
3
_
1 e
2
_3
2
t (4)
So, Equation 4 becomes
M
e
=
2
T
t (5)
With this, it can be seen that numerically solving the Keplers equation for each
time step and directly using an Euler forward rst-order time-stepping scheme is faster
than directly integrating Equation 1 using a Runge-Kutta (4th-order or higher) xed
time-step scheme. Therefore, this method is primarily used to propagate the satellite
constellation for this project.
2.3 Orbit Perturbations
Equation 1 assumes that the revolving body experiences forces only from the central
body. However, this assumption does not hold for actual orbits. In reality, there are
many perturbations to the orbital motion. The equation can then be generalized to
include all the other perturbation forces acting on the body as described by Equation
6.

r =

r
2
r
r
+a
p
(6)
As the Earth is not a perfect sphere, the gravitational potential will deviate from
that of a perfect sphere. Due to this eect, the orbital motion will encounter pertur-
bations. This project considers the eects of Earth Gravity Harmonics mainly, the
J
2
perturbation on the secular rates of the right ascension of ascending node (),
the argument of perigee (
p
), and a small correction to the mean motion of the orbit
(

M
e
). The rates are computed by the following equations:

=
3
2
J
2
R
2
E
p
2
ncos i (7)

p
=
3
2
J
2
R
2
E
p
2
n
_
2
5
2
sin
2
i
_
(8)

M
e
= n =
_

a
3
0
_
1 +
3
2
J
2
R
2
E
p
2
_
1
3
2
sin
2
i
_

1 e
2
_
(9)
4 K. Loh
with the symbols taking their usual notations. Atmospheric drag perturbation eects
was not included in the simulations since it is assumed that the main propulsion sys-
tem onboard a satellite will perform the necessary station-keeping maneuvers during
the mission.
For the special perturbation method employed in orbit.m, Cowells method [2]
was used. From Equation 2, it is seen that the gravitational force acting on the
revolving body can be described by the gradient of a potential function. From this,
the potential theory provides a clear picture of the gravity harmonics of the Earth.
The potential function for the primary body due to the J
2
perturbation is described
by Equation 10 in the inertial Cartesian coordinate frame (ECI).
=

r
1
2
J
2
_
R
E
r
_
2
_
1 3
_
z
r
_
2
_
(10)
where r =
_
x
2
+y
2
+z
2
It can then be derived from Equation 10 and Equation 2 that
a
J
2
=
3
2
J
2
R
2
E
r
4
_
x
r
_
1 5
_
z
r
_
2
_
,
y
r
_
1 5
_
z
r
_
2
_
,
z
r
_
3 5
_
z
r
_
2
__

(11)
The acceleration a
J
2
is then added to the right-hand side of Equation 6 as part of the
perturbation acceleration.
2.4 Earth Coverage
This section will only reference gures from the text used [3] to compute the Earth
target coverage by a satellite. Also, the programs assume that for coverage purposes,
the Earth is a perfect sphere. An area access coordinate frame on the Earths surface
is dened corresponding to the nadir coordinate frame of the spacecraft. Figure 9-2
of the text [3] shows the denition of the access area coordinate frame. The trans-
formation between latitude and longitude and access area coordinates is relatively
straightforward.
For the programs, the maximum Earth Central Angle is given by,

0
=
R
E
R
E
+H
(12)
However, to account for the optical sensor resolution, the Earth central angle is
limited by the resolution requirements (Figure 9-3 of the text [3]). From the Rayleigh
Walker-Delta Satellite Constellation for Earth Observation 5
criterion and cosine rule,
sin =
x
res
2h

= 1.220
WL
D
a
(13)
cos = 1
(h

)
2
h
2
2R
E
(R
E
+H)
(14)
The transformation from the latitude and longitude of the target (P) and sub-
satellite point (SSP) (Figure 9-4 of the text [3])are as follows:
Lat

P
= 90

Lat
P
, Lat

SSP
= 90

Lat
SSP
(15)
L = Long
SSP
Long
P
(16)

r
P
= cos
1
[cos Lat

P
cos Lat

SSP
+ sin Lat

P
sin Lat

SSP
cos L] (17)
The target will be observable from the spacecraft optical payload if
r
P
min(
0
, )
2.5 Walker-Delta constellation
The most symmetric of the satellite patterns is the Walker Constellation. The Walker-
Delta Pattern contains a total of T satellites with S satellites evenly distributed
in each of P orbit planes. All of the orbit planes are assumed to be at the same
inclination, i, relative to the equator. The ascending nodes of the P orbit planes
in Walker patterns are uniformly distributed around the equator at intervals of
360

P
.
Within each orbital plane, the S satellites are uniformly distributed at intervals of
360

S
. The phase dierence, , in a constellation is dened as the angle in the
direction of motion from the ascending node to the nearest satellite at a time when
a satellite in the next most westerly plane is at its ascending node. must be an
integral multiple, F, of
360

T
, where 0 < F P 1.
Figure 1 shows an example of a Walker-Delta pattern generated by the program
ground track2.m with the satellites at their initial positions. The green cube shows
the seed/rst satellite.
Figure 2 shows the classical orbit elements for the satellites in the 45:4/2/1 Walker-
Delta constellation as generated by the ground track2.m program.
3 General Description of the Codes
Several MATLAB functions (mainly, optim main.m, ground track2.m, orbit.m, and
test t.m) were written to analyze the behaviour of the Walker-Delta constellation.
6 K. Loh
Figure 1: Walker-Delta Pattern 45:4/2/1 at 10,000 km altitude
Figure 2: The initial Walker-Delta Pattern in Classical Orbit Elements
Walker-Delta Satellite Constellation for Earth Observation 7
The ground track2.m function was written as a tool to test the orbit propagator
(satellite constellation.m) for a perturbed Keplerian orbit used in optim main.m as
well as to visualize the nal design parameter chosen for the observation of KAIST.
In general the basic structure of the codes are similar to each other. The rst sec-
tion will be the inputs of the simulation, second function will initialize the spacecraft
state vectors in the ECI (SV) or in the Keplerian orbit elements (COE). After that
each satellite in the constellation are propagated along the same rates, hence, for
the constellation, only the initial condition is important. After the SVs or COEs
are known at each time step from the propagator, the position of the satellites
are then converted into the Earth latitude and longitude coordinate frame. In the
nd ra and dec function, this is being done. This makes it easier to obtain the
subsatellite point latitudes and longitudes for coverage determination.
The appendix lists the important codes. The optim main.m code is a wrapper code
which can brute force calculate potentially all the design variables of the Walker-Delta
constellation. It uses the MATLAB parallel toolbox for independent design variable
computations.
4 Validation cases
There currently are no test cases to be compared against. However, the author
believes that if the individual functions are working and tested, the results after the
functions are combined can be reliable.
4.1 Initialization and propagation
The rst test is the satellite initialization subroutine. The output in Figure 2 shows
that the initialization was done correctly. The orbit was also propagated correctly
with all the outputs showing that the satellites were uniformly propagated. The
nal output is shown in Figure 3. Table 1 presents the results of a single satellite
using dierent propagator models. One is run using orbit.m and the other using
ground track2.m. As can be seen, the results show not a large dierence between
the RK5 and Kepler propagator model, with the exception that the Kepler model is
faster to run but the disadvantage is that it assumes an elliptical orbit.
8 K. Loh
Figure 3: Final output
Table 1: Results of FOM for single satellite using dierent propagator models
Method of propagation Average altitude (km) % Coverage Mean response time (s)
RK5 - Newton 392.75 0.255 27987
Kepler 400 0.301 27667
4.2 Earth coverage
Unfortunately, there is no known standard test case to be run for the calculation of
the earth coverage. However, another MATLAB function orbit2.m was written to test
the earth coverage calculations and to verify them. The function orbit2.m calculates
the Earth central angle
rp
from target on the surface to the subsatellite point using
both Equations 17 and 18 for comparison. Equation 18 is dened as follows:
cos
rp
=
r r
t
|r||r
t
|
(18)
where r is the position vector of the satellite from the center of the Earth in ECI
frame, and r
t
is the position vector of the target on the surface of the Earth from the
center of the Earth in ECI frame.
Figure 4 shows the path of the target on the surface of the Earth with respect to
the ECI frame as calculated by the orbit2.m function. It shows the correct path as
predicted since the Earth is modelled to spin only in the z-axis. It is then assured
that the target position vector is correct. Also, it is known that the satellite position
vector is also correct. Therefore, Equation 18 can be used as a verication for the
method used in the other functions.
A test orbit is used for method validation. The target is KAIST (36.372 N, 127.363
E), with inclination of 60

, and altitude of 400 km. Figure 5 shows the observation

histogram as calculated for both methods using Equation 17 (Method 1) and Equation
18 (Method 2). The percentage coverage estimated from both methods are identical
at 0.486%, mean response time at 20575 s, and maximum response time at 53768 s. It
shows that they are identical as expected, hence, the Earth coverage functions used in
Walker-Delta Satellite Constellation for Earth Observation 9
Figure 4: Path of target with respect to ECI frame (White dotted line)
all the other functions (orbit.m, satellite constellation.m, and ground track2.m) are
trusted to be correct.
Figure 5: Observation histogram for both methods
4.3 Figure-of-Merit
The calculation of the coverage FOM by the subroutine test t.m is particularly im-
portant to establish the performance of the design of the constellation. Hence the
percentage coverage and mean response time calculations must be validated. For the
10 K. Loh
test case, the observation patterns for case A is reproduced with varying number of
timesteps for the same length. The rst is that shown by the text [3] with dt = 1
time unit. The second case gives a dt = 0.1 time unit. The third case gives dt = 0.01
time unit.
Figure 6: First case A (dt = 1)
The percent coverage and mean response time matches well with the FOM given
out by the text (Page 485 in [3]). When the time steps are rened, however, the
mean response time diers by quite a signicant value. The second case A with a
ner time step of 0.1 time unit, gives a mean response time of 0.32 time unit instead
of 0.5 as quoted by the text. Figures 6 to 8 show the case A observation histogram
with increasingly ner timesteps.
Figure 7: Second case A (dt = 0.1)
The third case A with the nest time step of 0.01 time unit gave a mean response
time of 0.302 time unit with a percent coverage of 60%. Figure 9 shows the results
Walker-Delta Satellite Constellation for Earth Observation 11
Figure 8: Third case A (dt = 0.01)
from the test t.m for the mean response time as a function of dt. It shows that
the actual mean response time as dt gets closer to 0 is 0.3 time unit instead of the
given 0.5 by the text. However, after reading up on the text and checking the code
several times, the author believes that the way to calculate the mean response time
by the subroutine test t.m was correct and attributes the dierence in the results by
temporal renement issues.
Figure 9: Mean response time as a function of dt
12 K. Loh
5 Results and Discussion
5.1 Case 1 - Observation of the city of Patna
Patna is the capital of the Indian state Bihar. It is situated at 25.611 N, and 85.144
E. It is a growing state in India therefore, there is a need for remote sensing from
space to help city planners in building the city. As such, this case serves as a rst
test for both the code and also the idea of having remote sensing satellites to aid in
city planning. The requirements are laid out such that a ground resolution is no more
than 5 m. Also, the number of satellites should be no more than 10.
The design chosen to test is a Walker-Delta constellation of 26:10/5/1. The al-
titude is given at 1000 km. The percent coverage obtained for a 24 hour simulation
was 70.03% and the mean response time was shown to be 81.99 s with a maximum
response time of 759.8 s. Figures 10 to 12 shows the constellation conguration along
with the orbit, ground track, and the observation histogram.
Figure 10: Case 1 - Initial constellation conguration
5.2 Case 2 - Observation of KAIST
KAIST is South Koreas rst research oriented science and engineering institution.
Naturally, as South Koreas foremost center for strategic research and development,
the safety of the institute is of paramount importance. In order to help safeguard this
institution, it is required that a constellation of remote sensing satellites be placed in
Walker-Delta Satellite Constellation for Earth Observation 13
Figure 11: Case 1 - Ground track for 1 revolution
Figure 12: Case 1 - Observation histogram for full 24 hours
14 K. Loh
orbit with a ground resolution of 1.1 m, and the camera payload aperture diameter
at 1 m. The inclination of the orbit has to be 60

and the total number of satellites

must not be more than 30 due to cost considerations. The percent coverage is of
utmost priority and has to be maximized for a duration of 1 day.
Figure 13: as a function of altitude, h, and ground resolution, x
res
Figure 13 shows the contour plot of the Earth Central Angle limitation imposed
by Rayleighs Criterion as a function of altitude and ground resolution requirements.
For the aperture diameter of 1 m, the limited Earth central angle is at 5

with an
altitude of 700 km. This also imposes an altitude limitation in which the satellite
altitudes must not go above 1000 km.
For the case of observing KAIST, the optim main.m code was used to determine
the design variables which gives the maximum percent coverage. The brute force
searching technique was used. Figure 14 provides an indication of the percent coverage
variation with respect to the total number of satellites (T) parameter and the altitude.
For this case, the number of planes (P) is set at the maximum which is equal to T.
The discontinuous looking pattern is very likely due to the discrete nature of the
variable T, and also, the simulation for observing the ground, the time step might be
too large in that it might have skipped a few points on the ground, thereby causing
the percent coverage plot to be rather irregular. However, a pattern can be seen,
which shows that the percent coverage increases with decreasing altitude. This is
probably due to the increasing limited Earth Central Angle imposed by the Rayleigh
criterion. Of particular interest is shown in Figure 15 that increasing or decreasing
the number of planes with a maximum number of satellites (T = 30) does not change
Walker-Delta Satellite Constellation for Earth Observation 15
the percent coverage by a signicant amount. What is clear from the plot is that the
largest dierence comes again yet from changing altitude. This can only be explained
by performance plateaus which is aected by altitude and inclination as described by
the text [3].
Figure 14: Contour map of percent coverage as a function of altitude and number of
satellites
An optimization wrapper, optim main2.m was used to obtain the maximum Earth
coverage with respect to the number of planes (P), relative phasing parameter (F),
and altitude. As expected the maximum percentage coverage happens in an altitude
of 100 km, however, this is too low when drag is considered. Hence, an altitude of
400 km is chosen from the data set and the maximum set of the p cover1 variable
corresponding to the altitude of 400 km is chosen. Using the information gained, the
constellation 60:30/30/15 is chosen as the nal design point to maximize the percent
coverage to observe KAIST. The constellation will have an altitude of 400 km and
the initial longitude of the rst satellite at t=0 is at 0

. The estimated performance

of this constellation design is a percentage coverage of 14.72%, a mean response time
of 509.26 s, and a maximum response time of 1239.71 s. Figures 16 to 18 shows the
conguration, ground track, and observation histogram respectively for the satellite
constellation for a time length of 1 orbital period.
The results show that the ground resolution, and inclination requirements have to
be relaxed to attain a more feasible solution. As an additional feature, an animation
was also produced of the constellation chosen. The animation MATLAB function
16 K. Loh
Figure 15: Contour map of percent coverage as a function of altitude and number of
planes
Figure 16: Initial satellite constellation (60:30/30/15 conguration)
Walker-Delta Satellite Constellation for Earth Observation 17
Figure 17: Ground track for the constellation after 1 period
Figure 18: Observation histogram for the constellation observing KAIST in 24 hours
18 K. Loh
is written as ground track anim.m. The animation can be accessed via a youtube
link http://youtu.be/w1eb37tOei4.
6 Summary
A suite of MATLAB functions was developed to analyze the behaviour of Walker-
Delta Constellations. Two propagator models were compared, and both models con-
sidered the oblateness of the Earth (J
2
perturbation eects). The functions were
validated and the functions were used to analyze the performance of such constella-
tions in observing primarily, KAIST, and Patna. The results show that the design
requirements are quite tight for the KAIST case while the imaginary Patna case seem
to be feasible.
Acknowledgments
The author is grateful to Prof. Ahn Jaemyung for his lectures and out-of-class tech-
nical discussion which were necessary to complete this project. The author acknowl-
edges the use of the algorithms as structured by one of the text used [1]. The author
would also like to thank Miss Shruti Pavagadhi for several out-of-class discussions
which gave the necessary motivations and support for the completion of this project.
References
[1] Curtis, H. Orbital Mechanics for Engineering Students, 2nd Ed. (Elsevier Ltd,
Oxford, UK, 2010).
[2] Chobotov, V. Orbital Mechanics, 2nd Ed. (AIAA Education Series, Reston, VA,
USA 1996).
[3] Wertz, J. Mission Geometry: Orbit and Constellation Design and Management,
1st Ed. (Space Technology Library, Microcosm Press, CA, USA, 2001).
A Listing of code - optim main2.m
1 cl ear al l ; cl ose al l ; cl c ;
2 matl abpool ( 4) ; %Comment out i f no p a r a l l e l t ool b ox i n s t a l l e d
Walker-Delta Satellite Constellation for Earth Observation 19
3 deg = pi /180;
4
5 mu = 398600;
6 J2 = 0. 00108263;
7 Re = 6378;
8 we = (2 pi + 2pi /365. 26) /(243600) ;
9 %Walker c o ns t e l l a t i o n des i gn v a r i a b l e s
10 i ncl w = 60deg ; %I nc l i na t i o n of t he pl anes
11 Nsat = 30; %t t o t a l number of s a t e l l i t e s Common sense t e l l s us t hat s at
+ 1 >= s at ! So max s at == max coverage v a r i a b l e
12 Nplane = f ac t or 2 ( Nsat ) ;
13
14 for i = 1: length( Nplane )
15 f 1 = 0:max( Nplane ) ; %f phase mu l t i p l i e r
16 f ( i , : ) = f 1 ;
17 end
18 %. . . Seed/ Fi r s t s a t e l l i t e paramet ers
19 Al t = 100: 50: 750;
20 hP1 = Al t ;
21 hA1 = Al t ;
22 TAo1 = 0deg ;
23 Woi = 0deg ;
24 wpo1 = 0deg ;
25
26 %Case Si mul at i on l e ng t h and t i mes t eps
27 dt = 20;
28 to = 0;
29 t f = 1243600; % Ts i n assi gnment paper
30
31 %Target l a t i t u d e and l ong i t ude
32 l atT = 36. 372 deg ; %36. 372 deg ;
33 longT = 127. 363 deg ; %127. 363 deg ;
34
35 %S a t e l l i t e payl oad paramet ers
36 Da = 1;
37 WL = 5e 7;
38 x r e s = 1 . 1 ;
39 KA = 20626. 4806; % f or IAA area i n deg 2
40
41 for k = 1: length( Nplane )
42 for i i = 1: length( f ( k , : ) )
43 i f ( f ( k , i i ) <= max( f ( k , : ) ) )
20 K. Loh
44 par f or i = 1: length( Al t ) %Change par f or t o f or i f no
p a r a l l e l t ool b ox i s i n s t a l l e d
45 [ t1 SEE1 ] = s a t e l l i t e c o n s t e l l a t i o n (mu, J2 , Re , we , i ncl w ,
Nsat , Nplane ( k) , f ( k , i i ) , hP1( i ) , hA1( i ) , TAo1, Woi , wpo1 ,
dt , to , t f , l atT , longT , Da,WL, x r es ,KA) ;
46 [ p cover 1 ( k , i i , i ) mean response1 ( k , i i , i ) max response1 ( k
, i i , i ) ] = t e s t t (SEE1 , t1 ) ;
47 drawnow;
48 end
49 el se cont i nue
50 end
51 end
52 end
53
54 matl abpool cl ose ; %Comment out i f no p a r a l l e l t ool b ox i n s t a l l e d
55
56 save r e s ul t s maxi s i ng P c ove r . mat
57
58 %Show i ndex f or maximum P cover
59 [ i x i y i z ] = i nd2sub ( si ze ( p cover 1 ) , find ( p cover 1 == max(max(max(
p cover 1 ) ) ) ) ) ;
60 f pri ntf ( The optimum s o l ut i o n i s %g:%d/%d/%d at an a l t i t ude of %g km
wi th a cover age of %g per cent \n , i ncl w/deg , Nsat , Nplane ( i x ) , f ( i x , i y )
, Al t ( i z ) , p cover 1 ( i x , i y , i z ) 100) ;
61
62 %Se l e c t a l t i t u d e at 400 km
63 [ Aind ] = i nd2sub ( si ze ( Al t ) , find ( Al t == 400) ) ;
64 [ i i x i i y ] = i nd2sub ( si ze ( p cover 1 ( : , : , Aind) ) , find ( p cover 1 ( : , : , Aind) ==
max(max( p cover 1 ( : , : , Aind) ) ) ) ) ;
65 f pri ntf ( The optimum s o l ut i o n i s %g:%d/%d/%d at an a l t i t ude of %g km
wi th a cover age of %g per cent \n , i ncl w/deg , Nsat , Nplane ( i i x ) , f ( i i x ,
i i y ) , Al t ( Aind) , p cover 1 ( i i x , i i y , Aind) 100) ;
66
67 % p l o t ( Nplane , p cover1 ( : , 1) 100 , ok ) ; x l a b e l ( Number of pl anes ) ; y l a b e l
( Percent age coverage %) ;
B Listing of code - satellite constellation.m
1 %
2 function [ t X] = s a t e l l i t e c o n s t e l l a t i o n (mu, J2 , Re , we , i ncl w , Nsat , Nplane , f
, hP1 , hA1 , TAo1, Woi , wpo1 , dt , to , t f , l atT , longT , Da,WL, x r es ,KA)
Walker-Delta Satellite Constellation for Earth Observation 21
3 %
4
5 %c l e ar a l l ; c l os e a l l ; c l c
6 global ra dec n cur ves RA Dec
7
8 deg = pi /180;
9 hours = 3600;
10 %Cal cul at e Pat t ern Unit f or Walker c o ns t e l l a t i o n
11 PU = 360 deg/Nsat ;
12 %Cal cul at e Number of s a t e l l i t e s i n a pl ane
13 NSsat = Nsat /Nplane ;
14 %Angl e d i v i s i o n bet ween s a t e l l i t e s i n t he pl ane
15 beta = PUNplane ;
16 %RAAN s paci ng
17 dRAAN = PUNSsat ;
18 %. . . Seed/ Fi r s t s a t e l l i t e paramet ers
19 hP( 1 : Nsat ) = hP1 ;
20 hA( 1 : Nsat ) = hA1 ;
21 rP1 = hP1 + Re ; rP( 1 : Nsat ) = rP1 ;
22 rA1 = hA1 + Re ; rA( 1 : Nsat ) = rA1 ;
23 i nc l 1 = i ncl w ; i n c l ( 1 : Nsat ) =i nc l 1 ;
24 TAo( 1 : Nsat ) = TAo1;
25 Wo1( 1 : Nsat ) = Woi ;
26 wp( 1 : Nsat ) = wpo1 ;
27
28 f pri ntf ( Sol vi ng EOM. . . )
29 for i =1: Nsat
30 pi ndex ( i ) = cei l ( i /NSsat ) ;
31 Spi ndex ( i ) = i ( pi ndex ( i ) 1) NSsat 1;
32 f s at i nde x ( i ) = abs(1sign( Spi ndex ( i ) ) ) ;
33 TAo( i ) = TAo1 + Spi ndex ( i ) beta + ( pi ndex ( i ) 1) f PU;
34 W( i ) = Wo1( i ) + ( pi ndex ( i ) 1) dRAAN;
35 end
36
37 %Cal cul at e Earth c e nt r al angl e from payl oad paramet ers
38 h pri me = Da x r e s 1e 3/(2. 440WL) ; %From Rayl ei gh o p t i c a l r e s o l ut i o n (
km)
39 %. . . End dat a de c l ar at i on
40
41 %. . . Compute t he r at e s of node r e g r e s s i on and pe r i g e e advance
42 a = ( rA + rP) /2; %1: Nsat
43 T = 2pi /sqrt (mu) . a . ( 3/2) ; %1: Nsat
22 K. Loh
44 e = ( rA rP) . / ( rA + rP) ; %1: Nsat
45 h = sqrt (mu. a . ( 1 e . 2) ) ; %1: Nsat
46 Eo = 2atan( tan(TAo/2) . sqrt ((1e ) . /(1+e ) ) ) ; %1: Nsat
47 M = Eo e . si n (Eo) ; %1: Nsat
48 P = a.(1 e . 2) ;
49 f ac = 3/2J2 (Re . /P) . 2 ;
50 Mdot = (2 pi . /T) . ( 1 + f ac . sqrt(1e . 2) . (1 (3/2) . si n ( i n c l ) . 2) ) ; %
From Chobotov
51 Wdot = Mdot . f ac . cos ( i n c l ) ;
52 wpdot = f ac . Mdot. ( 5/2 si n ( i n c l ) . 2 + 2) ;
53
54
55 % f or i =1: Nsat
56 % f p r i n t f ( \n Sat # = %d , h = %g km, i = %g deg , Wo = %g deg ,
wpo = %g deg , TAo = %g deg , Mo = %g deg \n , i , hP( i ) , i nc l ( i ) /deg ,
W( i ) /deg , wp( i ) /deg , TAo( i ) /deg , M( i ) /deg ) ;
57 % end
58 ti mes=to : dt : t f ;
59 ra ( 1 : length( ti mes ) , 1: Nsat ) = 0;
60 dec ( 1 : length( ti mes ) , 1: Nsat ) = 0;
61 t het a = 0;
62 r r e l mag ( 1 : length( ti mes ) , 1: Nsat ) = 0;
63 for i = 1: length( ti mes )
64 t ( i ) = ti mes ( i ) ;
65 for i i = 1: Nsat
66 M( i i ) = M( i i ) + Mdot( i i ) dt ;
67 E( i i ) = kepl er E ( e ( i i ) ,M( i i ) ) ;
68 TA( i i ) = 2atan( tan(E( i i ) /2) sqrt ((1+e ( i i ) ) /(1e ( i i ) ) ) ) ;
69 r = h( i i ) 2/mu/(1 + e ( i i ) cos (TA( i i ) ) ) [ cos (TA( i i ) ) si n (TA( i i ) ) 0 ] ;
70
71 W( i i ) = W( i i ) + Wdot( i i ) dt ;
72 wp( i i ) = wp( i i ) + wpdot ( i i ) dt ;
73 Rone = R3(W( i i ) ) ;
74 Rtwo = R1( i n c l ( i i ) ) ;
75 Rthree = R3(wp( i i ) ) ;
76 QxX = ( Rthree RtwoRone ) ;
77 R = QxX r ;
78
79 t het a = we( t ( i ) to ) ;
80 Q = R3( t het a ) ;
81 r r e l = QR;
82 r r e l mag ( i , i i ) = norm( r r e l ) ;
Walker-Delta Satellite Constellation for Earth Observation 23
83 [ al pha de l t a ] = r a and de c f r om r ( r r e l ) ;
84
85 ra ( i , i i ) = al pha ;
86 dec ( i , i i ) = de l t a ;
87 end
88 end
89
90 SEE( 1 : length( ti mes ) , 1: Nsat+1) = 0;
91 for i = 1: length( ti mes )
92 t i = ti mes ( i ) ;
93 for i i = 1: Nsat
94 a l t = r r e l mag ( i , i i ) Re ;
95 clambda = 1( h pri me 2 a l t 2) /(2Re(Re + a l t ) ) ; %Cos l ambda
96 IAA = KA. ( 1 clambda ) ;
97 l at P pr i me = pi /2 l atT ; %Transf ormat i on from l a t l ong t o
acces s area coor di nat es ( Sect i on 9. 1 Mi ssi on Geometry (
Wertz ) )
98 l atSSP pri me = pi /2 dec ( i , i i ) deg ;
99 del t a L = ra ( i , i i ) deg longT ;
100 lambda ( i , i i ) = acosd ( clambda ) ;
101 lambda0 ( i , i i ) = acosd (Re/(Re + a l t ) ) ;
102 rlambda ( i , i i ) = acosd ( cos ( l at P pr i me ) cos ( l atSSP pri me ) + si n (
l at P pr i me ) si n ( l atSSP pri me ) cos ( del t a L ) ) ; %Angul ar
Di st ance from t a r g e t t o s u b s a t e l l i t e poi nt , Di t t o Eqn 99
103 i f ( rlambda ( i , i i ) > 180)
104 rlambda ( i , i i ) = 360 rlambda ( i , i i ) ;
105 end
106 i f ( rlambda ( i , i i ) <= min( lambda ( i , i i ) , lambda0 ( i , i i ) ) )
107 SEE( i , i i ) = 1; %I f t a r g e t i s wi t hi n t he minimum of e i t h e r
l ambda or lambda0 , then , i t i s covered by t he i i t h
s a t e l l i t e
108 end
109 end
110 SEE( i , Nsat+1) = sign(sum(SEE( i , 1 : Nsat ) ) ) ; %As l ong as t her e i s one
s a t e l l i t e covered , t hen t he t a r g e t i s covered by t he
c o ns t e l l a t i o n
111 end
112 % f or m s epar at e cur ves
113 % pl ot g r ound t r ac k
114 % p r i n t o r b i t a l d a t a
115 % f i g ur e (3)
116 % p l o t ( t /hours , SEE( : , Nsat +1) , k ) ;
24 K. Loh
117 X = SEE( : , Nsat+1) ;
118 %[ p cover mean response max response ] = t e s t t (SEE( : , Nsat +1) , t ) ;
119
120 % f i g ur e (4)
121 % p l o t ( t /hours , al t , k ) ;
122 %ret urn
123 f pri ntf ( Done ! \n )
124 end %gr ound t r ack
125 %
C Listing of code - ground track2.m
1 %
2 function gr ound t r ack2
3 %
4
5 cl ear al l ; cl ose al l ; cl c
6 global ra dec n cur ves RA Dec
7 %mat l abpool (4) %Only used i f p a r a l l e l t ool b ox i s i n s t a l l e d
8 deg = pi /180;
9 hours = 3600;
10 mu = 398600;
11 J2 = 0. 00108263;
12 Re = 6378;
13 we = (2 pi + 2pi /365. 26) /(243600) ;
14 %Walker c o ns t e l l a t i o n des i gn v a r i a b l e s
15 i ncl w = 60deg ; %I nc l i na t i o n of t he pl anes
16 Nsat = 30; %t t o t a l number of s a t e l l i t e s
17 Nplane = 30; %p t o t a l number of pl anes
18 f = 15; %f phase mu l t i p l i e r
19 %Cal cul at e Pat t ern Unit f or Walker c o ns t e l l a t i o n
20 PU = 360 deg/Nsat ;
21 %Cal cul at e Number of s a t e l l i t e s i n a pl ane
22 NSsat = Nsat /Nplane ;
23 %Angl e d i v i s i o n bet ween s a t e l l i t e s i n t he pl ane
24 beta = PUNplane ;
25 %RAAN s paci ng
26 dRAAN = PUNSsat ;
27 %Al t i t ude ( i n km)
28 H = 400;
29 %Longi t ude of t he f i r s t s a t e l l i t e i n t he f i r s t pl ane at t =0
Walker-Delta Satellite Constellation for Earth Observation 25
30 Lon0 = 0;
31 %. . . Seed/ Fi r s t s a t e l l i t e paramet ers
32 hP1 = H; hP( 1 : Nsat ) = hP1 ;
33 hA1 = H; hA( 1 : Nsat ) = hA1 ;
34 rP1 = hP1 + Re ; rP( 1 : Nsat ) = rP1 ;
35 rA1 = hA1 + Re ; rA( 1 : Nsat ) = rA1 ;
36 i nc l 1 = i ncl w ; i n c l ( 1 : Nsat ) =i nc l 1 ;
37 TAo1 = 0deg ; TAo( 1 : Nsat ) = TAo1;
38 Wo1 = Lon0deg ; Wo1( 1 : Nsat ) = Wo1;
39 wpo1 = 0deg ; wp( 1 : Nsat ) = wpo1 ;
40 i n i t i a l i z e wa l k e r
41 dt = 20;
42 %Target l a t i t u d e and l ong i t ude ( on t he ground )
43 l atT = 36. 372 deg ; %36. 372 deg ;%25. 611 deg ;
44 longT = 127. 363 deg ;%85. 144 deg;%
45 %S a t e l l i t e payl oad paramet ers
46 Da = 1;
47 WL = 5e 7;
48 x r e s = 1 . 1 ;
49 KA = 20626. 4806; % f or IAA area i n deg 2
50 %Cal cul at e Earth c e nt r al angl e from payl oad paramet ers
51 h pri me = Da x r e s 1e 3/(2. 440WL) ; %From Rayl ei gh o p t i c a l r e s o l ut i o n (
km)
52 %. . . End dat a de c l ar at i on
53
54 %. . . Compute t he r at e s of node r e g r e s s i on and pe r i g e e advance
55 a = ( rA + rP) /2; %1: Nsat
56 T = 2pi /sqrt (mu) . a . ( 3/2) ; %1: Nsat
57 e = ( rA rP) . / ( rA + rP) ; %1: Nsat
58 h = sqrt (mu. a . ( 1 e . 2) ) ; %1: Nsat
59 Eo = 2atan( tan(TAo/2) . sqrt ((1e ) . /(1+e ) ) ) ; %1: Nsat
60 M = Eo e . si n (Eo) ; %1: Nsat
61 coe0 = [ h e W / deg i nc l / deg wp / deg TAo / deg ] ;
62 to = 0;
63 t f = T; %1243600;
64 ti mes=to : dt : t f ;
65 P = a.(1 e . 2) ;
66 f ac = 3/2J2 (Re . /P) . 2 ;
67 Mdot = (2 pi . /T) . ( 1 + f ac . sqrt(1e . 2) . (1 (3/2) . si n ( i n c l ) . 2) ) ;
68 Wdot = Mdot . f ac . cos ( i n c l ) ;
69 wpdot = f ac . Mdot. ( 5/2 si n ( i n c l ) . 2 + 2) ;
70
26 K. Loh
71 for i =1: Nsat
72 [ R0 V0] = s v f r om coe ( coe0 ( i , : ) ,mu) ;
73 Rpl ot ( 1 , : , i ) = [ R0 V0 ] ;
74 f pri ntf ( \n Sat # = %d , h = %g km, i = %g deg , Wo = %g deg , wpo = %g
deg , TAo = %g deg , Mo = %g deg \n , i , hP( i ) , i n c l ( i ) /deg , W( i ) /
deg , wp( i ) /deg , TAo( i ) /deg , M( i ) /deg ) ;
75 end
76 f i nd r a and de c
77 fi gure ( 2)
78 ground map
79 for i i = 1: Nsat
80 f or m s e par at e c ur ve s
81 pl ot gr ound t r ac k
82 end
83 % p r i n t o r b i t a l d a t a
84 fi gure ( 1)
85 output
86 fi gure ( 3)
87 plot ( t 60/ hours , SEE( : , Nsat+1) , k ) ;
88 axis ( [ t ( 1) 60/ hours t (end) 60/ hours 0 1 . 1 ] ) ;
89 xlabel ( Time ( Mins ) ) ; ylabel ( Observati on ) ;
90 f pri ntf ( Fi nal o r bi t a l el ements ) ;
91 for i =1: Nsat
92 f pri ntf ( \n Sat # = %d , h = %g km, i = %g deg , Wo = %g deg , wpo = %g
deg , TAo = %g deg , Mo = %g deg \n , i , hP( i ) , i n c l ( i ) /deg , W( i ) /
deg , wp( i ) /deg , TA( i ) /deg , M( i ) /deg ) ;
93 end
94 [ p cover mean response max response ] = t e s t t (SEE( : , Nsat+1) , t ) ;
95 f pri ntf ( \n Percentage coverage = %g per cent | Mean r es pons e ti me = %g
seconds | Max r es pons e ti me = %g seconds \n , p cover 100 ,
mean response , max response )
96 % f i g ur e (4)
97 % p l o t ( t /hours , al t , k ) ;
98 return
99
100 function i n i t i a l i z e wa l k e r
101 f pri ntf ( I n i t i a l s a t e l l i t e c o ns t e l l a t i o n pr o pe r t i e s \n )
102 for i =1: Nsat
103 pi ndex ( i ) = cei l ( i /NSsat ) ;
104 Spi ndex ( i ) = i ( pi ndex ( i ) 1) NSsat 1;
105 f s at i nde x ( i ) = abs(1sign( Spi ndex ( i ) ) ) ;
106 TAo( i ) = TAo1 + Spi ndex ( i ) beta + ( pi ndex ( i ) 1) f PU;
Walker-Delta Satellite Constellation for Earth Observation 27
107 W( i ) = Wo1( i ) + ( pi ndex ( i ) 1) dRAAN;
108 end
109 end
110
111 function f i nd r a and de c %Get s u b s a t e l l i t e poi nt and al s o det ermi ne i f
s a t e l l i t e can see t he t a r g e t
112 ra ( 1 : length( ti mes ) , 1: Nsat ) = 0;
113 dec ( 1 : length( ti mes ) , 1: Nsat ) = 0;
114 t het a = 0;
115 r r e l mag ( 1 : length( ti mes ) , 1: Nsat ) = 0;
116 for i = 1: length( ti mes )
117 t ( i ) = ti mes ( i ) ;
118 for i i = 1: Nsat
119 M( i i ) = M( i i ) + Mdot( i i ) dt ;
120 E( i i ) = kepl er E ( e ( i i ) ,M( i i ) ) ;
121 TA( i i ) = 2atan( tan(E( i i ) /2) sqrt ((1+e ( i i ) ) /(1e ( i i ) ) ) ) ;
122 r = h( i i ) 2/mu/(1 + e ( i i ) cos (TA( i i ) ) ) [ cos (TA( i i ) ) si n (TA( i i ) ) 0 ] ;
123
124 W( i i ) = W( i i ) + Wdot( i i ) dt ;
125 wp( i i ) = wp( i i ) + wpdot ( i i ) dt ;
126 coe0 = [ h( i i ) e ( i i ) W( i i ) i n c l ( i i ) wp( i i ) TA( i i ) ] ;
127 [ R0 V0] = s v f r om coe ( coe0 , mu) ; %Obtai n R, V f or p l o t t i n g
128 Rpl ot ( i , : , i i ) = [ R0 V0 ] ;
129
130 Rone = R3(W( i i ) ) ;
131 Rtwo = R1( i n c l ( i i ) ) ;
132 Rthree = R3(wp( i i ) ) ;
133 QxX = ( Rthree RtwoRone ) ;
134 R = QxX r ;
135
136 t het a = we( t ( i ) to ) ;
137 Q = R3( t het a ) ; %Transform wi t h r e s pe c t t o s pi n angul ar r ot at i on of
ear t h
138 r r e l = QR;
139 r r e l mag ( i , i i ) = norm( r r e l ) ;
140 [ al pha de l t a ] = r a and de c f r om r ( r r e l ) ;
141
142 ra ( i , i i ) = al pha ;
143 dec ( i , i i ) = de l t a ;
144 end
145 end
146
28 K. Loh
147 SEE( 1 : length( ti mes ) , 1: Nsat+1) = 0;
148 for i = 1: length( ti mes )
149 t i = ti mes ( i ) ;
150 for i i = 1: Nsat
151 a l t = r r e l mag ( i , i i ) Re ;
152 clambda = 1( h pri me 2 a l t 2) /(2Re(Re + a l t ) ) ; %Cos l ambda
153 IAA = KA. ( 1 clambda ) ;
154 l at P pr i me = pi /2 l atT ; %Transf ormat i on from l a t l ong t o
acces s area coor di nat es ( Sect i on 9. 1 Mi ssi on Geometry (
Wertz )
155 l atSSP pri me = pi /2 dec ( i , i i ) deg ;
156 del t a L = (ra ( i , i i ) deg longT) ;
157 lambda ( i , i i ) = acosd ( clambda ) ;
158 lambda0 ( i , i i ) = acosd (Re/(Re + a l t ) ) ;
159 rlambda ( i , i i ) = acosd ( cos ( l at P pr i me ) cos ( l atSSP pri me ) + si n (
l at P pr i me ) si n ( l atSSP pri me ) cos ( del t a L ) ) ; %Angul ar
Di st ance from t a r g e t t o s u b s a t e l l i t e poi nt , Di t t o Eqn 99
160 i f ( rlambda ( i , i i ) > 180)
161 rlambda ( i , i i ) = 360 rlambda ( i , i i ) ;
162 end
163 i f ( rlambda ( i , i i ) <= min( lambda ( i , i i ) , lambda0 ( i , i i ) ) )
164 SEE( i , i i ) = 1; %I f t a r g e t i s wi t hi n t he minimum of e i t h e r
l ambda or lambda0 , then , i t i s covered by t he i i t h
s a t e l l i t e
165 end
166 end
167 SEE( i , Nsat+1) = sign(sum(SEE( i , 1 : Nsat ) ) ) ; %As l ong as t her e i s one
s a t e l l i t e covered , t hen t he t a r g e t i s covered by t he
c o ns t e l l a t i o n
168 end
169
170 end %f i nd r a and de c
171
172 %
173 function f or m s e par at e c ur ve s
174 %
175 % Breaks t he ground t r ack up i nt o s epar at e curves which s t a r t
176 % and t ermi nat e at r i g h t as cens i ons i n t he range [ 0 , 360 deg ] .
177 %
178 t o l = 100;
179 curve no = 1;
180 n cur ves = 1;
Walker-Delta Satellite Constellation for Earth Observation 29
181 k = 0;
182 r a pr ev = ra ( 1 , i i ) ;
183 for i = 1: length( ra )
184 i f abs ( ra ( i ) r a pr ev ) > t o l
185 curve no = curve no + 1;
186 n cur ves = n cur ves + 1;
187 k = 0;
188 end
189 k = k + 1;
190 RA{ curve no }( k) = ra ( i , i i ) ;
191 Dec{ curve no }( k) = dec ( i , i i ) ;
192 r a pr ev = ra ( i , i i ) ;
193 end
194 end %f or m s epar at e cur ves
195
196 %
197 function pl ot gr ound t r ac k
198 %
199 for i = 1: n cur ves
200 plot (RA{ i } , Dec{ i } , . r , MarkerSi ze , 4)
201 end
202 plot ( longT/deg , l atT/deg , ow ) ;
203 % ax i s ( [ 0 360 90 90] )
204 text ( Pos i t i on , [ ra ( 1 , i i ) dec ( 1 , i i ) ] , St r i ng , [ o St ar t num2str( i i )
] )
205 text ( Pos i t i on , [ ra (end, i i ) dec (end, i i ) ] , St r i ng , [ o Fi ni s h num2str
( i i ) ] )
206 text ( Pos i t i on , [ longT/deg , l atT/deg ] , St r i ng , Target , FontSi ze ,
12 , Col or , whi te )
207 l i ne ( [ min( ra ( : , i i ) ) max( ra ( : , i i ) ) ] , [ 0 0] , Col or , k ) %t he equat or
208 end %pl ot g r ound t r ac k
209
210 function ground map
211 load( topo . mat , topo , topomap1 ) ;
212 contour ( 0: 359 , 89: 90 , topo , [ 0 0] , b )
213 axis equal
214 box on
215 set ( gca , XLim , [ 0 360] , YLim , [ 90 90] , . . .
216 XTick , [ 0 60 120 180 240 300 360] , . . .
217 Yti ck , [ 90 60 30 0 30 60 90] ) ;
218 hold on
219 image ( [ 0 360] , [ 90 90] , topo , CDataMapping , s c al e d ) ;
30 K. Loh
220 colormap( topomap1 ) ;
221 xlabel ( East l ongi t ude ( degr ees ) )
222 ylabel ( Lat i t ude ( degr ees ) )
223 axis equal
224 grid on
225 end %ground map
226
227 function output
228 %. . . Pl ot t he r e s u l t s :
229 % Draw t he pl anet
230 load( topo . mat , topo , topomap1 ) ;
231 [ xx , yy , zz ] = sphere ( 100) ;
232 cl a reset
233 props . Ambi entStrength = 0 . 1 ;
234 props . Di f f us e St r e ngt h = 1;
235 props . Specul ar Col or Ref l ect ance = . 5 ;
236 props . Specul arExponent = 20;
237 props . Specul ar St r engt h = 1;
238 props . FaceCol or= t ext ur e ;
239 props . EdgeCol or = none ;
240 props . FaceLi ghti ng = phong ;
241 props . Cdata = topo ;
242 s f = surface (Rexx , Reyy , Rezz , props ) ;% f acecol or , texturemap , cdata
, t opo ) ;
243 for r t = 1: 100 %Reor i e nt t opo map t o l ong i t ude eas t 0 at GMT
244 rotate ( s f , [ 0 , 0 , 1 ] , 4 5 ) ;
245 end
246 % s f ;
247 %col ormap ( l i g h t g r a y )
248 caxis ([ Re/10 Re /10] )
249 %shadi ng i nt e r p
250 % Draw and l a b e l t he X, Y and Z axes
251 l i ne ( [ 0 2Re ] , [ 0 0] , [ 0 0 ] ) ; text (2Re , 0 , 0 , X )
252 l i ne ( [ 0 0] , [ 0 2Re ] , [ 0 0 ] ) ; text ( 0 , 2Re , 0 , Y )
253 l i ne ( [ 0 0] , [ 0 0] , [ 0 2Re ] ) ; text ( 0 , 0 , 2Re , Z )
254 % Pl ot t he or bi t , draw a r a d i a l t o t he s t a r t i ng poi nt
255 % and l a b e l t he s t a r t i ng poi nt ( o) and t he f i n a l poi nt ( f )
256 hold on
257 for i i =1: Nsat
258 plot3 ( Rpl ot ( : , 1 , i i ) , Rpl ot ( : , 2 , i i ) , Rpl ot ( : , 3 , i i ) , k )
259 % l i ne ( [ 0 r0 (1) ] , [ 0 r0 (2) ] , [ 0 r0 (3) ] )
260 i f ( i i ==1)
Walker-Delta Satellite Constellation for Earth Observation 31
261 draw sat ( Rpl ot ( 1 , 1: 3 , i i ) , [ 350 350 350] , g , 1) ;
262 el se
263 draw sat ( Rpl ot ( 1 , 1: 3 , i i ) , [ 350 350 350] , r , 1) ;
264 end
265 end
266 % Se l e c t a vi ew di r e c t i on ( a vect or di r e c t e d outward from t he or i g i n )
267 view( [ 1 , 1 , . 4 ] )
268 % Speci f y some pr ope r t i e s of t he graph
269 grid on
270 axis equal
271 xlabel ( km )
272 ylabel ( km )
273 zl abel ( km )
274 end %out put
275
276 %
277 function pr i nt o r bi t a l da t a
278 %
279 coe = [ h e Wo i n c l wpo TAo ] ;
280 [ ro , vo ] = s v f r om coe ( coe , mu) ;
281 f pri ntf ( \n \n )
282 f pri ntf ( \n Angul ar momentum = %g km2/ s , h)
283 f pri ntf ( \n Ec c e nt r i c i t y = %g , e )
284 f pri ntf ( \n Semimajor axi s = %g km , a )
285 f pri ntf ( \n Per i gee r adi us = %g km , rP)
286 f pri ntf ( \n Apogee r adi us = %g km , rA)
287 f pri ntf ( \n Peri od = %g hours , T/3600)
288 f pri ntf ( \n I nc l i na t i o n = %g deg , i n c l /deg )
289 f pri ntf ( \n I n i t i a l t r ue anomaly = %g deg , TAo/deg )
290 f pri ntf ( \n I n i t i a l RA = %g deg , Wo/deg )
291 f pri ntf ( \n RA dot = %g deg/day , Wdot/deg ( t f to ) )
292 f pri ntf ( \n I n i t i a l wp = %g deg , wpo/deg )
293 f pri ntf ( \n wp dot = %g deg/ per i od , wpdot/degT)
294 f pri ntf ( \n )
295 f pri ntf ( \n r0 = [%12g , %12g , %12g ] (km) , ro ( 1) , ro ( 2) , ro ( 3) )
296 f pri ntf ( \n magnitude = %g km\n , norm( ro ) )
297 f pri ntf ( \n v0 = [%12g , %12g , %12g ] (km) , vo ( 1) , vo ( 2) , vo ( 3) )
298 f pri ntf ( \n magnitude = %g km/ s \n , norm( vo ) )
299 f pri ntf ( \n \n )
300
301 end %p r i n t o r b i t a l d a t a
302 end %gr ound t r ack
32 K. Loh
303 %
D Listing of code - orbit.m
1 function or bi t %s i mi l ar t o groundt rack .m except t hat s t a t e vect or i s
propagat ed i ns t e ad usi ng Newton s l aw and RungeKutta 5 t h order
2 cl c ; cl ose al l ; cl ear al l
3 global ra dec n cur ves RA Dec
4 hours = 3600;
5 deg = pi /180;
6 %. . . Input dat a :
7 % Earth :
8 R = 6378;
9 mu = 398600;
10 we = (2 pi + 2pi /365. 26) /(243600) ;
11 %Phys i cal model i nput
12 J2 i nc = 1; %I ncl ude J2 pe r t ur b at i on e f f e c t s ? 1 yes , 0 f or no
13 l ow t hr us t = 0; %I ncl ude l ow t hr us t t a ng e nt i a l t o o r b i t pat h? 1 yes ,
0 no
14 at t hr us t = 5e 6;
15 Atmos drag = 0; %I ncl ude at mospheri c drag ? 1 yes , 0 no
16 b a l l c o e f f = 150; % Ba l l i s t i c c o e f f i c i e n t Used f or at mospheri c drag
c a l c ul a t i o ns
17 %Si mul at i on Parameters
18 t0 = 0;
19 t f = 124 hours ;
20 dt = 20; %i n seconds f or ode5 r out i ne
21 %I n i t i a l Orbi t paramet ers
22 hP = 1000;
23 hA = 1000;
24 TAo = 0deg ;
25 Wo = 0deg ;
26 i n c l = 60 deg ;
27 wpo = 0deg ;
28 %Target l a t i t u d e and l ong i t ude
29 l atT = 36. 372 deg ; %36. 372 deg ;
30 longT = 127. 363 deg ;%127. 363 deg ;
31 %S a t e l l i t e payl oad paramet ers
32 Da = 1;
33 WL = 5e 7;
34 x r e s = 1 . 1 ;
Walker-Delta Satellite Constellation for Earth Observation 33
35 KA = 20626. 4806; % f or IAA area i n deg 2
36 %. . . End i nput dat a
37 %Cal cul at e Earth c e nt r al angl e from payl oad paramet ers
38 h pri me = Da x r e s 1e 3/(2. 440WL) ; %From Rayl ei gh o p t i c a l r e s o l ut i o n
(km)
39
40 %Numerical c ondi t i ons
41 r kr 45 i nt = 0;
42 % Obtai n R0 and V0 v e c t or s from coe dat a
43 rP = hP + R; rA = hA + R;
44 a = ( rA + rP) /2;
45 T = 2pi /sqrt (mu) a ( 3/2) ;
46 % t f = 0. 5T;
47 e = ( rA rP) /( rA + rP) ;
48 h = sqrt (mua (1 e 2) ) ;
49
50 coe = [ h e Wo i n c l wpo TAo ] ;
51 coe0 = coe ;
52 [ r0 v0 ] = s v f r om coe ( coe , mu) ;
53 %. . . Numerical i nt e g r at i on :
54 % mu = G(m1 + m2) ;
55 t het a = 0; %Earth r e v ol ut i on
56 y0 = [ r0 v0 t het a ] ;
57 i f ( r kr 4 5 i nt == 1)
58 [ t , y ] = r kf 45 ( @rates , [ t0 t f ] , y0 ) ;
59 el se
60 t = t0 : dt : t f ;
61 [ y ] = ode5 ( @rates , t , y0 ) ;
62 end
63 coe1 = coe f r om s v ( y(end, 1 : 3 ) , y(end, 4 : 6 ) ,mu) ;
64 %. . . Output t he r e s u l t s :
65 fi gure ( 1)
66 output
67 f i nd r a and de c
68 fi gure ( 2)
69 f or m s e par at e c ur ve s
70 pl ot gr ound t r ac k
71 pr i nt o r bi t a l da t a
72 fi gure ( 3)
73 hold on
74 plot ( t 60/ hours , SEE, k ) ;
75 axis ( [ t ( 1) 60/ hours t (end) 60/ hours 0 1 . 1 ] ) ;
34 K. Loh
76 xlabel ( Time ( Hours ) ) ; ylabel ( Observati on ) ;
77 fi gure ( 5)
78 plot ( t 60/ hours , rlambda , .k ) ;
79 fi gure ( 4)
80 plot ( t 60/ hours , al t , k ) ;
81
82 [ p cover mean response max response ] = t e s t t (SEE, t ) ;
83 mean al t = mean( a l t ) ;
84 f pri ntf ( \n Al t i t ude = %g km | Percentage coverage = %g per cent | Mean
r es pons e ti me = %g seconds | Max r es pons e ti me = %g seconds \n ,
mean al t , p cover 100 , mean response , max response )
85 % ax i s ( [ 0 t f /hours 0 1 . 1 ] ) ;
86 return
87
88 function dydt = r at e s ( t , f )
89 x = f ( 1) ;
90 y = f ( 2) ;
91 z = f ( 3) ;
92 vx = f ( 4) ;
93 vy = f ( 5) ;
94 vz = f ( 6) ;
95 % we = 0; %(2 pi + 2 pi /365. 26) /(243600) ;
96 r = norm( [ x y z ] ) ;
97 i f ( ( J2 i nc == 1) && ( i n c l = 0) )
98 %To c a l c u l a t e J2 o r b i t pe r t ur b at i on from Chobotov Quest i on 10. 5 J2
99 %pe r t ur b at i on i n Cart esi an coor di nat es
100 J2 = 0. 00108263;
101 f ac = (mu/ r 2) ( 3/2) J2 (R/ r ) 2;
102 px = f ac ( x/ r ) (15( z/ r ) 2) ;
103 py = f ac ( y/ r ) (15( z/ r ) 2) ;
104 pz = f ac ( z/ r ) (35( z/ r ) 2) ;
105 el se
106 P( 1 : 3 ) = 0;
107 px = 0;
108 py = 0;
109 pz = 0;
110 end
111
112 i f ( l ow t hr us t == 1)
113 T vec = [ vx vy vx ] /norm( [ vx vy vz ] ) ;
114 H vec = cross ( [ x y z ] , [ vx vy vz ] ) ;
115 W vec = H vec/norm( H vec ) ;
Walker-Delta Satellite Constellation for Earth Observation 35
116 N vec = cross ( T vec , W vec ) ;
117 Qmat = [ N vec T vec W vec ] ;
118 at = at t hr us t ;
119 AT = at [ 0 ; 1 ; 0 ] ;
120 AT = QmatAT;
121 el se
122 AT( 1 : 3 ) = 0;
123 end
124
125 ax = mux/ r 3 + px + AT( 1) ;
126 ay = muy/ r 3 + py + AT( 2) ;
127 az = muz/ r 3 + pz + AT( 3) ;
128 dydt = [ vx vy vz ax ay az we ] ;
129 end %r at e s
130
131 function f i nd r a and de c %Get s u b s a t e l l i t e poi nt and al s o det ermi ne i f
s a t e l l i t e can see t he t a r g e t
132 ti mes=t ;
133 ra = [ ] ;
134 dec = [ ] ;
135 t het a = 0;
136 for i = 1: length( ti mes )
137 t i = ti mes ( i ) ;
138 r = [ y( i , 1 ) y( i , 2 ) y( i , 3 ) ] ;
139 v = [ y( i , 4 ) y( i , 5 ) y( i , 6 ) ] ;
140 coe1 = coe f r om s v ( r , v , mu) ;
141 Rr = r ;
142 Q = R3( y( i , 7 ) ) ;
143 r r e l = QRr ;
144 r r e l mag ( i ) = norm( r r e l ) ;
145 [ al pha de l t a ] = r a and de c f r om r ( r r e l ) ;
146
147 ra = [ ra ; al pha ] ;
148 dec = [ dec ; de l t a ] ;
149 end
150 SEE( 1 : length( ti mes ) ) = 0;
151 for i i = 1: length( ti mes )
152 t i = ti mes ( i i ) ;
153 r = [ y( i i , 1 ) y( i i , 2 ) y( i i , 3 ) ] ;
154 a l t ( i i ) = norm( r r e l mag ( i i ) ) R;
155 clambda ( i i ) = 1( h pri me 2 a l t ( i i ) 2) /(2R(R + a l t ( i i ) ) ) ;
156 IAA = KA. ( 1 clambda ) ;
36 K. Loh
157 l at P pr i me = pi /2 l atT ;
158 l atSSP pri me = pi /2 dec ( i i ) deg ;
159 del t a L = ra ( i i ) deg longT ;
160 lambda ( i i ) = 90; %acosd ( cl ambda ( i i ) ) ;
161 lambda0 ( i i ) = acosd (R/(R + a l t ( i i ) ) ) ;
162 rlambda ( i i ) = acosd ( cos ( l at P pr i me ) cos ( l atSSP pri me ) + si n (
l at P pr i me ) si n ( l atSSP pri me ) cos ( del t a L ) ) ;
163 i f ( rlambda ( i i ) > 180)
164 rlambda ( i i ) = 360 rlambda ( i i ) ;
165 end
166 i f ( rlambda ( i i ) <= min( lambda ( i i ) , lambda0 ( i i ) ) )
167 SEE( i i ) = 1;
168 end
169 end
170
171 end %f i nd r a and de c
172
173 function output
174 for i = 1: length( t )
175 r ( i ) = norm( [ y( i , 1 ) y( i , 2 ) y( i , 3 ) ] ) ;
176 end
177 [ rmax imax ] = max( r ) ;
178 [ rmin i mi n ] = min( r ) ;
179 v at rmax = norm( [ y( imax , 4 ) y( imax , 5 ) y( imax , 6 ) ] ) ;
180 v at r mi n = norm( [ y( imin , 4 ) y( imin , 5 ) y( imin , 6 ) ] ) ;
181 f pri ntf ( \n\n\
n )
182 f pri ntf ( \n Earth Orbi t \n )
183 f pri ntf ( %s \n , dat e s t r (now) )
184 f pri ntf ( \n The i n i t i a l po s i t i on i s [%g , %g , %g ] (km) . , . . .
185 r0 ( 1) , r0 ( 2) ,
r0 ( 3) )
186 f pri ntf ( \n Magnitude = %g km\n , norm( r0 ) )
187 f pri ntf ( \n The i n i t i a l ve l o c i t y i s [%g , %g , %g ] (km/ s ) . , . . .
188 v0 ( 1) , v0 ( 2) ,
v0 ( 3) )
189 f pri ntf ( \n Magnitude = %g km/ s \n , norm( v0 ) )
190 f pri ntf ( \n I n i t i a l ti me = %g h. \ n Fi nal ti me = %g h. \ n , 0 , t f / hours )
191 f pri ntf ( \n The minimum a l t i t ude i s %g km at ti me = %g h . , . . .
192 rminR, t ( i mi n ) /
hours )
193 f pri ntf ( \n The speed at that poi nt i s %g km/ s . \ n , v at r mi n )
Walker-Delta Satellite Constellation for Earth Observation 37
194 f pri ntf ( \n The maximum a l t i t ude i s %g km at ti me = %g h . , . . .
195 rmaxR, t ( imax) /
hours )
196 f pri ntf ( \n The speed at that poi nt i s %g km/ s \n , v at rmax )
197 f pri ntf ( \n\n\
n )
198 %. . . Pl ot t he r e s u l t s :
199 % Draw t he pl anet
200 load( topo . mat , topo , topomap1 ) ;
201 [ xx , yy , zz ] = sphere ( 100) ;
202 cl a reset
203 props . Ambi entStrength = 0 . 1 ;
204 props . Di f f us e St r e ngt h = 1;
205 props . Specul ar Col or Ref l ect ance = . 5 ;
206 props . Specul arExponent = 20;
207 props . Specul ar St r engt h = 1;
208 props . FaceCol or= t ext ur e ;
209 props . EdgeCol or = none ;
210 props . FaceLi ghti ng = phong ;
211 props . Cdata = topo ;
212 for j j =1: length( t )
213 c l f
214 f i l ename1 =[ . /Anim/ t e s t , num2str( j j ) , . png ] ;
215 s f = surface (Rxx , Ryy , Rzz , props ) ;% f acecol or , texturemap , cdata ,
t opo ) ;
216 for r t = 1: 100 %Reor i e nt t opo map t o l ong i t ude eas t 0 at GMT
217 rotate ( s f , [ 0 , 0 , 1 ] , 4 5 + y( j j , 7 ) , [ 0 0 0 ] ) ;
218 end
219 caxis ([ R/5 R/5] )
220 % Draw and l a b e l t he X, Y and Z axes
221 l i ne ( [ 0 2R] , [ 0 0] , [ 0 0 ] ) ; text (2R, 0 , 0 , X )
222 l i ne ( [ 0 0] , [ 0 2R] , [ 0 0 ] ) ; text ( 0 , 2R, 0 , Y )
223 l i ne ( [ 0 0] , [ 0 0] , [ 0 2R] ) ; text ( 0 , 0 , 2R, Z )
224 % Pl ot t he or bi t , draw a r a d i a l t o t he s t a r t i ng poi nt
225 % and l a b e l t he s t a r t i ng poi nt ( o) and t he f i n a l poi nt ( f )
226 hold on
227 view( [ 1 , 1 , 0 . 5 ] )
228 % pl ot 3 ( y ( : , 1) , y ( : , 2) , y ( : , 3) , k )
229 draw sat ( [ y( j j , 1 ) y( j j , 2 ) y( j j , 3 ) ] , [ 3 5 0 350 350] , g , 1) ;
230 % l i ne ( [ 0 r0 (1) ] , [ 0 r0 (2) ] , [ 0 r0 (3) ] )
231 % t e x t ( y ( 1 , 1) , y ( 1 , 2) , y ( 1 , 3) , o )
232 % t e x t ( y ( end , 1) , y ( end , 2) , y ( end , 3) , f )
38 K. Loh
233 % Se l e c t a vi ew di r e c t i on ( a vect or di r e c t e d outward from t he or i g i n )
234
235 % Speci f y some pr ope r t i e s of t he graph
236 grid on
237 axis ([ 1. 5R 1. 5R 1.5R 1. 5R 1.5R 1. 5R] ) ;
238 xlabel ( km )
239 ylabel ( km )
240 zl abel ( km )
241 drawnow
242
243 print dpng f i l ename
244 FILE=i mread ( f i l ename , png ) ;
245 i mwri te ( FILE, f i l ename1 , png ) ;
246 end
247 %
248 end %out put
249 %
250
251 %
252 function f or m s e par at e c ur ve s
253 %
254 % Breaks t he ground t r ack up i nt o s epar at e curves which s t a r t
255 % and t ermi nat e at r i g h t as cens i ons i n t he range [ 0 , 360 deg ] .
256 %
257 t o l = 100;
258 curve no = 1;
259 n cur ves = 1;
260 k = 0;
261 r a pr ev = ra ( 1) ;
262 for i = 1: length( ra )
263 i f abs ( ra ( i ) r a pr ev ) > t o l
264 curve no = curve no + 1;
265 n cur ves = n cur ves + 1;
266 k = 0;
267 end
268 k = k + 1;
269 RA{ curve no }( k) = ra ( i ) ;
270 Dec{ curve no }( k) = dec ( i ) ;
271 r a pr ev = ra ( i ) ;
272 end
273 end %f or m s epar at e cur ves
274
Walker-Delta Satellite Constellation for Earth Observation 39
275 %
276 function pl ot gr ound t r ac k
277 %
278 load( topo . mat , topo , topomap1 ) ;
279 contour ( 0: 359 , 89: 90 , topo , [ 0 0] , b )
280 axis equal
281 box on
282 set ( gca , XLim , [ 0 360] , YLim , [ 90 90] , . . .
283 XTick , [ 0 60 120 180 240 300 360] , . . .
284 Yti ck , [ 90 60 30 0 30 60 90] ) ;
285 hold on
286 image ( [ 0 360] , [ 90 90] , topo , CDataMapping , s c al e d ) ;
287 colormap( topomap1 ) ;
288 xlabel ( East l ongi t ude ( degr ees ) )
289 ylabel ( Lat i t ude ( degr ees ) )
290 axis equal
291 grid on
292 for i = 1: n cur ves
293 plot (RA{ i } , Dec{ i } , r )
294 end
295 plot ( longT/deg , l atT/deg , or ) ;
296 % ax i s ( [ 0 360 90 90] )
297 text ( ra ( 1) , dec ( 1) , o St ar t )
298 text ( ra (end) , dec (end) , o Fi ni s h )
299 text ( longT/deg , l atT/deg , Target )
300 l i ne ( [ min( ra ) max( ra ) ] , [ 0 0] , Col or , k ) %t he equat or
301 end %pl ot g r ound t r ac k
302
303 %
304 function pr i nt o r bi t a l da t a
305 %
306 coe = [ h e Wo i n c l wpo TAo ] ;
307 [ ro , vo ] = s v f r om coe ( coe , mu) ;
308 f pri ntf ( \n \n )
309 f pri ntf ( \n Angul ar momentum = %g km2/ s , h)
310 f pri ntf ( \n Ec c e nt r i c i t y = %g , e )
311 f pri ntf ( \n Semimajor axi s = %g km , a )
312 f pri ntf ( \n Per i gee r adi us = %g km , rP)
313 f pri ntf ( \n Apogee r adi us = %g km , rA)
314 f pri ntf ( \n Peri od = %g hours , T/3600)
315 f pri ntf ( \n I nc l i na t i o n = %g deg , i n c l /deg )
316 f pri ntf ( \n I n i t i a l t r ue anomaly = %g deg , TAo/deg )
40 K. Loh
317 % f p r i n t f ( \n Time s i nce pe r i g e e = %g hours , t o /3600)
318 f pri ntf ( \n I n i t i a l RA = %g deg , Wo/deg )
319 % f p r i n t f ( \n RA dot = %g deg/ peri od , Wdot/deg T)
320 f pri ntf ( \n I n i t i a l wp = %g deg , wpo/deg )
321 % f p r i n t f ( \n wp dot = %g deg/ peri od , wpdot /deg T)
322 f pri ntf ( \n )
323 f pri ntf ( \n r0 = [%12g , %12g , %12g ] (km) , ro ( 1) , ro ( 2) , ro ( 3) )
324 f pri ntf ( \n magnitude = %g km\n , norm( ro ) )
325 f pri ntf ( \n v0 = [%12g , %12g , %12g ] (km) , vo ( 1) , vo ( 2) , vo ( 3) )
326 f pri ntf ( \n magnitude = %g km\n , norm( vo ) )
327 f pri ntf ( \n \n )
328
329 end %p r i n t o r b i t a l d a t a
330 end %o r b i t
331 %
E Listing of code - test t.m
1 function [ p cover mean response max response mi n response ] = t e s t t ( x , t )
2 dt = ( t (end) t ( 1) ) /( length( t ) ) ;
3 p cover = sum( x) /length( x) ;
4 t r e s pons e ( 1 : length( t ) ) = 0;
5 for i =1: length( t )
6 i f ( x( i ) > 0)
7 t r e s pons e ( i ) = 0;
8 el s e i f ( i =length( t ) )
9 i i = i ;
10 while ( ( x( i i ) <= x( i i +1) )&&(x( i i ) < 1) ) %End of Gap reached when
x ( i i ) == 1
11 t r e s pons e ( i ) = t r e s pons e ( i ) + dt ;
12 i i = i i + 1;
13 i f ( i i == length( t ) )
14 break
15 end
16 end
17 el se
18 t r e s pons e ( i ) = t r e s pons e ( i ) + dt ; %Al ready at t he end , and i t
i s a gap , so , add dt
19 end
20 end
21
Walker-Delta Satellite Constellation for Earth Observation 41
22 nt r e s pons e = nonzeros ( t r e s pons e ) ;
23 t o t a l r e s po ns e = sum( t r e s pons e ) ;
24 max response = max( nt r e s pons e ) ;
25 mi n response = min( nt r e s pons e ) ;
26 mean response = t o t a l r e s po ns e /length( t ) ;