Ch 1 Introduction Integrated Aircraft Navigation by James

Ch. 1 Introduction: Integrated Aircraft Navigation by James L. Farrell • Navigation is performed by avionics system. Navigation is the science of determining/estimating/measuring position and velocity of the flight vehicle at any instant. • 1. 1 Basic Motion patterns Rigid body dynamic involves § rotation of an object about its mass center § translation of that mass center § absolute (wrt a fixed ref. ) § relative (wrt a changing ref. ) The motion may or may not be constrained.

• Velocity and Position on a Track Fixed Ref. Point § Translation at a fixed speed Xc (t) = Xc (t₀) + Vc (t-t₀) (1. 1) § Translation at a variable speed t Xc (t) = Xc(t₀) + ∫ Vc dτ (1. 2) t₀ • Egs. (1. 1 -1. 2) can be used for translation restricted to curved line paths (e. g. a meridian). In this case, Xc represents are length, and Vc tangential speed. See the figure below. Fig. 1. 2 Translation along a circular are. Instantaneous velocity Vc is tangent to the are

§ In fig. 1. 2, if Vc (t), velocity time history, is known, only the initial Xc (t₀) is needed to determine Xc (t). For constant radius ac the angle traversed is λ = [Xc – Xc (t₀)]/ac (1. 3)

Motion in Two and Three Dimensions • Motion in Two – dimension: When the third dimension of an excursion is: a) a know profile function of the position on the surface, or b) Unimportant, then Eq. (1. 2) can be applied to two directions t Xc₁ = Xc₁(t₀) + ∫ Vc₁(τ)dτ (1. 4 a) t₀ t Xc₂ = Xc₂(t₀) + ∫ Vc₂(τ)dτ t₀ (1. 4 b) Fig. 1. 3 Translation with Two DOF

• Fig. 1. 3 can also represent the ground track (the projection of flight on a level plane) for an aircraft with instantaneous horizontal position as given by (1. 4) • Motion in 3. dire at Xc₃ represent vertical position and velocity in all three directions can be express as (with components measured along fixed axes) Rc = Xc₁ Vc = Vc₁ Xc₂ Vc₂ (1. 5) Xc₃ Vc₃ t t₀ In vector notations: Rc = Rc(t₀) + ∫ Vc dτ (1. 6)

§ Let Ac = non gravitational acceleration of the vehicle Gc = gravitational acceleration so the instantaneous velocity would τ be t₀ Vc (τ) = Vc (t₀) + ∫ (Ac+Gc)dτ τ τ and position would be t₀ t₀ Rc(t) = Rc (t₀)+Vc(t₀)(t-t₀)+∫ ∫(Ac+Gc)dτ’dτ (1. 7)

Which may be moving: § Relative Position: Rr = Rc – RT (1. 9 a) § Relative velocity: Vr = Vc – VT (1. 9 b) t § Rr = Rr (t₀) + ∫ Vr dτ (1. 8) t₀ § Let AT be the total acceleration of the tracked object; then, relative acceleration is Ar = Ac + Gc – AT Then, analogous to Eg. (1. 7) t τ Rr (t) = Rr(t₀) + Vr(t₀) (t-t₀) + ∫t₀ ∫ t₀ (Ac+Gc-AT)dτ’dτ(1. 10)

Motion Over a Spherical solid • Consider the motion of an airplane on a sphere of R = R E+ n (1. 11) As h varies, R varies. The airplane Is located at latitude λ(t), longitude Φ (t), and height h(t). The Φ and λ are measured relative to the GM and equator, which define the earth – centered – earth – fixed (ECEF) frame IE JE KE as shown.

• • IE as shown KE = axis of rotations JE = KE X IE Radius of the latitude circle: • The three orthogonal components of the velocity of the aircraft relative to the rotation earth, are VN = northern velocity VE = east velocity VV = vertical, down velocity So rate of change of latitude, longitude, and height would be λ = VN/Q Φ = VE/(R cos λ) h = -VV

otion over a spherical or spheroidal solid: Equation in Inertial frame Φ = anticlockwise, but λ = clockwise (to north) • Latitude λ = VN/Q t So λ = λ (t₀) + ∫ VN dt t₀ R (1. 12)

