Modeling And Simulation Of Prolate Dual-Spin Satellite Dynamics In An Inclined Elliptical Orbit: Case Study Of Palapa B2R Satellite J. Muliadi* , S.D. Jenie† and A. Budiyono‡ Abstract— In response to the interest to re-use Palapa B2R satellite nearing its End of Life (EOL) time, an idea to incline the satellite orbit in order to cover a new region has emerged in the recent years. As a prolate dual-spin vehicle, Palapa B2R has to be stabilized against its internal energy dissipation effect. This work is focused on analyzing the dynamics of the reusable satellite in its inclined orbit. The study discusses in particular the stability of the prolate dual-spin satellite under the effect of perturbed field of gravitation due to the inclination of its elliptical orbit. Palapa B2R physical data was substituted into the dual-spin’s equation of motion. The coefficient of zonal harmonics J2 was induced into the gravity-gradient moment term that affects the satellite attitude. The satellite’s motion and attitude were then simulated in the perturbed gravitational field by J2, with the variation of orbit’s eccentricity and inclination. The analysis of the satellite dynamics and its stability was conducted for designing a control system for the vehicle in its new inclined orbit. Keywords—dual spin, prolate dual-spin, satellite dynamics, simulation Abstrak— Sebagai reaksi atas adanya minat dari beberapa pihak untuk mengoperasikan kembali satelit Palapa B2R yang mendekati masa akhir operasinya (End of Life), suatu ide untuk menginklinasi orbit satelit telah mengemuka pada beberapa tahun terakhir. Sebagai sebuah wahana dual-spin prolate, Palapa B2R harus distabilkan terhadap efek disipasi energi internal. Paper ini berkonsentrasi pada analisis dinamik dari satelit pada orbit barunya yang terinklinasi. Studi yang dilakukan mendiskusikan secara khusus kestabilan dari satelit dual-spin prolate dalam pengaruh efek dari medan gravitasi terganggu akibat inklinasi dari orbit eliptiknya. Data fisik Palapa B2R disubstitusikan ke dalam persamaan gerak dual-spin. Koefisien harmonik zonal J2 diinduksikan ke persamaan momen gradient gravitasi yang mempengaruhi sikap satelit. Gerak dan sikap satelit kemudian disimulasikan dengan variasi eksentrisitas dan inklinasi orbitnya. Analisis dinamika dan kestabilan satelit dilakukan untuk keperluan perancangan sistem kendali wahana pada orbit barunya yang terinklinasi. Kata kunci—dual-spin, dual-spin prolate, dinamika satelit, simulasi     * Department of Aeronautics and Astronautics, ITB, mie_pn@yahoo.com † Department of Aeronautics and Astronautics, ITB, sdjenie@ae.itb.ac.id ‡ Formerly with Department of Aeronautics and Astronautics, ITB, presently with Department of Aerospace-IT Fusion Engineering, Konkuk University, Korea, the author to whom all correspondence to be addressed budiyono@alum.mit.edu   1. INTRODUCTION A semi-rigid body is stable only when spinning about its major axis. In a related study, Bracewell and Garriott (Ref. [8] pp. 62-64) concluded that the four turnstile wire antennae of Explorer I were dissipating energy; thus, causing a transfer of body spin axis from the minimum inertia (prolate) to a transverse axis of maximum inertia (oblate). To meet this stability criterion, most of early dual-spin vehicles were designed in an oblate configuration. As a case of study, this work used Palapa B2R physical data to analyze the dynamics of the vehicle. Palapa B2R is a communication satellite of Indonesia. In its orbit, it was operated by PT Telkom (Indonesian State Telecommunication Company). Near the satellite’s End of Life (EOL) time, several Africans and Polynesians countries have shown interest to buy and re-use Palapa B2R. Because of those countries’ location in the southern latitudes, an idea emerged to incline the satellite’s orbit. The current paper elaborates the analysis of the vehicle dynamics in its inclined orbit. 2. REFERENCE COORDINATE SYSTEM 2.1. Body Reference Coordinate System (Body Axes) B X B Y B Z p x = ω q y = ω r z = ω Sattelite's Symmetrical Plane Fig. 1 Platform Axis Components The definition of Platform and Body axes is well-defined in the literature. Fig. 1 illustrated those axes with their origin at the satellite’s c.g. while Fig. 2 showed the axes in the space. In this paper, the Platform Axis Components will be identified as Body Reference Coordinate System or Body axis. Fig. 2 Body Reference Coordinate System 2.2. Stability Reference Coordinate System (Stability Axes) Fig. 3 Stability Reference Coordinate System Stability Reference Coordinate System (Stability Axes) was defined as a set of local horizon axes for the satellite. It is a target axes for the satellite’s Body Axes to point its antennae to the Earth.. All these axes are presented in Fig. 3. 2.3. Inertial Reference Coordinate System (Inertial Axes) Fig. 4 Inertial Reference Coordinate System Inertial Reference Coordinate System (Inertial Axes) is defined as a geocentric non-rotating equatorial reference frame with ZI axis which coincides with the rotation axis of the Earth and points to the North Pole; the XI axis lies in equatorial plane and points towards the vernal equinox. The YI axis completes a right- handed Cartesian frame of reference. In this Inertial Axis, Newton’s laws of motion are valid for the satellite’s translation and rotation. 3. Euler Angles (Orientation Angles) 3.1. Orientation of Body Axes in Inertial Axes Fig. 5 Orientation of Body Axes in Inertial Axes In order to describe the attitude of the satellite with respect to the Inertial Axes (Fig. 5), the Euler angles are used. Fig. 6 Euler Angles of Body Axis in Inertial Axis The yaw angle ψ, pitch angle θ and roll angle φ, respectively defines the angle of rotation in Z-, Y- and X- axis of the Body frame with respect to its nominal condition in Inertial Axes. These angles are shown in Fig. 6. 3.2. Orientation of Body Axes in Stability Axes Fig. 7 Deviation of Body Axes from Stability Axes Euler Angles of Body Axes with respect to Stability Axes are defined to describe the attitude perturbation from its local horizon (its stationary or nominal condition). If the satellite deviates to large, the antennae will point away from the Earth. The Euler Angles are: Fig. 8 Euler Angles of Body Axis with Respect to Stability Axis Yaw, Pitch and Roll Angle Deviation ( ψS, θS, and φS), which respectively denotes the angle of perturbation because of the rotation in Z-,Y-, and X-axis of the Body frame with respect to Stability Axes. These angles are shown in Fig. 8. 4. GRAVITY GRADIENT MOMENT Agrawal derived an expression for Gravity Gradient Moment in Axially Symmetric Spacecraft at Equatorial Circular Orbit (Ref [1] pp. 131-133) with this equation,       ⋅ − ⋅ − ⋅ = 0 ) ( ) ( µ 3 S B X B Z S B Y B Z 3 0 E G θ φ I I I I R M Eq. 1 Fig. 9 Concept of Gravity Gradient Moment It this paper the spacecraft will be treated as a non-axially symmetric vehicle. In addition, the satellite will be operated in inclined elliptical orbit, which means non-equatorial and non-circular orbit. So, by following and combining Kaplan’s (Ref. [8] pp. 199-204) and Agrawal’s techniques in deriving the equation of Gravity Gradient Moment for Spacecraft, the authors derived the equation of Gravity Gradient Moment for Prolate Dual-Spin Satellite in its inclined elliptical orbit. The gravitational force (dFG) corresponding to a differential element of mass, dm, shown in Fig. 9, is m d d G ⋅ = g F Eq. 2 where g denotes gravity vector at dm. The unsymmetrical mass distribution of the Earth induced a zonal harmonic coefficient (Jk, Ref [8] pp. 273-282) that perturbs the homogenous of the Earth gravity’s field. In this work, the oblateness of the Earth induced the zonal harmonic coefficient that is limited to second order, J2. The gravity vector equation at a point in space is R g ⋅               − ⋅ ⋅     ⋅ ⋅ − ⋅ − = 1 5 Re J 2 3 1 µ 2 2 E Z 2 2 3 E R R R R Eq. 3 With µE = Earth gravitational parameter; R and R = the distance from the center of the Earth in scalar and vector notation; RZE = the height measured perpendicular from the Earth equatorial plane; Re = the radius of the Earth equator. Continuing to the moment equation, ∫ ∫ ⋅ × = × = B B m d d G G g r F r M Eq. 4 while the position of differential element of mass, dm, in Fig. 9 is r R R + = m Eq. 5 Therefore, substituting Eq. 5 to replace R in Eq. 3, then inserting the result to Eq. 4, tranforms Eq. 4 into: ( ) m R B d 1 5 Re J 2 3 1 µ - 2 2 E Z 2 2 2 3 E G ⋅    + ⋅             − + ⋅        ⋅ + ⋅ ⋅ − ⋅ + × = ∫ r R r R r R r R r M Eq. 6 Fig. 10 Deviation from Steady Attitude The use of Agrawal’s technique to work with Euler Angles between Body Axis and Stability Axis see also Fig. 8), yields (eq. 3.101, Ref. [1] p. 132), ) cos (cos ) cos (sin sin S S Z S S Y S X θ φ θ φ θ ⋅ ⋅ − = ⋅ ⋅ − = ⋅ = R R R R R R Eq. 7 Binomial Expansion is applied into Eq. 6 and then Eq. 7 was inserted into the result to give Gravity Gradient Moment Equation. By linearizing the equations, the Linearized Gravity Moment Equations read, S Z Z G S Y Y G S X X G θ θ φ ⋅ = ⋅ = ⋅ = G M G M G M Eq. 8 Where coefficients of S φ and S θ are: B YZ Z B X B Z Y B Y B Z X ) ( ) ( I g G I I g G I I g G • = − • = − • = µ µ µ Eq. 9 And coefficient gµ is, Re J µ 2 15 Re J µ 2 105 µ 3 2 2 5 E 2 I Z 2 2 7 E 3 E        ⋅ ⋅ ⋅ −        ⋅ ⋅ ⋅ ⋅ + ⋅ = R R R R gµ Eq. 10 The RZE variable is the component of satellite’s position at Z-Axis in Earth. The factor of inclination was induced in RZI variable, because the height of the satellite will be varying in the inclined orbit. The factor of eccentricity was induced in R variable. For e>0, the value of R is varying along the orbit. 5. ORBITAL MOTIONS 5.1. Parameters of Keplerian Orbit Astronomy defines 6 quantities to describe the orbit and position of heavenly body, namely a, e, i, Ω, ω, and τ. The definition of those parameters can be found in many textbook in orbital mechanics. Fig. 12 describes the geometry of the orbital parameters. Fig. 11 Orbital Plane Orientation with respect to Inertial Axes 5.2. Orbital Motion Equations Simulation of the orbital motion uses scalar differentials equation in XI-axis, YI- axis and ZI-axis as follows: I Z 3 E I Z I Y 3 E I Y I X 3 E I X µ µ µ R R R R R R R R R ⋅ = ⋅ = ⋅ =    Eq. 11 The equation was derived by relating Newton’s Second Law of Motion and Newton’s Law for Gravitation. To sketch the satellite’s orbit, Eq. 11 were integrated two times with 6 initial values, initial velocity I X R , I Y R , I Z R , and 3 initial position I X R , I Y R , and I Z R , when t=0. From Jenie (Ref. [6]), initial velocity can be expressed by Eq. 12 and initial position can be expressed by Eq. 13 with transformation matrix Orb I C , and its elements as follows:           = 33 32 31 23 22 21 13 12 11 Orb I C C C C C C C C C C ) sin (sin ) sin cos cos cos sin ( ) sin cos sin cos (cos 13 12 11 i C i C i C Ω = Ω − Ω − = Ω − Ω = ω ω ω ω ) cos sin ( ) cos cos cos sin sin ( ) cos cos sin sin (cos 23 22 21 Ω − = Ω + Ω − = Ω + Ω = i C i C i C ω ω ω ω ) (cos ) sin (cos ) sin (sin 33 32 31 i C i C i C = = = ω ω When simulating the satellite’s orbit, the orbital parameters are 0 0 0 ~ = θ , Ω0 = 0o, ω0 = 0o a0 = 8078.14 km for circular orbit R = 1700 km above the sea level; for elliptic orbit with e=0.2, perigee, Rperigee = 85 km above the sea level. Fig. 102 Velocities in Orbit In the state space model for attitude dynamics, the satellite transversal velocity that perpendicular to its radius to the Earth, Vθ, will be needed. A. E. Roy (Ref. [9] p. 81) relates ϕ~ sin in Fig. 10 with eccentricity and semimajor axis as follows: 2 1 2 2 ) 2 ( ) 1( ~ sin     − ⋅ ⋅ − ⋅ = R a R e a ϕ Eq. 14 With Eq. 14, the component of velocity in theta direction can be stated as ) 2 ( ) 1( 2 2 θ R a R e a V V − ⋅ ⋅ − ⋅ ⋅ = Eq. 15 6. DYNAMICS OF PROLATE DUAL SPIN SATELLITE Palapa B2R is a prolate configuration Communication Satellite (HS 376). In order to stabilize its attitude and pointing direction, Palapa B2R uses its rotor spinning. Control moments were produced by the angular acceleration and deceleration of the rotor’s spin. Because of the rotor’s spinning motion and the configuration of satellites inertia, the satellite’s motion in the yaw mode is coupled with its roll mode. In addition, by the imbalance from the satellite’s antennae reflector configuration (IYZ), the satellite’s motion in the pitch mode is coupled with its yaw mode. 6.1. Dynamic Equations of Motion Bryson has shown the equation of spacecraft motion in De-Spin Active Nutation Damping in Ref. [2 pp. 62-68]. In this work, the author will insert gravity gradient moment as an external moment that perturbs the satellite’s attitude. Using Bryson’s technique in deriving the De- Spin Active Nutation Damping equations, the author adds several modifications. The results are: The Dynamics Equation at X-axis, X 0 R S T X ) ( ) ( M r I p I I = ⋅ Ω ⋅ − ⋅ +  Eq. 16 The Dynamics Equation at Y-axis, Y . . S YZ . . V S Y . . V YZ Y 1 M e R N r I I c R K N q I I c R K N r I q I c d c d c d = ⋅ + ⋅      + ⋅ + ⋅      + ⋅      + ⋅ + ⋅ + ⋅ δ   Eq. 17 The Dynamics Equation at Z-axis, Z 0 R S T Z YZ ) ( ) ( M p I r I I q I = ⋅ Ω ⋅ + ⋅ + + ⋅   Eq. 18 6.2. Dynamic Equations with Gravity Gradient Moments Substituting MX in Eq. 16 with MGX in Eq. 8 yields the Dynamics Equation at X- axis in Inclined Elliptical Orbit, 0 ) ( ) ( S X 0 R S T X = ⋅ − ⋅ Ω ⋅ − ⋅ + φ G r I p I I  Eq. 19 Substituting MY in Eq. 17 with MGY in Eq. 8 yields the Dynamics Equation with Gravity Gradient Moments at Y-axis in Inclined Elliptical Orbit, 0 1 . . S Y S YZ . . V S Y . . V YZ Y = ⋅ + ⋅ − ⋅      + ⋅ + ⋅      + ⋅      + ⋅ + ⋅ + ⋅ e R N G r I I c R K N q I I c R K N r I q I c d c d c d δ θ   Eq. 20 Finally, substituting MZ in Eq. 18 with MGZ in Eq. 8 yields the Dynamics Equation at Z-axis in Inclined Elliptical Orbit, 0 ) ( ) ( S X 0 R S T X = ⋅ − ⋅ Ω ⋅ − ⋅ + φ G r I p I I  Eq. 21 6.3. Kinematics Equations of Motion 6.3.1. Kinematics Equations in Inertial Axes Fig. 113 Rotation of Local Horizon Axes (Stability Axes), n Bryson has shown that the equation of kinematics in relationship between Euler Angle Rates and Angular Velocity (Ref. [2] pp. 8-9) can be approximated as φ ψ θ ψ φ ⋅ − ≅ + ≅ ⋅ + ≅ n r n q n p    Eq. 22 where n = orbital angular velocity. These kinematics equations are derived in Inertial Axis by Newton’s first and second law of motion. 6.3.2. Kinematics Equations in Inertial Axes In the beginning of simulation, the satellite’s antennae still points to the Earth for a moment. By the inertia (the Newton first law of motion) the satellite’s attitude will drift at rotor spin axis with drift rate equals to orbital angular rate, θ~ , or equals to initial value of n, n0. Therefore, to measure the deviation in Stability Axes, the kinematics equation needs to be reduced by the effect of Local Horizon Axes (Stability Axes) initial angular drift with respect to Inertial Axes, n0. Following Byson’s techniques in deriving the kinematics equation, the Kinematics Equation of Dual Spin Satellite in Stability Axis can be written as follows S S S S S φ δ ψ δ θ ψ δ φ ⋅ − ≅ + ≅ ⋅ + ≅ n r n q n p    Eq. 23 where ) ( 0 n n n − = δ . 6.4. State Space Model of Prolate Dual-Spin Combining The Dynamics Equation in Inclined Elliptical Orbit and Kinematics Equation of Dual Spin Satellite in Stability Axis yields State Space Model for Dual-Spin Satellite in Stability Axes as follows, [ ] [ ]       • +                   • =                   n e r q p r q p δ δ ψ θ φ ψ θ φ B A S S S S S S       Eq. 24 where the [A] and [B] matrices are: [ ]                     − = 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 A 0 A A A 0 A 0 A A A 0 0 A A 0 0 35 33 32 31 25 23 22 21 14 13 n n δ δ A [ ]                     = 0 0 1 0 0 0 0 B 0 B 0 0 31 21 B By defining ∆I as follows, ) ( 2 YZ T Y Z Y I I I I I I − + = ∆ the [A] matrix elements in 1st line are: ) ( ) ( 1 A 0 R S T X 13 Ω ⋅ ⋅ + − = I I I X T X 14 ) ( 1 A G I I ⋅ + − = the [A] matrix elements in 2nd line are: ) ( A 0 R S YZ 21 Ω ⋅ ⋅ ∆ − = I I I       + ⋅      + ⋅ ⋅ ∆ + = 1 A S Y . . V T Z 22 I I c R K N I I c d I      ⋅      + ⋅ ⋅ ∆ + = S YZ . . V T Z 23 A I I c R K N I I c d I Z YZ Y T Z 25 ) ( A G I G I I I I ⋅ ∆ + ⋅ ∆ + − = the [A] matrix elements in 3rd line are: ) ( A 0 R S Y 31 Ω ⋅ ⋅ ∆ = I I I       + ⋅      + ⋅ ⋅ ∆ − = 1 A S Y . . V YZ 32 I I c R K N I c d I      ⋅      + ⋅ ⋅ ∆ − = S . . V 2 YZ 33 1 A I c R K N I c d I Z Y Y YZ 35 A G I G I I I ⋅ ∆ − − ⋅ ∆ = the [B] matrix elements in 1st column are:      ⋅ ∆ + = . . T Z 21 B c d I R N I I      ⋅ ∆ − = . . YZ 31 B c d I R N I and R V n θ θ − = − = ~ R V R V n n n ) ( ) ( ~ ~ ) ( θ 0 0 θ 0 0 − = − = − = θ θ δ   7. RESULTS OF OPEN LOOP SIMULATION AND INTERPRETATIONS The values of [A]’s and [B]’s elements are shown in Eq. 25. The value of A14, A25, A35 and δn will be time-varying for elliptical or inclined orbit. However, the elements of [A] are constant value only for circular orbit at equatorial plane.                     − × × × − × − = − − − − 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 A 0 10 1.1912 - 10 3636 .3 4.0326 - 0 A 0 10 4402 .3 10 7138 .9 49773 .0 0 0 A 3.7113 0 0 35 6 5 25 5 4 14 n n δ δ A                     × × − = − − 0 0 1 0 0 0 0 10 7735 .1 0 10 1218 .5 0 0 5 4 B Eq. 25 7.1. Simulation in Longitudinal Mode 7.1.1. Effect of Eccentricity in Inclined Orbit (i = 30°) Plot of θS because of impulsive input δeref Fig. 124Plot of θS; input δeref ; i=30deg Plot of θS due to elliptic orbital drift input δn Fig. 135 Plot of θS and δn; e=0 & e=0.2; 2 Orbital Periode; i=30deg After impulsive disturbance were applied, Gravity Gradient Moment induced subsidence mode in circular orbit. While the increment of eccentricity, drove that subsidence mode into long periodic oscillation mode, which had the same period with orbital period. This long periodic oscillation mode oscillates from – 20o to +90o as superposition of many oscillation mode. 7.1.2. Effect of Inclination in Elliptic Orbit (e = 0.2) Plot of θS due to impulsive input δeref Fig. 146 Plot of θS i=0deg, 30deg dan 60deg; e=0.2 Plot of θS due to elliptic orbital drift input δn Fig. 157 Plot of θS; i=0deg, 30deg dan 60deg;1Orbit and 2 Orbits; e=0.2 Inclination increment effect in longitudinal motion mode was negligible for Dual-Spin Satellite. For inclination at 0°, 30°, and 60°, the graphic curves for θS were almost aligned. 7.2. Simulation in Lateral Mode 7.2.1. Effect of Eccentricity in Inclined Orbit (i = 30°) Plot of φS due to impulsive input δeref Fig. 18 Plot of φS; e=0 and e=0.2; i=30deg Plot of φS due to elliptic orbital drift input δn Fig. 16 Plot of φS and δn; e=0 & e=0.2; 4 Orbits; i=30deg After impulsive disturbance were applied, Gravity Gradient Moment induced roll-librations mode in for φS, which oscillates from -1×10-3 to +1×10-3 deg. In circular orbit, this roll-librations mode in for φS was damped after ±400 sec. While in elliptic orbit, the eccentricity induced long periodic oscillation mode, which oscillates from –4×10-3 to +4×10-3 deg. 7.2.2. Effect of Inclination in Elliptic Orbit (e = 0.2) Plot of φS due to impulsive input δeref Fig. 170 Plot of φS; i=0deg, 30deg dan 60deg; e=0.2 Plot of φS due to elliptic orbital drift input δn Fig. 18 Plot of φS and δn; i=0deg, 30deg and 60deg; 2 Orbits e=0.2 Inclination increment effect in longitudinal motion mode was negligible for Dual-Spin Satellite. For inclination at 0°, 30°, and 60°, the graphic curves for φS were almost aligned. 7.3. Simulation in Directional Mode 7.3.1. Effect of Eccentricity in Inclined Orbit (i = 30°) Plot of ψS due to impulsive input δeref Fig. 19 Plot of ψS; e=0 and e=0.2; i=30deg Plot of ψS due to elliptic orbital drift input δn Fig. 20 Plot of ψS and δn; e=0 & e=0.2; 4 Orbits i=30deg After impulsive disturbance were applied, Gravity Gradient Moment induced yaw-librations mode in for φS. In circular orbit, this yaw -librations mode in for ψS was damped after ±400 sec. In the elliptic orbit, the eccentricity induces long periodic oscillation mode, which oscillates from –5.5×10-3 to +1.5×10-3 (deg). 7.3.2. Effect of Inclination in Elliptic Orbit (e = 0.2) Plot of ψS due to impulsive input δeref Fig. 21 Plot of ψS; i=0deg, 30deg and 60deg; e=0.2 Plot of ψS due to elliptic orbital drift input δn Fig. 22 Plot of ψS; i=0deg, 30deg and 60deg;1 Orbit and 2 Orbits; e=0.2 Inclination increment effect in longitudinal motion mode was negligible for Dual-Spin Satellite. For inclination at 0°, 30°, and 60°, the graphic curves for ψS were almost aligned. 8. CONCLUDING REMARKS In Open-Loop simulation, the librations and subsidence mode are present in longitudinal, lateral and directional motion. However, in lateral and directional motion, the Gravity Gradient moment induced a divergence mode if the simulation were run for more than 4 orbital periods (1 Orbital Periods is 7225 sec.). The effect of Gravity Gradient Moments is destabilizing the lateral and directional motion for the dual-spin satellite. In reverse, the effect of Gravity Gradient Moments is stabilizing longitudinal motion. The effect of increasing the orbital eccentricity, e, is the presence of the long period oscillation mode in the longitudinal and lateral directional motion. The effect of increasing the inclination of Orbital Plane, i, can be neglected in the longitudinal motion. In lateral directional motion, increasing the inclination will reduce steady state error of φS and ψS. REFERENCES [1] B. N. Agrawal, Design of Geosynchronous Spacecraft, Prentice-Hall Inc., Englewood Cliffs (New Jersey), 1986. [2] A. E. Bryson, Control of Spacecraft and Aircraft, Princeton University Press, Princeton (New Jersey), 1994. [3] Cornelisse, J. W. with H. F. R. Schöyer and K.F. Wakker, Rocket Propulsion and Spaceflight Dynamics, Pitman Publishing Ltd., London, 1979. [4] Hughes Aircraft Company, Palapa- B System Summary, Hughes Aircraft Company, 1982. [5] Hughes Aircraft Company, Palapa- B2R System Summary HS-376E- 22868, Hughes Aircraft Company, 1990. [6] Jenie, Said D., Astrodynamics 1 (Lecture Notes PN-350), In Indonesian. Department of Aeronautics and Astronautics, Bandung Institute of Technology, Bandung, 2002. [7] Jenie, Said D., Flight Control (Lecture Notes PN-4041), In Indonesian. Department of Aeronautics and Astronautics, Bandung Institute of Technology, Bandung, 2003. [8] M. H. Kaplan, Modern Spacecraft Dynamics & Control, John Wiley & Sons, New York, 1976. [9] Roy, Archie E., Orbital Motion, Adam Hilger Ltd., Bristol (England), 1982.