• Longitude: Φ(t₀) t VE____ R cos λ Φ = Φ(t₀) + ∫ VE t₀ R cosλ dt • Longitude loser its meaning when λ = 90º t • Height h(t) = h(t₀) t₀ - ∫ VV dt (1. 13) (1. 14) • For nonspherical earth, we define § means sea level § Ellipsoid of revolution about KE § Ellipsoid used to define geodetic lat-Lon § Gravity force is normal to the ellipsoid § Local vertical doesn’t quite pass three the earths mass center.

• Equation of motion in the fixed, that is, inertial frame: Applying Newton’s laws of motion: • d. Ve = Vc₁ = d²Rc = Xc₁ dt Vc₂ dt² Yc₂ Vc₃ Zc₁ = Ac + Gc (1. 15) For notations, see P. 2 • For more details about nonspherical earth, see sec. 9. 4. 3, Earth models in Grow as. M. S. , weill, L. R. , and Andrews, A. P. , Global positioning systems, Inertial Navigation, and Integration, Ed. II, wiley also see, Rogers, R. M. , Applied math In Integrated Nav systems, 3 rd Ed. , AZAA 2007.

• Eqs. Of motion (1. 15) are not much useful because Vc is defined in an inertial frame, not In the ECEE (earth – centered – earth – fixed) frame IF JE KE, or better yet, in the geographic frame IG JG KG defined locally.

Navigation coordinate frames §Sec. 1. 2. 1 Relation of ℐE with ℐZ KE = KI = Earth’s spin axis is not precessing or nutating The unit vector IE, JE rotate about 2π rad in 24 hrs, so the points on the equator translator 2πReq in 24 hrs. so their speed is § There ate three frames: Inertial frame IF JE KE ℐI ECEE frame : IEJEKE : ℐE Geographic frame at the vehicle CM or nav center: 2π X 6378 km/24 hrs = 1669 km/hr =0. 46 km/s [1 maut. Mile = 1. 852 km] = 901 naut. Miles/hr = 901 knots = 0. 25 naut. Miles/sec

§Clearly, Ie coincides with I Ionce a day. §∑ = Θ defined earlier, and IE cos ∑ sin∑ JE = -sin∑ (1. 22) §Sidereal rate of the earth about K₁ = KE axis: Ws = 360. 9854º = 0. 004178 deg/s 24 X 3600 = 7. 292115 rad/s (1. 20) The earth’s gravitational accn. Including the earth’s rotation is [Meriam and kraig p. 248] grel = g – Req W²s cos²λ = 9. 825 – 0. 03382 cos² λ m/s² KE 0 0 II cos∑ 0 JI 0 1 KI

1. 2. 2 Relation of Geographic Axes to ECEF Frame IE JE KE KE Φ IE’ JE’ KE’ J’E λ Clockwise IE”JE”KE” IG = K”E JG = J”E (1) KG = I”E §The longitude Φ is anticlockwise about the KE – axis (north pole axis) §The latitude λ is clockwise about the local east J’E – axis §Note that, at the intersection of the GM and equator : IE = local up, JE = local east, KE = local north Hence, after the rotations Φ and λ, the resulting I J K frame will still be UEN (up, east, north), but what we need is NED (north, east, down) frame. So some relieving of the axes will be required.

§ The coordinate transformation matrices are: • • • I”E = Cλ 0 Sλ CΦ SΦ 0 I”E J”E = 0 1 0 -SΦ CΦ 0 J”E K”E = -Sλ 0 c λ 0 0 1 K”E where C(. ) = cos (. ) clockwise s(. ) = sin (. ) = CΦCλ CλSΦ -SΦ CΦ Sλ I”E 0 J”E s λSΦ -s λSΦ Cλ K”E § Hence, reshuffling according to Eq. (1), we have I”E = -SλCΦ J”E K”E -SΦ CλCΦ - SλSΦ CΦ -CλSΦ Cλ I”E (1. 26) 0 J”E (1. 27) -Sλ K”E (1. 28)

§ Earth rate of rotation: WI KI = Ws KE = WS [Cλ IG – Sλ KG] § Velocity of the point on the earth’s surface directly below the aircraft, due to the sidereal rate: Ws [Cλ IG – Sλ KG] X [-R earth KG] = Ws R earth Cλ IG = easterly velocity = Ws R earth Cλ [-SΦ IE + CΦ JE] (1. 29 a) (1. 29 b)

Start unit vector in the ECI and ECEF frames: SHA = Sidereal Hour Angle S = declination For Polaris star: SHA = 2 h 31 m 49. 09 s S = 89º 15’ 50. 8” Both angles are clockwise. CE = start sightline unit vector in the ℐE CI= start sightline unit vector in ℐI Coordinate frames: XI YI ZI SHA X I’ Ѕ ZI Y I’ Z I’ YZ XЅ YЅ ZЅ

S ₀ CE = 1 0 0 = [-Ѕ]Y [-SHA]Z II JI KI = [Ѕ]Y C(SHA) – S(SHA) 0 + s(SHA) C(SHA) 0 JI 0 0 1 KI II = CS 0 + SS C (SHA)II – S(SHA)JI 0 1 0 S (SHA)II + C(SHA)JI -SS 0 CS KI CS[C(SHA)II – S (SHA) JI] + SS KI …… ……

ℐ Hence C E I = CS C(SHA) II – CS S(SHA) JI + SS KI ℐZ = CSC (SHA) -CSS(SHA) SS E ⇒ ℐE E ℐECEF C = CSS (SHA + ∑) -CSS (SHA + ∑) SS (1. 35)

Motion of the Geographic Frame FG: IG JG KG • The earth is rotating, so the frame FG is rotating with it also. The rate of earths spin is WS KZ = Ws KE = Ws [Cλ IG – Sλ KG] • The frame FG locates the airplane whosse velocity in this frame, relative to the ground, is VG = VN IG + VE JG + VV KG • This velocity changes the lat-lan of the plane thus lat : λ = VN / (RE + h) lan : Φ = VE / [(RE + h) Cλ These rates, in vector form, written together in FE and FG frames are: Φ KE – λG Λ is clockwise, so “_” Φ (Cλ IG – Sλ KG) – λ JG The origin of the frame FG translates with the velocity WSR earth Cλ JG [easterly Vel. ]

• Thus the total inertial, absolute rate of roatation of the FG frame is, denoted as WG, WG = Ws [ Cλ IG – Sλ KG] + Φ [Cλ IG – Sλ KG] – λ JG = (Ws + Φ) [Cλ IG –SλKG] – λ JG ℐ WG = (Ws + Φ)CλG -λ -(Ws + Φ) Sλ The inertial velocity of the frame FG is, denoted UG for later analysis, UG = UN IG + (VE + Ws R earth Cλ)JG + Vv KG

Jan. 18, 2012 Wed. Jan. 6, 2011 wed. 10: 15 AM 1. 2. 3 Aircraft and platform Axes : The platform is a north – slaved locally level, gimbals mechanization. The platform tilt angles relative to the geographic frame: XN = orientation error about The frame IP JP KP, rigidly attached to North, IG – axis the stabilized frame, is continuously = KG. JP driven in an attempt to maintain coincidence with the geographic XE = IG. KP X = JG. IP (1. 30)frame IG JG KG.

The plat form is a north – slaved locally level, gimbals mechanization Platform orientation Errors

The angles of the aircraft frame relative to the plat form frame IP JP KP [ ideally this is IG JG KG] Heading X KP (1. 33) IA 2 JA 2 KA 2 Pitch: Φ JA 2 IA 1 JA 1 KA 1 IA 2 Roll : Φ IA 2 JA 2 (1. 32) I A JA K A KA 2 IA 1 IA 2 = CX SY JA 2 SX CX 0 0 IP JP KA 2 1 KP IA 1 = 0 0 CΘ 0 -SΘ JA 1 0 KA 1 SΘ 0 CΘ • For a strap IA down = 1 INS, 0 the 0 gyro/accelerometer frame would be, ideally, aligned with the ℒA: IA JA KA and IP JA be the 0 mathematical CΘ SΘ nav. JP KP would

Inertial frame • Basically an inertial frame is one that does not move, or it is frozen, or it moves with a constant translational velocity. But frequently, depending on the duration or interest, this definition is mildly violated. Two inertial frame

Figures courtesy of: Ref. cohobtor, v. , orbital mechanics, AIAA Education series, 2002. Right – handed system: 1. The coordinate frame XE YE ZE (E = ecliptic) XE = vernal equinoz ZE YE spring equal day and night = to Ecliptic north pole = completes a right – handed frame 2. Earth – centered equatorial inertial frame (ECI frame) : XECI YECI ZECI XECI = XE = vernal equinox ZECI = to Geographic or celestial north pole XECI – YECI : Equatorial plane YECI completes a right – handed frame

Relationship of Equatorial Frame with Ecliptic Frame XE YE ZE (Ecliptic frame) XE E (clockwise) XECI YECI ZECI (Equatorial frame from the figure: by inspection XECI = XE YECI = YE cos. E – ZE sin E ZECI = YE sin. E+ ZE cos E ZECI = YE sin E + ZE cos E or CECI, E XECI = 1 0 0 X YECI 0 CE -SE YE ZECI 0 S E CE ZE Where CE = cos E SE = sin E

• Likewise, we can also express XE YE ZE in terms of XECI YECI ZECI: XE = XECI YE = YECI cos. E + ZECI sin E ZE= - YECI sin. E + ZECI cos. E or XE YE ZE 1 o o 0 C E SE CE, CEI 0 -SE CE XECI YECI ZECI Observe that 1) CE, ECI = C ECI, E 2) XE = CE, ECI YE ZE XECI = YECI ZECI CE, ECI CECI, E Hence : CE, ECI CECI, E 13 x 3 XE (3 x 3 identity matrix) YE E So. ZCE, ECI = C-1 ECI, E = CT ECI, E

Earth – sun unit vector in the ECI frame U sun = unit vector from the earth to the sun In the ecliptic inertial frame uℐE sun = XE = cos 2 YE = sin 2 ZE = 0

In the earth – centered inertial frame ECI: ECI ℒ u sun XECI YECI ZECI = 1 0 0 0 CE SE 0 -SE CE C 2 S 2 0 = C 2 CE S 2 SE S 2

Angular coordinates of the sun in the ECI frame The sun is located relative to the earth equator in terms of Right ascension angle œs Declination angle Ss The axed of rotation of these angles are: 1) XECI 2) Fig. 2. 3 Motion of the sun in the ECI reference frame (from Ref. 3) YECI ZECI (anticlockwise) X’ Y’ ZECI œs Z

X’ Cœs Sœs Y’ = -Sœs Cœs 0 Z’ 0 0 0 1 XECI YECI ZECI 2) X’ Y’ Z’ Y’ Xs Xs = Ys Zs œs (clockwise) Ys Zs (Ys = Y’) (Zs = ZE) CSs 0 0 + SSs X’ 1 0 Y’ CSs Z’ -SSs 0

Determination of Angular Coordinates of the Sun in the ECI frame Objective: Determine the right ascension angle œs and declination Ss in terms of known angles U and E Xs Ys Zs ℒ Solution : The Earth-Sun unit vector Us in the frame Xs, Ys, Zs is: 1 , i. e. Xs = 1 ECI We. ECI need to transform this to the frame ℒ where 0 Ys U 0 ℒsu 0 Zs 0 n is known Diversion : 1) we saw that Xs = CSs 0 Ys 0 1 Zs -SSs 0 Transpo se Then : X’ Y’ Z’ Xs = C-1 s, int Ys = CTs, int Zs Xs Ys Zs SSs 0 CSs X’ Y’ Z’ ≜ Cs, int X’ Y’ note the Z’ sequence ≜ Xs C int, s Ys Zs

Note : Cs, int = the transformation matrix from an intermediary frame to the s frame Cint, s = C-1 s, int = CT s, int These transformation matrices are orthogonal , i. e. , orthogonal + normal Orthogonal means : a raw vector dot multiplied with any other raw is zero a column vector dot multiplied with any other raw is zero Normal means: a raw (column) vector dot multiplied with it self is equal to 1. -1 T

Determination of Angular coordinates of the sun in the ECI (cont’d) Then ℒX’Y’ = Us Z’ by visual CSs 0 -SSs 1 0 0 0 Fig. 2. 3, by -sun X’Y’Z’ ℒ CSs = SSs X’Y’Z 0 ℒ CSs ’ 0 ECI SSs ℒ This can be verified inspection of projecting the Earth vector on the axes ZECI 2) Next transform Us to Us From X’ = Cœs Sœs 0 XECI = Cint, ECI Y’ -Sœs Cœs 0 YECI Z’ 0 0 1 Z ECI We get : XECI = CECI, int X’ YECI Y’ ZECI Z’ = Cœs -Sœs 0 sœs Cœs 0 0 0 1 X’ and Z’ = : XECI YECI ZECI X’ Y’ Z’

ECI Then ℒ = XECI Us. C 2 YECI CE S 2 (1) ZECI SE S 2 = Cœs - Sœs 0 0 Cœs 0 ℒ 1 X’Y’Z’ CSs ℒ ECI = Cœs CSs 0 Sœs CSs SSs ECI ℒ = From the three component equations of (1), determine Ss and œs : (a) Sin Ss = sin E sin 2 where – π/2 < S₁<π/2 ; E = 23. 440 ; v = 0 2π Where CSs> 0 because - π/2 <Ss < π/2 So Ss = sin-1 (sin. E sin 2) (b) Sœs CSs = CE S 2 Cœs CSs = C 2 tan œs = CE tan 2 ⇒ œs = atan 2 (CE sin 2), œs and v are in the same

Fig. 2. 4 Solar declination and right ascension VS date (from Ref. 3)

Rotation of the Greenwich meridian relative to J 2000 Frame q. The Earth takes 365. 25 days to go around the sun. so the mean motion of the Earth around the sun is Fig. 2. 8 b sidereal day and mean solar day [Ref. Moral and Bouquet, satellite communication systems, Wiley, 2006] q. The right ascension angle of the Greenwich meridian is: J 2000 Θ = 280. 46060 + 360. 985 60 where d = the time in days = Greenwich Hour Angle = the time in seconds/86, 400 s

q In one orbital period T, Θ changes by 360. 9856 o/ day x. T in day = 360. 9856 X T in minutes 24 x 60 minutes = 0. 2507 0/x T in minutes = 0. 2070/minute x T in minutes q Geographic longitude λ = œ – Θ (2. 8) so, because of – Θ, λ shifts westward by – 0. 25070/minute. T (orbital period in minutes) Homework: Show the ground track of an INSAT transfer orbit and an IRS polar orbit for three orbital periods.

Greenwich Hour Angle at Jan. 1. 2000 Noon 12: 00 h (J 2000) q The source of the angle Θ₀= 280. 4606º at 12 th noon on Jan 1. 2000 Review the heliocentric ecliptic coordinate frame Below, the angle Θ₀ is shown to be 279. 05º, using preliminary calculations, up to Jan. 1, 2000, 12: 00 noon. Autumnal Equinox, Sept. 23 rd 1999, 0 h (midnight)

Greenwich Hour Angle at Jan. 1. 2000, 12: 00 hr. (Noon) (cont’d) • To identify the source of remaining difference of ~1. 40, we examine relative motion of the sun more closely. • The Earth’s orbit around the sun is not circular. Indeed, it is elliptical with e = 0. 01673 As shown here. For or elliptical orbit, The true anomaly rate is V = √ 1 + e cosv r Where v is measured from the perihelion. Because r = p / (1 + ec 2) p = a (1 -e 2) n= =mean rate we get v = n ( 1 + ec 2)2 Figure 2. 5 Orbit of the earth round the sun.

Greenwich Hour Angle at Jan. 1, 2000, 12. 00 hr (noon) (cont’d) • To calculate Θ₀, earlier we used 0. 9856º/day which is the mean rate, n. for more precise calculations, one should use v and integrate it from Sept. 23: 99 to Jan. 1, 2000. 0 h. the sun’s motion in the