Cbt),\t»- MUb ^^l A UNITED STATES DEPARTMENT OF COMMERCE PUBLICATION **» TO 'ceL VVr '^r.s o« ' NOAA TR NOS 42 NOAA Technical Report NOS 42 U.S. DEPARTMENT OF COMMERCE National Oceanic and Atmospheric Administration National Ocean Survey Computational Procedures for the Determination of a Simple Layer Model of the Geopotential From Doppler Observations BERTOLD U. WITTE ROCKVILLE. MD. April 1971 NOAA TECHNICAL REPORTS National Ocean Survey Series The National Ocean Survey (NOS) provides charts and related information for the safe navigation of marine and air commerce. The Survey also furnishes other earth science data—from geodetic, hydrographic, oceanographic, geomagnetic, seismologic, gravimetric, and astronomic surveys , observations, investigations, and measurements — to protect life and property and to meet the needs of engineering, scientific, defense, commercial, and industrial interests. Because many of these reports deal with new practices and techniques, the views expressed are those of the authors, and do not necessarily represent final Survey policy. NOS series of NOAA Technical Reports is a continuation of, and retains the consecutive numbering sequence of, the former series, ESSA Technical Report Coast and Geodetic Survey (CSGS) , and the earlier series, C§GS Technical Bulletin. Reports in the series are available through the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402. Those publications marked by an asterisk are out of print. Coast and Geodetic Survey Technical Bulletins *C§GS 1. Aerotriangulation Adjustment of Instrument Data by Computational Methods. William D. Harris, January 1958 (Superseded by No. 23.) *CSGS 2. Tellurometer Traverse Surveys. Lt. Hal P. Demuth, March 1958. *C£GS 3. Recent Increases in Coastal Water Temperature and Sea Level—California to Alaska. H. B. Stewart, Jr., B. D. Zetler, and C. B. Taylor, May 1958. *C5GS 4. Radio Telemetry Applied to Survey Problems. Richard R. Ross, February 1959. *CSGS 5. Raydist on Georges Bank. Capt. Gilbert R. Fish, April 1959. CSGS 6. The Tsunami of March 9, 1957, as Recorded at Tide Stations. Garrett G. Salsman, July 1959. Price $0.15. *C$GS 7. Pantograph Adjustment. G. C. Tewinkel, July 1959. *C$GS 8. Mathematical Basis of Analytic Aerotriangulation. G. C. Tewinkel, August 1959 (Superseded by No. 21.) *C$GS 9. Gravity Measurement Operations in the Field. Lt. Comdr. Hal P. Demuth, September 1959. *C5GS 10. Vertical Adjustment of Instrument Aerotriangulation by Computational Methods. William B. Harris, September 1959 (Superseded by No. 23.) *C5GS 11. Use of Near-Earth Satellite Orbits for Geodetic Information. Paul D. Thomas, January 1960. *C$GS 12. Use of Artificial Satellites for Navigation and Oceanographic Surveys. Paul D. Thomas, July 1960. C§GS 13. A Singular Geodetic Survey. Lansing G. Simmons, September 1960. Price $0.15. *C£GS 14. Film Distortion Compensation for Photogrammetric Use. G. C. Tewinkel, September 1960. *C$GS 15. Transformation of Rectangular Space Coordinates. Erwin Schmid, December 1960. *C§GS 16. Erosion and Sedimentation— Eastern Chesapeake Bay at the Choptank River. G. F. Jordan, January 1961. CSGS 17. On the Time Interval Between Two Consecutive Earthquakes. Tokuji Utsu, February 1962. Price $0.10. (Continued on inside back cover) jjtl * C v U.S. DEPARTMENT OF COMMERCE Maurice H. Stans, Secretary NATIONAL OCEANIC AND ATMOSPHERIC ADMINISTRATION Robert M. White, Administrator NATIONAL OCEAN SURVEY Don A. Jones, Director T( ecnnicai Ke Computational Procedures for Determination of a Simple Layer Model of the Geopotential From Doppler Observations BERTOLD U. WITTE a o U o c- a .1; a c/S 3 ROCKVILLE, MD. APRIL 1971 UDC 528.27:528.375:621.3.093.2:629.783 528 Geodesy .27 Gravity measurement .375 Geopotential height measurement 621.3 Electrical engineering .093 Radio-frequency phenomena .2 Doppler effect 629.78 Astronautics - Space vehicles .783 Artificial earth satellites For sale by the Superintendent of Documents, U.S. Government Printing Office Washington, D.C. 20402 - Price 65 cents 11 CONTENTS Page Abstract 1 1. Introduction 2 2. Observational Equations 2 3. Determination of Partial Derivatives of the Doppler-Equation 5 3.1 Evaluation of Derivatives of Doppler Frequency With Respect to Station Coordinates 5 3.2 Evaluation of Derivatives of Doppler Frequency With Respect to Parameters Influencing the Orbit 8 3.3 Evaluation of the Derivative for the Bias Parameter 9 4. Tropospheric Refraction Model 10 5. Representation of the Geopotential 14 6. Partial Derivatives of the Potential W 17 6.1 Gradient of the Potential U 18 6.2 Second Derivatives of the Potential U 19 7. Variational Equations 24 8. Orbit Integration Procedures 31 8.1 Runge-Kutta-Gill Method 3 3 8.2 The Adams-Bashforth Method 34 8.3 The Adams -Moulton Method 36 8.4 The Stormer Method 37 8.5 The Cowell Method 37 8.6 Integrating Procedure for Variational Equations 3 8 8.7 The Lagrangian Interpolation Method 39 in Page 9. Forces Besides the Earth's Gravity Field Acting Upon a Satellite 41 9.1 Attraction of the Sun and the Moon 41 9.2 Direct Solar Radiation Pressure 46 9.3 Atmospheric Drag 4 8 10. Adjustment of Observations 52 Acknowledgment 6 References 61 IV Computational Procedures for the Determination of a Simple Layer Model of the Geopotential from Doppler Observations Bertold U. Witte ABSTRACT. The geopotential of the earth is rep- resented by the potential of a simple layer dis- tributed over the earth's surface. Density val- ues of this layer for 104 surface elements have been determined from Doppler observations. Com- putational procedures using this method are pre- sented here. Numerical integration is employed to compute satellite orbits with time sequenced Doppler observations given as input. The geo- potential forces acting upon the orbit are de- scribed as well as the attraction of the sun and moon, solar radiation pressure, and air drag. The results of the orbit fits are combined in an ad- justment to determine the density values of the 104 surface elements. 1 . Introduction The determination of the earth's gravity field from Doppler observations using a simple layer model follows an approach outlined by Koch [1968]. The first results ob- tained by this method, from optical satellite observations, have been published by Koch and Morrison [19 70]. These results showed good agreement with existing solutions , so this approach was applied to other available data. Many stations of the worldwide satellite triangulation network [Schmid, 1969] are close to the Doppler tracking sites of the U. S. Navy Doppler Tracking Network (TRANET) [Anderle 1965] and connected by local survey to the BC-4 stations of the worldwide triangulation network. In order to combine the results of both systems the Doppler observa- tions are as a first step processed by means of the simple layer potential. Results are given by Koch and Witte [19 70 ], whereas the computational procedures are described here . 2 . Observational Equations Let r be the vector of the apparent position of the satellite and r = [x,y,z] (2-1) where x, y, z are the geocentric rectangular coordinates of the satellite; the origin is the center of mass of the earth; the z-axis is identical to the instantaneous axis of the earth, and the x-axis points at an angle east of the true vernal equinox which equals the precession and nutation in right ascension since 1950. Since the orbit computations extend over arcs not longer than 7 days this coordinate system is a good approximation to an inertial system. The rectangular coordinates of the station vector (u ,v ,w ) are given in a geocentric coordinate system whose o o o orientation coincides with the worldwide satellite triangula- tion network; that is, the w-axis points towards the mean pole 1900-1905 and the u-axis towards the intersection of the Greenwich meridian (i.e. the zero meridian of the Bureau International de l'Heure UTl-System) with the equator. By m<5ans of the polar motion as determined by the Bureau de l'Heure and the sidereal angle as defined by the Smithsonian Institution [1966] the station coordinates are rotated into the x,y,z system. x y, = s u v w S = cos9 -sinB -a cos6 - b sin sin6 cos9 -a sin6 + b cos a -b 1 (2-2) where a, b are the components of the polar motion, is the sidereal angle, a simple linear expression of UT1 , or = 0. rev 277987616 + 1 . rev 0027 3781191 (MJD-32 28 2 . ) , MJD is the Modified Julian Day Number as defined by the Smithsonian Institution. The vector r of the apparent position at a time t is a function of the vector e of the orbital elements at the epoch t , the vector x of the parameters of the gravity field, the vector of the station coordinates r , the s radiation pressure parameter K R over the observation interval and the drag parameter CL over the same time interval. r = ?(e(t ),x, ? s , K R (t), C D (t)) (2-3) For the formation of the observation equations the following expression is used, which is derived from the Doppler equation. f, , f = f K - — -Sr ( |? - ? | ) + 6f . (2-4) b c dt ' s ' tro Here f is the calculated frequency, f, the base frequency, c the velocity of light, which is 299,792,500 m/sec [Fricke et al. 1965], -rr-lr - r ! the rate of change of the distance ' dt ' s ' & between the satellite and the observation station, and <$f the tropospheric refraction correction. A detailed de- scription of the model used for tropospheric refraction is given by Witte [1970], but a condensed version is included as chapter 4. The observed Doppler frequency f is obtained by comparing the incoming satellite signal with a reference signal generated from a local precision oscillator and measuring the time required to count a preset number of beats. In order to get the calculated frequency f(eq (2-1)) the base frequency f , , which is the es- timated satellite oscillator frequency, is used. The assumed value of f, may contain an unknown systematic error which in this report is considered constant over one pass. The observed frequencies are influenced by this systematic error, which is introduced into the adjustment as a bias parameter. r. ■*■ ■*■ "*■ Approximate values for e, x-> v •> K„ , C~ are available so that we obtain by means of a Taylor series expansion the follow- ing observation equation, where a bias parameter b for the frequency offset of every base frequency for each pass is added Af. = I (|£) Ae v + 1 E 4 (||) A X , + I l^-) Ar 1 k=l V3e7 ik k ° £ = l V ^i£ l m3l V8r s ; . s mq lmq ^ i o i i denotes here the ith observation j denotes here the jth pass q denotes here the qth station (total of n stations) k denotes here the kth orbital element I denotes here the ilth surface density element m denotes here the mth coordinate of r g o denotes here the oth orbit Af = observed frequency f minus calculated frequency f. 3 . Determination of Partial Derivatives of the Doppler-Equation 3.1 Evaluation of Derivatives of Doppler Frequency With Respect to Station Coordinates. 8 f For the partial derivative ~ using equation (2-4) sm 9f b 3(|r - r R J * ) (3-D Here 9r c 8r sm sm d( r - r ) -> -*■ 1 • s ' r - r s ' dt and the slant range is I? - r I = p = /(x - x ) 2 + (y - y ) 2 + (z - z ) 2 (3-2) 1 s ' s J J s s The differentiation of the slant range with respect to the time yields d/Cx - x ) 2 + (y - y ) 2 + (z - z ) 2 s J J s s dt rtx dv dz dx 9 dy q dz c= 2/x-x ) 2 + (y-y ) 2 + (z-z ) 2 s J -'s s (r - r ) • (r - r , s s ^_ (3-3) l"*" ■*" l r - r 1 s ' Introducing p as a unit vector the change in the slant range yields with equation (3-3) the following relation r - r I" = p-(r - r ) (3-4) s ' s From this we get for equation (3-1) 3f 3r sm 3p 3r sm ,1" -t V (r - r ) - p 3r 3r sm where 3r 3( r - sp- ar sm 3r sm (? - 3r ? ) s sm r - 3r 3p 3r sm 3r sm r - r 3r 3r • c? sm ? ) s C? - r ) s r - r (3-5) with 3? 3x 3r c 3? 3z and 3r c 3x~ -uy f L 3? c (OX ' 3r < 3z~ where co is the angular velocity of the earth. Finally we get 3f 3r sm 3r _ 3? f, r/ 3r D C sm .- + 3r sm • (r - r ) s (r - r s ) L r - r r - r -).Cr - r )-p. ^^J sm (3-6) 3.2 Evaluation of Derivatives of Doppler Frequency with Respect to Parameters Influencing the Orbit. The derivatives with respect to the parameters e , Xo» ^ and C , as used in equation (2-5), are now derived. We set 3f 3p K m , 3X 3p = 1 m r 3X 3 a . 3X m „ 3f __m m=l 3X 3p K (3-7) K m where p„ stands for the above mentioned parameters e ^> Xn? Kt, and C~; and X for the coordinates x, y, z of the K x, R D m ' J apparent position of the satellite. The partial derivatives 3X 3 X ~ and g m are computed with the help of variational equations 3p K 3p K * M (section 7 ) . The derivation of the Doppler frequency in equation (3-7) with respect to the apparent position of the satellite leads to r D 3f 3X f 3(|? m 3X (3-8) m Applying (3-4) 3f f, 3(p • (r - r ) ) b s 3X m 3X m (3-9) 3f 3X m ■-*( b ( 3p ft t x . . Cr - r ) + 3r 3X m w m Analogous to equation (3-5) and with 3r 3X we get 3f 3X m m 9?_ 3?_. (?-?)(?- r ) 3X ax V s ; ^ r s ; m m (3-10) r - r r - r (r - r ) s (3-11) The differentiation of the Doppler frequency with respect to the velocity X , where X stands for the coordinates x, y,z J m m leads to 3f 9X m f ^ d ( I r - r I * ) b ' s ' m f. 9(p • (? - r )) D S 9X m 8? ax (3-12) (3-13) (3-14) m dr 9r For the components of the derivative -^— and ~ — we use m 9X m 9r 8x 3r ' 8y -*- 1 8r 9 z (3-15) and 3x 8r 97 8r 31 (3-16) 3.3 Evaluation of the Derivative for the Bias Parameter From equation (2-4) we get for the partial derivative of the Doppler frequency with respect to the base frequency f. 3f 1 ,■> -*■ ttp- = 1 - — r - r df, c ' s (3-17) where |r - r | is obtained from equation (3-3) 10 4 . Tropospheric Refraction Model For the computation of the correction term <5f in equation (2-4) a tropospheric refraction model derived by the Naval Weapons Laboratory is used. It is preferred to models developed by Hopfield [1963 and 1969] because the mean error of unit weight resulting from the use of this formula, which will be derived below, is somewhat less than that ob- tained with the Hopfield models [Witte 1971D The range error As caused by the tropospheric refraction may be described by the following equation As = Jn ds - / ds (4-1) tro g where n is the index of refraction for radio waves and varies along the signal path. The first integral is taken along the signal path through the troposphere and the second along the geometric path. If both ways are set equal, only an error of second order occurs As = / (n-1) ds (4-2) g The contribution of the tropospheric refraction to the Doppler shift of the signal received from a passing satellite is given by (compare also Hopfield 1963) = _ £b d(As_). (l+ _ 3) tro c dt 11 In equation (4-3) As is differentiated with respect to the time t. For the approximation of a flat earth we get with equation (4-2) and (4-3) H f J(n-l) dh 5f = - — — (° tro c dt cos Z ) (4-4) The angle Z is the zenith distance and H the height of the satellite above the earth. After integration and differentia- tion we get for (4-4) f x(pxH) px(pxH) q = P XrpXH) (4-7) where |px(p xH) | = sin Z (4-8) ->- p in equation (4-6) can be rewritten as -> -f -f P = r - r (4-9) and p = A (u - io) Inserting these expressions (4-7), (4-8), (4-9), (4-10) in equation (4-6) leads to ->- cos Z s sin Z and substituting px(pxH) = (p-H)p - (p-p)H = p cos Z - H into (4-11), equation (4-5) becomes f Sf. = - -£ (r - r )-(p L) (n - 1) (4-12) tro c s H cos Z We replace (n - 1) by ., tro AN, - N . , n - 1 = -4- / (n-l)dh * — ^- ^^ (4-13) H h . . (r - r )-H 1 ' stat s 13 h , . = height of station stat & h^ = height of troposphere tro to r AN = total change in refractive index from sea level |! ' ; "' up to the top of the troposphere scaled to height. AN = change in refractive index from sea level up to the height of the station scaled to height. AN , . is determined by means of an exponential model s ~ca~c similar to that described by [Bean et. al . 1966]. C AN . . = C (l-e" C l h stat) + — (4-14) stat ° cos 2 Z C 2 In order to get flat observations corrected, the term cos 2 Z is added. Observations with Z>85° are deleted, because the computed corrections indicate that for this data the model used is overcorrecting. Combining equation (4-13) and (4-12) leads to the results given in the following expression f. ( _ r-r \ AN^_ - AN . . of. = - — 1 r - r - — [ = (4-15) tro c ' s ' !■*■ i r, [ p cos Z I r cos Z J 1 s ' For the constants C , C, , C_ and AN, the following numerical o 1 2 tro 6 values, which are used by the Naval Weapons Laboratory (NWL), have shown good results C = 0.0025525 o C 1 = 0.1238 C 2 = 0.657X10" 5 AN^ = 0.0 02 3 tro These numbers were used for the computation of the tropospheric correction model as given by equation (4-15). "Private communication from R. J. Anderle [1970]. 14 5 . Representation of the Geopotential The gravitational potential W of the earth is divided into the normal potential U, which is taken as known, and the disturbing potential T, which will be determined by the approach used by Koch [1968], W = U + T (5-1) with 7 n U = »i r 1 + E E C-) n P (sine})) • ( C cos mX n n r nm \ nm n=2 m=0 + S sin mX nm + i oo 2 r 2 cos 2 (5-2) r, , A are spherical coordinates in an earth-fixed coordinate system, k is the gravitational constant, M the mass of the earth, a the mean equatorial radius, P the fully normalized ' H ' nm J associated Legendre function of degree m and order n, and co the angular velocity of the earth. C and S are fully nm nm normalized harmonic coefficients. Their values are taken from the solution of Anderle [1965]. The coefficients C2 1 and S2,i have been deleted because of the chosen coordinate system, defined above. The centrifugal term in (5-2) is omitted for points outside the earth. The disturbing potential in equation (5-1) is represented as the potential of a simple layer with the density x C <|) , X ) 15 distributed over the surface E of the earth T -_ jj Xiiiil dE (5-3) E S s being the distance between the point at which T is com- puted and a variable point on E. To evaluate this integral, the earth's surface is divided into elements on which the density of the simple layer is assumed to be constant so that the integral in (5-3) can be replaced by a summation. 104 surface elements AE, are chosen here. These are bordered by parallels and meridians and approximate the size of a 20° x 20 area at the equator. 104 AT T = I Xo // ~ (5-4) £=1 AE £ I The integral over the surface elements AE p in (5-4) is computed numerically by splitting AE p into 4 subdivisions and assuming the kernel of the integral computed for the midpoint of the subdivision to be constant over the subdivision To compute the distance s between a satellite and the midpoint of the subdivision of AE. it is necessary to determine the coordinates of these midpoints. For this purpose the equi- potential surface U(eq 5-2) is defined as an ellipsoid U . U is computed from the equation which holds for the surface 16 of a level ellipsoid kM U, f 2 ,2,1/2 (a - b ) arc tan (a 2 - b 2 ) 1/2 + | co 2 a 2 (5-5) All the parameters of (5-5) are available from (5-2) except b. We obtain b from b = a(l-f), and f from the value of C20 in (5-2) in the following manner. The relation between C 2 o and f is given by [Heiskanen and Moritz 1967, p. 110]. C20 = - J2 / *$ C5-6) t 1 e* 2 J 2 ~- 3 1+e '2 , _2_ me ' 15 q D with and 1 q Q = 2 1 + — rsrj arc "t an e ■1 m = 2 2 (0 a kM/l+e' 2 ' 2 - 1 " (1-f) 2 - 1 e' is obtained from (5-6) as a function of J 2 by the Newton-Raphson method of successive approximation. For dJ : de this, the derivative -3 — f is needed dJ 2 1 / > „'K '), ■/< de' " 3 l 2e' 3 nrv2 2 T 2 15 qoAi + e .2 Cl + e' 34 )^" 45 IT? 7 " 2 ^ q kMv / (1 + e ,2)3 m e' 2 a) 2 a 3 me 2q, ((-i^)r^ 6 4. 1 x 3 pr-3- arc tan e ? + p-jr (5-7) 17 Now all parameters for the computation of U ( eq (5-5)) are known and therefore the radius vector of the reference sur- face can be computed in setting U = U. Solving equation (5-2) for r we get r = kM 7 n 1 + Z Z (-) n P (sin*) x _ _ n r nm T ■*■ n=2 m=0 (C cos mX + S sin mX) nm nm co 2 r 2 cos 2 (j) 2U. (5-8) Equation (5-8) is determined for the earth-fixed co- ordinates of each midpoint by direct iteration for the reference surface, where the initial value for r is obtained for a given latitude (J) on an ellipsoid of revolution with semi-axes a and b. To these values of the radius vector the topographic heights above sea level, published by Kaula et al . [1966], were added. Employing the given values , X of the midpoints, we can compute the corresponding u, v, w coordinates, and then the x, y, z coordinates of the midpoints by means of the rotation matrix (2-2) 6 . Partial Derivatives of the Potential W The derivatives of the force acting upon the satellite with respect to the position of the satellite are needed in the solution of the variational equations (Sec. 7) for the orbit computation. These derivatives are 8 2 W 3 2 W 8 2 W 8 2 W 8 2 W 3x 2 ' 3x3y ' 3x3z ' By 7 ' 8y3z , and 3 2 W The gradient of the potential U is used for the computation of the reference orbit. The approximate density values are set equal to zero. Therefore the gradient of the potential T will not be computed here. 6.1 Gradient of the Potential U Using the chain rule method for the partial derivatives of U, where U = U(r, cf>jA), we find 3x 3r 3x 8cj) 3x 3X 3x (6-1) M £uto + 3U8£3U8x 9y ' 8r 8y 8(f) 3y 9X 8y (6-2) m min + mii + MiA 3z ' 3r 3z 3 3z 3X 3z (6-3) The differentiation of equation (5-2) with respect to r yields 3U 3r kM r " " r^L 7 n 1 + Z Z (n+1) (-) P (sin) x o n r nm n=2 m=0 (C cos mX + S sin mX ) nm nm + w 2 r cos 2 cJ> (6-4) d P (sincf>) With nm d(J> = -m tan d> P (sin d>) + P ,, (sin ) Y nm Y n,m+l 19 we obtain 1 r\ |x = - — [ E 2 C-) (m tan cj> P (sin cf>) - P ... (sin <(>)) ^ r L n =2 m=0 r nm n > m+1 (C cos mX + S sin mX ) - w 2 r 2 sin d> cos d> (6-5) nm nm J and £E = - — [ E S m (-) n P (sin cf>)(C sin mA - S cos mX) 8X r o „ r nm T nm nm J L n=2 m=0 (6-6) As mentioned earlier the centrifugal terms in C6-4) and (6-5)' should be omitted for points outside the earth. The other partial derivatives — equations C6-1) - C6-3)--are found with r = A 2 +y 2 + z 2 , tan <£ = . = , and /x 2 + y 2 tan X = J 8x " ' r* 8y " r s 3z " r ; M xz 8jb_ = -yz 8_£ . /x 2 + y 2 ,. . 3 * " " Wx 2 J y 2 ' K " Wx 2 + y 2 ' 8z ' r * C ^ |4 = — ^ , |i = ^ , |i = (6-9) 3x x 2 + y 2 9y x 2 + y 2 8z 6.2 Second Derivatives of the Potential U The second derivatives with respect to x, y z of the 20 potential U yield (using equations (6-1) - (6-3)) the following expressions 3 2 U 3 2 U /9rV 9x- 3r 2 + 2 3x 3 Z U /8(f) 3<}> 3x 8_nj 3A 2 3x 3 2 U 3r 3_£ + 3 2 U 3r 3A_ + 3 2 U _3J>_ 3_A 3r • 3 • 3 A 3x 3x + 1H ilz + 1M ill + 8U ilA 3r - 2 3x 3 9 rs 3x 3 A _ 3x C6-10) 3 2 U 3. 2 U 3r 3r 3^U 3^ 3^ 3JU 3jV_ _3J. 3x3y j, 2 3x 3y 3c}) 2 3x 3y ' -.2 3x 3y 3 2 U 3r3(J> 3 2 U 3A3c}) /3r ii + 3r 3£\ + 3 2 U /3r _3_A + 3r 3_A\ \3x 3y 3y 3x/ 3r3A\Jx 3y ' 3y 3x / / 3A_ ii + il d±\ \ 3x 3y 3y 3x / + 3U 3 2 r + 3U 3 2 (j) + 3U 3 2 A 3r 3x3y 3$ 3x3y 3A 3x3y (6-11) 3 2 U . 3^U 3r 3r _3_^U 3£ 3£ + _3_fU 3_A 3_A 3x3z 3r 2 3x 3z 3<1) 2 3x 3z 3A 2 3x 3z 3 2 U 3r3(}) 3 2 U 3A3(j) / 3r 2i + 9r 3 $ \ 3 2 U /3r W , 3r 3_A \ ^ 3x 3z 3z 3x / 3r3A^3x 3z 3z 3x ) ( 3_A d± . _3A Z±\ \3x 3z 3z 3x/ , 3U 3 2 r . 3U 3 2 4> , 3U 3 2 A f 6 -12) 3r "SxTz 3^ ?x3~z "3~X "3~xTz 3 2 U 3y : 3 2 U / 3r X2 3r- 3y 3 2 U 3<}) 2 ii 3y 3 2 U 3A 3A^ 2 3y, + 9 2u lz ii + 9 2 u 9r li + 32u 2* ii 3r3(J) 3y 3y 3r3A 3y 3y 3X3$ 3y 3y 21 + £U lf£ + £U 3J^_ + 3U 3J\X " " ' H 3y 2 3r 3y 2 3X 3y 2 C6-13) 3 2 U _ ^U 3r 3r SJ^U 3f 3| + S^U 3X. 3X 3y3z " 3r 2 3y 3z 2 3y 3z 3X 2 3y 3z ^!£_ f i£ i£ + 9r 3^\ 3 2 U / 3r 3A , 3r 3X\ 3r3(J> \3y 3z 9z 3y / 3r3X \ 3y 3z 3z 3yJ 3 2 U ( IX d± . d\ d±\ i 3U 3 2 r + 3U 3 2 ({) + 3U 3 2 X (6-14) 9A9 \ 3y 3z 3z 3yJ 3r 3y3z 3$ 3y3z 3X 3y3z 3 Z U 3 2 U 3z ; 3r' + 2 3 2 U 3r 3 3r3<£ 3z 3z ^ + 2 3 2 U 3r 3X ~ ^ + 2 3r3X 3z 3z 3 2 U 3_X 3^ 3X3(b 3z 3z 3U 3 2 r 3U 3^1 3JJ 3 2 X 3r 3z 2 3d) . 2 ' 3X _ 2 T 3z 3z (6-15) Now the partial derivatives in equations (6-10) - (6-15) are derived. With equation (6-4) we find 3 2 U , 2kM 3r 2 r 3 1 n 1 -: ;, Z (n + l)(n+2) (-)" P (sin 4) 2 - r nm r n=2 m=0 (C cos mX + S sin mX ) nm nm (6-16) I^U = _ kM (n+1) ( a 5 ( p- (sin 4,) 3r3d; 2 o -n r n,m+l T r Ln= 2 m=0 m tan d> P (sin 6) ) T nm r (C cos mX + S sin mX ) nm nm (6-17) 22 9r9X 9 2 U vm r 7 n ^ n - - SI E E (n + 1) (|) P^ (sin E E (-)"( - - -- P (sin cj)) + m tan (P m+1 (sin (J)) - m tan P (sin )) + (m+1) tan P« mil ( sin ) " ^ xo ^ sin ))CC cos mX n,m+l T n,m+2 r J nm + S sin mX) nm (6-19) 9 2 U 9 P (sin ) o n 2? nm v _n=2 m=0 - P ,, (sin d)))(-C sin mX + S cos mX) n,m+l r nm nm (6-20) 3 2 U 9X 2 kM r " 7 n n E E m 2 (§) P_ (sin <|>)(C,._ -n=2 m=0 r nm nm + S sin mX ) nm ] (6-21) And for the other terms in equations (6-10) - (6-15) we get 9 2 r _ 1 x 2 . 9 2 r _ xy . 9 2 r _ xz „ 2 r 3 ' 9x9y 3 ' 9x9z 3 9x v J r r (6-22) 9 9x 2 r 2 /x 2 +y = /_ 1 + 2xi + _xf_\ +v 2 V r 2 x 2 +v 2 ' r>z x z +y 23 9 2 (j) _ xyz / ,2 + 1 \ 8 2 (J) _ -x + 2xz 2 (6-23) 3x3y r 2/ x 2 +y 2 \ T 2 x 2 +y 2) ' 9X8Z r 2/ x 2 +y 2 T » /^T^T ^ - 2 ^ ; *ii = y 2 -* 2 ; ^fA = (6 _ 24 ) 9x 2 (x 2 + y 2 ) 2 dxdy (x 2 + y 2 ) 2 9x8z 9 2 r 1 y^ 9 2 r yz ^2 ^ 3 ' 9y9z 3 9y r J r 9 2 4> : z / 2y 2 _ . + y 2 (6-25) 2 r 2 /? = /2Yl _ 1 + _r_\ + v 2 V r 2 x 2 +v 2 / 8y^ r"/x"+y" V r 2 x z +y 9 2 . _ Y + 2yz 2 8y3z r 2 /x 2 +y 2 r*^ 1 ^ (6-26) ill ■ 2xy__ . _^A = Q (6 _ 27) ^2 / 2 i 2\2 9y 9 z 9y (x +y ) 9 2 r 1 _ zj_ 9j^_ _ 2/x 2 + y 2 z 9 2 A _ (6-28) 9z 2 r r 3 ' 9z 2 r" 8z 2 This method of computing the partials of U was found by Gulick [1970] to be the most efficient compared to two other alternatives. Subroutines which compute these partials based on this and other methods may be found in this reference along with a de- tailed discussion of comparative efficiencies and a complete derivation of expressions given in equations (6-10) - (6-28). 24 7 , Variational Equations The position of a satellite can be described as a function of the parameters of the earth's gravity field plus non-conservative forces F c r = VW(r,x,t) + F (r,r) (7-1) For the time being the non-conservative forces F will not & c be considered here so that we have r = VW(r,x,t) (7-2) For the solution of this differential equation a numerical integration procedure, as described in the next section is needed. In order to do an orbit adjustment, additional equations, called variational equations, are necessary for 3X 9X the determination of the partial derivatives ~ and ~ 9p k 8p K given in equation (3-7). To make the notation more concise, * 8t = (?,?) and r 0st = Cr 0J ? , we form state vectors r__,_ = (r,r) and r = (r ,r ) . r consists of the quantities X , X . r n constitutes a subset n mm °st of the parameters p, . Here r is the initial position vector ? ... -> -f and r the vector of the initial velocity; r and r are equivalent to the six orbital elements represented by e, in equation (2-5). Using the state vector notation we get for eq (7-2) ? st = (r\VW) = S 8 (r 8t ,x»t) (7-3) 25 The position of a satellite is now given as a solution to the differential equation (7-3) with the initial state vector as boundary condition. If we form the time derivative of the partial of r . with respect to the parameters p. r St -^ K we get J , d [~5l V% 1 ^st (7-4) dt " .-*• The numerical integration procedure will be applied to (7-4) in order to obtain the required partial derivatives. First we will determine the above mentioned time derivative of the partial of r with respect to the initial state vector and get, with eq (7-3), V3r °2i / = 8p »t 3r st , 3 *st 3x , 3 *st at dt a + -HJ + --* 9t 3r . 9r„ , 3y 9r . dr , st °st A °st °st Note that the last two terms are zero, because x an ^ t do not depend on the initial state vector. So we finally get / 9? st \ \ o st /_ = st st (y _ 5) dt St °St 3r_ 9r , 26 Using the rectangular coordinate system x, y, z, defined in section 2, we get in matrix notation for the first factor in eq (7-5) 9r st 9r st 9x 9x 9x 9x 9x 9x 9x 9y 9z 9x 9y 9z 9? 9x 9y 9y 9y 9z 9y 9x 9y 9y 9y 9z 9z 9z 9z 9z 9z 9z 9x 9y 9z 9x 9y 9z 9x 9x 9x 95^ 9x 9x 9x 9y 9z 9x 9y 9z 3£ 9x 9y" 9y 9z 12 9x 9£ 9y 9£ 9z 9z 9z 9z 9z 9z 9z 9x 9y 9z 9x 9y 9z (7-6) If we partition this matrix, we will get the following structure 9r st 9r st 9 2 W 9X m 9X. 9X. ; ^ (7-7) where is a 3^3 zero matrix and I is a 3x3 identity matrix -*■ "*■ (since r and r are the independent parameters of the differential equation (7- 3)), the lower left hand portion arises from differentiation of (7-5) with 9 2 W as a shorthand notation 9X.9X. i ] 27 for this part, and the lower right hand portion is if drag is not present or is neglected; terms from drag would occur in the entire lower half of the matrix if drag were present. If a simple radiation pressure model is included, a Dirac delta function term occurs in the lower left portion of the matrix; for a very complex model with re-radiation an impulselike component would still appear. In setting up the variational equations for the solution of x an< 3 r , the nonconservative forces F are not being st -" <^ 9X applied. Hence, we set the term ^ £ n equation (7-7) 9X m to zero, and get instead of (7-7) - 3r st I 3? . st 9 2 W 9X.9X. (7-8) The differentiation of the state vector r . with respect to st r the initial state vector r yields st 9x 9x, dx ax 9x 9y< 9z, dx, dx 9y c 9x 9x 9y Q 9z 9x 3y c 9z 9y_ 9x 9y ay 9y 9z 9y 9x 9y 3y Q 9y 9z 3? t st 9z dx 9z 9y 9z 9z 9z ax 9z 3y c 9z 9z st dx 9x dx 3y Q 9x 9z 9x 9x 9x 8y G 9x 9z 9y 9x 3y 9y D 9y 9z 3y 3x 9y 3y 3y 3z 9z d± 9z 9z 9z 9z 9z, (7-9) 28 For the start of the integration procedure this matrix has the following structure 3r st 3?, st (7-10) For the time change of the derivatives of the state vector r with respect to x we get, analogous to (7-5), st 3r st X dt ^st 9? st 3? 3r st st ax ax (7-11) with 3? st 3X 3x 3X 3x 9-Xi 3X2 3y 3Xi 3y 3X2 3z 3z 3Xi 3X2 3x 3x 3Xi 3X2 By 3Xi 3y 3X2 3z dz 3X2 3x 3Xi 04 3Xi 04 3z 3Xi o k 3x 3X 10 4 *± 3Xi 04 3z 3Xi 04 (7-12) 29 3r st w o 3x 3Xi 3y 3Xi 3z 3Xi . . . . . . . . . . . . 3x 3X2 ' * 3x 3Xi ot ay 3X2 ' • 3'y 3Xi o 4 3z 9z 3X 3X 1 h (7-13) 3r and st supplied by (7-8 ) 3r st 3r In (7-11), st 3X represents the explicit instantaneous de- pendence of velocity and acceleration on the gravity field. Since instantaneous velocity is an independent variable, representing, with the instantaneous, position the instantaneous state of the satellite, the top half of (7-13) is null. On the other hand, 3r st 3X represents the dependence of the state vector (position and velocity) on the gravity field, considering the state vector as obtained by integrating the accleration over time; hence the terms of (7-12) are not necessarily zero, although the lower half of (7-12) literally appears to be the same as the upper half of (7-13). Initial conditions for the integration are st i = I (7-14) o L -i For the elements of matrix (7-13) we find 3X, 3x 3X £ 3y 3X 3z 3X c3?) 3X 3* •I 3 cli) 3x 3y ■i (P) 3X^ 3z (7-15) (7-16) (7-17) 30 where T is the disturbing potential as defined in eq. (5-3). Using the distance s. = Ax - x^) 2 + (y - yJ 2 + (z - z.) 2 between the satellite and the midpoint of the £th surface element we get for the partial derivatives of T with respect to x, y, z from eq (5-4) 104 S £=1 AE, |£ = I X, !! " ~^dE (7-18) 9T 'I 104 *~ x z dE s 3 y-y* dE s 3 Z " z £ dE s 3 3v S X £ // - -5T* dE C7-19) dy £ = 1 * AE £ ~ T 104 f^ S Xo // " "TT 1 dE (7-20) Z £=1 * AE £ For the equations (7-15) - (7-17) we now get x-x. af ( ^ } = // " ~^ dE (7 " 21) dX £ dX AE £ s if- ( H ) = // " ^ dE (7 " 22) dX £ dy AE £ S -JL (II) = // jA dE (7-23) d *£ dZ AE £ 31 8 . Orbi t Integration Procedures The equations of motion and the variational equations of the satellite expressed in the rectangular coordinate system defined in sec. 2 are integrated numerically with an 0.8 minute time step. The twelfth-order Cowell-Stormer multistep process is used for the positions and the tenth-order predictor corrector formulas of Adams-Bashforth and Adams-Moulton re- spectively are used for the velocities [Henrici 1962, pp. 191- 198, 291-294], The accompanying flow-chart shows the struc- ture of the above mentioned procedures. The derivatives with respect to the initial state vector and to the unknown density values (i.e.? eq (7-5) and (7-11)) are obtained numerically by integrating the variational equations. In this integration the corrector is solved ex- plicitly without predicting, using the linearity of the variational equations [Riley et al. 1967]. These procedures require the use of a starting process, because the first time steps cannot be computed without knowing previous values. Gill's modification of the Runge-Kutta in- tegration method [Romanelli 196 0] is therefore employed for the computation at these time steps. The methods mentioned are now explained. 32 PREDICT ? p+ i ADAMS- BASHPORTH PREDICT ? p+1 STQRMER 1 = COMPUTE ?p+i CORRECT ? p+1 BASED ON CURRENT m-* ADAMS -MOULTON ? p+l ? p+l >< COMPUTE CORRECT r p+1 Vp +1 SETS — p COWELL yPp+1 Tp+1 A>+1 V p + 1, Figure 2. Integration process for the computation of position, velocity, and acceleration 33 8.1 Runge-Kutta-Gill Method The totality of the variational equations is, as has been shown in section 7, a 6x(6+104) matrix with a=aCr . ,x), where a is a subset of parameters of p. 3(r st> 8(? st ) 8(a) 8(r 0st ,x) 3r 3r . st st 3r, st 3X (8-1) We now apply the Runge-Kutta-Gill integration formulas to determine 9 (r ,)/3(a) from its time derivatives, given by (7-5) and (7-11) K = t & Y dt p R x = $ (K - 6 Q ) P P > Y _, n = Y + R .- p+1 p p+1 Q 4. n = Q + 3R ,. + e K X D + 1 X D D+l D •p + 1 P + 1 P P J P = 0,1,2,3 (8-2) After each cycle an intermediate value for d/dt[3Cr , )/3(a)] J st is computed from (7-5) and (7-11) using Y for 3(r . )/3(.a) ° p st (i.e. for 3(r , )/3(r .) in (7-5) and 3 (r J .)/3(x> in (7-11)). st u st st A The value of 3 (r . )/3(a) for t + At is then Y u . The initial st H value for Q = o is a 6x (6+104) zero-matrix and for Y = Y p u p ° the analytical derivatives at the starting time t as given in section 7. 34 5 , , , 5 , e are p+1' p' p p b p 6 P E P 1/2 2 -1/2 1 1 - A/2 1 A/2 J 1 2 1 + A/2 1 -1 - A/2 3 1/6 2 -1/2 (8-3) This method is also used for the computation of the starting values for the components of the velocity and position-vectors with Y now as the initial state vector as obtained for the different satellites from the Smithsonian Institution. Q again is a zero vector. For both computations a time step At of 0.1 minute is used only in the starting procedure 2 The Adams-Bashforth Method Here we have t x L *-k+l 'P+1 = / P(t)dt t P = At X E Y V y m=0 •• LzJ (8-4) where the constants y are independent and will be calculated m r numerically below. At is the time step; q is the number of steps, which in this case is 10; and P(t) is the interpolating polynomial over the interval [t ,t ]. 35 From the Runge-Kutta-Gill method initial values for the state vector are supplied. Using these values X X X y 5 y 5 y .z_ _z_ _z J q-i q-2 x y l_2J are computed from the equation of motion so that the differences V m x y LZ, can be calculated from the relations ,m X X X y = y - y _z_ p _z_ P _z_ X / X y = v| V- 1 y z p \ z p-1 (3-5) The expression on the right hand side of equation (S-M-) is thus known, and can be calculated for p > q. After p+1 using the corrector formula, explained in sec. 8.3, the index p is increased by 1, and the same formula is used to calculate x y , etc. With the help of the following recurrence L 4 J p+2 formula, which is derived in Henrici [1962, p. 193], numerical values for y were calculated 'm i i : Y m + 2 Y m-1 + 3 Y m-2 + * * * + m+1 'm-q Y T = 1, m = 0,1,2, . 36 8.3 The Adams -Moulton Method Whereas the Adams -Bashforth method is a predictor formula the Adams -Moulton method is a corrector formula. Here we have x y , -- j p+i "- p + l q = J + P(t)dt = At I yZ V t m=0 m "P x y L_Z_I (8-6) p + l where P(t) is the interpolating polynomial over the interval A [t ,. . t ,,]. The recurrence formula for y is similar to p-q + 1 p+l 'm that one for y and can also be found in Henrici [1962] on m page 194. Formula (8-6) is used like the Adams-Bashforth formula (8-4) except that now only the values are known. It is employed to redetermine x y z L. _jp p-q + 1 an iterative procedure, where an approximate value x y * LJp + 1 p-1 in x has been obtained using the results of the predictor formulas (Adams-Bashforth and Cowell-Stormer ) (0) x -,m (1) x y L*Jp + l Calculating P (1) p+l + At E y v m=0 X y l-Z_ (8-7) p + l and re-evaluating the differences, a better values is then obtained for p+l 37 8 . 4 The Stormer Method By integrating a differential equation of the second order y"(x) = f(x,y,(x)) twice, we obtain finally the formulas for the computation of the position vector - 2 p+1 L J p-1 = At 2 E a V n m m= x y L-Z_ (8-8) where again the above used symbols are defined as before. a is obtained with the help of the following recurrence m to formula 2 2 = 1 - x h„ a .. - Tr- h a - m 3 2 m-1 4 3 m-2 with h =l+i+...+i;m=l,2 m 2 m ' 2 h . , a n ;m=l , 2 m+2 m+1 °' Formula C8-8) is used in the same manner as the Adams-Bashforth method. 8.5 The Cowell Method Here we have X X X y - 2 y + y z z -r-, -\ z p-1 = At p-2 1 a" V m m=0 m where a m 2 .'• 2 * 3 2 m-1 4 3 m-2 X y 2 P 2 m+2 - h m (8-9) + 1 cr , m = 1,2, 38 8.6 Integrating Procedure for Variational Equations The system of variational equations (7-5), (7-11) can be written as R = SR + T where R = 8r , 3r . st st (8-10) (8-11) and is a 6x(6+104) matrix; 3r S = st 8? , and is 6x6 st and T dt st 3X and is also 6*(6+104). The equations (8-10) are integrated using the qth order Adams-Moulton formula. At the (p+1) step, r st r q * • R xn = R + At E 8 R ... P + l P p -o qp P+l-P (8-12) where 3 = (-1) qp C>> Y* ♦ C%h Y * p+1 + + <3> y! p = 0,1 . . . , q. q = 0,1 ... . Assuming that R and R up to the pth step are known from previous applications of this procedure (or at the start from initial values obtained by the Runge-Kutta-Gill method), 39 we see that the only quantity not known on the right-hand side is R ,,. Let us rewrite (8-12) as p + 1 q R .. = R + At E 3* R xn + At 3* R in p + 1 p 1 qp p + l-p qo P + 1 Substituting for R ... from (8-10), we obtain ° p+1 q R .- = R + At t 3* R xn + At 3* SR xn +At 3* T , p + 1 p qp p + l-p qo p + 1 qo and solving for R ,, , & p + 1 R xn = [i-At 3" s] ( R +At ("3" T + £ 3* R x . 1) p + 1 L qo J I p L qo p=1 qp p + l-pj[ The matrix in brackets to be inverted is 6x6, and the matrix in braces is 6x(6+10U). Following the computation of R .-,, ^ D + i can be obtained directly from (8-10) for use in the next step. The actual computational procedure breaks the matrix R into components as shown in (8-11) and follows the outline given by Riley et al . [1967]. 8.7 The lagrangian Interpolation Method. To interpolate between the time steps, Lagrange's interpolation for equidistant abscissas as given by Henrici [1964] on pp. 201 - 202 is applied. Because an equal time 40 step At is used, we can assume that the points, where the values (f, ) of the positions, velocities and variational equations are given, are equally spaced. Here we have n p(s) = E £ (s) f k=m where n l f O = IT s 9, V S ; J1 k-q q=m ^ q*k This representation of the Lagrangian polynomial p(s) is independent of the size of At. The functions )L (s), which Henrici calls the normalized Lagrangian interpolation co- efficients, depend only on s (the relative location of the observation time t with respect to next lowest time t, ) and on the integers m and n, which are the bounds of the set of interpolating points. It turned out that the results for n=4- and m = -3 were sufficiently accurate. (tj) -3 -2 -1 1 2 3 1 Using this method only the values f, for eight time steps and the value t, need be saved. The value t, is incremented by At for each step in the integration process, and the time- sequenced observations are processed as soon as the time of 41 the next observation (t , ) is exceeded or equaled by the obs current t, . We then have s = 1 - (t, - t, )At 1 obs The f, are stored in a circular fashion so that the latest values replace the former values for f_„, and we need only monitor the index of the current location of f_„ in the eight place array. ^ • Force s Besides the Earth's Grav ity F ield Acting Upon a Satellite In addition to the perturbations caused by the varia- tions of the earth's gravitational field as observed in the previous chapters, satellite orbits will also be perturbed appreciably by the gravitational attractions of the sun and moon, the radiation pressure of the sun, and the drag due to the atmosphere. 9.1 Attraction of the Sun and the Moon All the effects of lunar and solar gravitation upon an artificial satellite may be expressed by adding additional terms for the potential field through which the satellite travels. If r~ denotes the position of the sun or moon and 42 M„ its mass, the disturbing function may be written [Brouwer and Clemence, 1961, p. 465] AF = GM 2 ( 1 - - T -^-\ (9-1) \ |r 2 - r| r 2 / The first term gives the potential due to the disturbing body. The second term arises from choosing the center of mass of the earth rather than the center of mass of the entire system of bodies as the origin [Danby 1962]. If we differentiate (9-1) with respect to X , where X stands r m' m for the coordinates x,y,z, we get ||£ = GM,/- - Xm - + C ( ; 1 ; i- A (9-2) 3X m 2 \ |? 2 - ?| 3 m |? 2 - ?| 3 r 2 3 / w th ith Tr, - r„ (£,£;,£;), where L 5 Co lie in the equator wi C-, pointing toward the vernal equinox. If — 5- is factored r 2 out we get GM r U£ = -~ [-5 + C5 - x ) - 2 - y ] (9-3) 8X 3 m m m 1 -> -* 1 3 m r 2 |r« - r I with 3 „-v -> 3 P 2 r 2 2r * r 2 r ~2 = [1 + (- ) Z - - (£-)] ^ • (9-4) ■* -*-i 3 r>„ r r r-, r rt - r 2 2 2 43 -> ->■ If (— ) - ( — ) < 10" , equation (9-4) will then be r 2 r r 2 r 2 — ^ developed in a series so that we get P 2 * n 3 r ,r N 2 2r ' r, 2 , r Nn , 15 r ,r N 2 2r ' r, 2 ,r N -,2 1 . 4 [(JL)^ - -^ (= £)] + J£[(JL)' £ (JL)] Ip _ pi 3 2 r 2 r r 2 r 2 8 r 2 r r 2 r 2 . 35 [(JL) 2 . ^2 jj. ]3 + 315 t ±2 _ £^2 _r ,-,4 (9 . 5) 16 r 2 r r r 2 128 r 2 r r 2 r 2 The quantity r-rJlr r„) is the product of the normal vectors pointing to the satellite and perturbing body. The scale of the perturbing body's orbit is given by a 2 and is needed in a few of the small terms, where r 2 = a~ (1 - e cosE). The series (9-5) is then applied to (9-3) to obtain the difference 3 r 2 -5+5. 'm m i -*■ ■*■ i 3 -p _ -p |r 2 r | 3 The coefficient GM/r 2 can be deduced from Kepler's third law. Applying this law to the orbit of the moon and the sun, where the subscript E indicates the earth, M the moon, and S the sun we find G(M E + M M > = n^ G(M E + M M ) = n^(l - e M cosE M )" 3 (9,6) with n„, e M , E M being elements of the moon's orbit where n = mean motion, e = eccentricity, E = eccentric anomaly. r„ stands for the radius of the moon's orbit. Hence, GM M V 1 = M— +V " M (1 - e M cosE M )- 3 (9-7) r M E M 44 For the sun GM M S 2 -3 "T = M Q + M r n S (1 " e S cosE S ) ^ ; r s b h (9-8) n^, e„, Eqj Po are the corresponding quantities for the sun's orbit referred to the earth; the ratio M^/ (M„ + M^) is very nearly 1. In order to get the inertial coordinates £_,£_,£,- (eq (9-3)) for the moon and sun, respectively, 3 rotations have to be applied 5- with P 1 C-e)P 3 (-n)P 1 C-i)P 3 (-u) (cosE - e)/(l-e cosE) /l-e 2 sinE/(l-e cosE) P 1 (-e) = P 3 (-fi) = 10 cose -sine sine cose cosfi - sinfi sinfi cosfi 1 •■ 1 P ] _(-i) = cos i -sin i sin i cos i COSU) -sinaj P 3 (-oo) = sinoo COSU) 1 (9-9) 45 £ = 23° 26' 44'.'84 fi,oj, i are given for the moon by O- = 12° 06' 46V05 - 0? 0529 5 386 52T u = 196° 43' 52V316 + . 164 358 002 5T cos i M = 0. 995970322 sin i M = 0. 089683648 with T = MJD = 3328 2.0 For the sun, i„ = 0. Hence, P, (-i) = I, and P~(-ft) and P„(-co) may be combined into one transformation employing co = Q~ + U}„. Q + cu s = Sg = 282°04' 45'.'92 + 0° 0000470684T For the solution of equations (9-7), (9-8) values for e and E are also needed: e g = 0.01677194 e„ = 0. 054900489 M With the help of m = E - e sin E values for E are computed by iteration. The mean anomaly m is found by m = 358° 0' 2V42 + 0° 98 56002647T m M = 215° 31' 53V26 + 13° 0649924490T The orbital elements given here for the sun and the moon have been computed for 1950.0 from values given in Explanatory Supplement to the Astronomical Ephemeris and the American Ephemeris and Nautical Almanac [I960]. Only secular terms have been considered here, i.e., the orbits of the sun and moon have been approximated by rotating ellipses. 46 9.2 Direct Solar Radiation Pressure The force exerted by sunlight pressure is proportional to the solar flux and to the area of the satellite as pro- jected on a plane perpendicular to the direction of the flux. According to Shapiro [1963] the acceleration is K A r Q = Cjcr-) (b — (9-10) x sp M c r r sa S where r q is the vector from the earth to the sun, as in section 9.1 , I is known as the solar flux in the vicinity of the earth, c is the speed of light, (A/M ) is the area-to-mass ratio sa of the satellite, and K R is a function lC(r,r«) whose value is or 2 depending on the location of the satellite in its orbit. In the space cylinder formed by the earth's shadow there is a sharp cutoff of the force for any point within the cylinder. Here the value K = is used. For the other part of the orbit the value 2, which expresses the reflection characteristics of the satellite, is used. Since the earth's orbit is elliptical, the solar flux will vary during the year: a 9 I = I (— ) (9-11) E 47 where a is the semimajor axis of the earth's orbit, and r_ e L is the earth-sun distance. The solar constant, I , which represents the mean rate at which enery is received at a point 2 on the earth, is approximately 2 cal/cm sec. and remains quite constant. Because the arcs used are not longer than 7 days equation (9-11) was not applied. To decide whether the satellite is in sunlight or in shadow the inner product of the position vectors of the sat- ellite and the sun is computed. £ is the angle between these two position vectors. The satellite is in shadow, if the following conditions are fulfilled (see also fig. 3) 1) cos L - J of the earth's omputation of f f ects 48 The components of eq (9-10) are transformed to the £-. 9 £ 9 j£;_ coordinate system as defined in section 9.1, i.e., the same rotations are used as for the computations of the attraction 3f of the sun. The partial K T as necessary for eq (2-5) is numerically integrated in the same manner as the variational equations . 9.3 Atmospheric Drag Most discussions of air drag in the current literature on artificial satellites begin with the equation for the magnitude of the drag force (e.g., Shapiro 1963 ) ■ ,2 with Q nd = " 2 C D P C R A •) sa D (~ — ) = area-to-mass ratio of the satellite sa C~ = dimensionless drag coefficient here approximately 2 and a variable in the adjustment p = air density |r n | = relative velocity of the satellite with D respect to the atmosphere " x + coy r* D i •*• r D = y - cox = y D z -1 Lz D J where co is the average angular velocity of the earth and |? D | = /(*g + y£ I z^ (9-17) C9-18) 49 The drag resistance is then i nd 2 C D p ( R— } ' r D sa (9-19) The air density p is given as a function decreasing' with height p = p .exp [-(h - 8-10 5 )Ap] , h is the height of the satellite in meters obtained from the following formula (9-20) h = r 1 - R . ma 3 Ar 2 + E z 2 ) s (9-21) with F . = semimaior axis ma j J 2 2 R . - R . and E = _I3^L SiH , R 2 . mm R . = minor axis , mm ' p in equation (9-20) is a value for the air density at a height of 80 0km and is taken from the U.S. Standard Atmosphere Supplements, 1966. 5 The value 8-10 is incorporated because of the chosen p and the unit of meters for h. A is the inverse of a quantity known as the "scale height" and is computed by Ap = In (-1)/(2.10 5 ) Po (9-22) 50 wh ere p, is a value for the air density at a height of 1000km and is also taken from the U. S. Standard Atmosphere Supple- 5 ments. The denominator 2*10 gives the height difference between p and p, in meters. In order to get the values for p and p, the exospheric temperature for the observation time has to be determined. For this purpose the solar flux, the monthly solar flux, and the geomagnetic index A for this period are taken from the World Data Center at Boulder, Colorado. The formulas used here can be found in the U. S. Standard Atmosphere Supplements. First the variation in the solar cycle is taken into account T = 362 + 3.60 F 1Q ? (9-23) — . -22 where F, n 7 is the 10.7cm solar flux in units of 10 2 watts/m /cycle/sec averaged over three solar rotations. The variation within one solar rotation yields T' = T + 1.8 F 10 _ ? - F 10 _ 7 ) C9-24) where F, n 7 is the daily solar flux averaged over the ob- servation time of one week. Now the semiannual variation is supplied T = TJ + f(d) F Q (9-25) 51 where f(d) = (0.37 + 0.14 sin(2^ d 3 6 5 51 )) sin Utt^^-) with d = number of days elapsed since January 1 of each year. The correction for the diurnal variation was also supplied T = 1.1 x T This expression is a very simplified version of the equation for the diurnal variation found in the U. S. Standard Atmosphere Supplements . Finally the variation with geomagnetic activity yields T=T+AT (9-26) 00 where AT = a + ICQ Cl - exp C-0.08 A 3) P P The geomagnetic index A is averaged over the observation time. With the value for the exospheric temperature T determined from equation (9-26), values for the air densities p and p., are taken out of the tables 6.1, 6.2, 6.3 given in the U. S. Standard Atmosphere Supplements. Values for the air density computed with equation (9-20) are compared with values given by Jacchia [19 70] and found to be acceptably in agreement over the range of the height of a satellite orbit. 52 In order to obtain the partial tttt- in eq (2-5) the 9C D partial derivative of Q , with respect to C- is numerically integrated in the same manner as the variational equations. 1 . Adjustment of Observations Rewriting equation (2-5) and using a vector-matrix notation we get A? = A Ae + B Ay + D Ar + E Ao + F Ab (10-1) ^ s with Ae = Ax Ay Az Ax Ay Az *■ H k = 1,6 ik Ae being a (6x1) vector and A a matrix of Cg x 6) elements, where g indicates the total number of observation equations per arc . ^Axi I = 1,104 AX A ^ AX 1 o «t J B = i i£ "The program provides an option to do the adjustment with the parameters for air drag and radiation pressure CCp. and K R ) as constants. In that case equation (10-1) is solved without Ao. 53 with B being a matrix of (g x 104) elements Ar = s Ar si 1 Ar Sl 2 Ar Sl 3 Ar s 2 1 Ar Ar Ar ni n2 n3 Ax si Ay sl Az s 1 Ax S2 Ax Ay< Az n n n D = 9f ia a = l,3n D being a matrix of (g x 3n) elements with a=3(q-l)+m where, as explained on page 5, n stands for the total number of stations, m denotes the mth coordinate of r , and q indicates the qth station. Note that this matrix has blocks of zeros. Ao = AK R AC D E = V3K R ) > \3C D ) with E a (g x 2) matrix Ab = 4b x • • Ab. • 1 Ab, *■ n j = l,t ID Ab being a vector of t elements, where t stands for the total number of passes per arc and Fa (g x t) matrix. Note that F 54 has blocks of zeros. Equation (lQ-1) may be written as [F, G, H] Ab Av = t + v C10-2) Here I is a vector equivalent to Af in equation (10-1) and v the vector of the residuals where G = [A,E] , H = [D,B] Av = Ao , and Ari = Ar By means of the covariance matrix S. associated with the I observations we get the normal equations Ab T -1 T -1 T -1 F Z £ X F F Z £ G F E £ H T -1 T -1 T -1 G Z £ F G E £ G G E £ ■ L H T -1 T -1 T -1 HE. F H E £ X G H E £ H Av An T -1 f e a T -1 G h * T -1 H E £ I (10-3) In the formation and solution of the least squares normal equations, the unknowns fall into three groups. The three groups are: bias parameters (one Ab for each pass), orbital parameters (orbital elements e, plus air drag C^, and radiation pressure iC,), and surface parameters (station K coordinates r and density values x)« 55 The normal equations formed from all observation aquations having the form of (2-5) contain all three groups of unknowns; a bias parameter for every pass, a set of orbital parameters for every orbit, and one set of station coordinates and densities. Rather than directly form and solve these large normal equations, we eliminate the bias and orbital parameters at the earliest possible stage. Only the reduced normal equations associated with the station coordinates and densities, the parameters of main interest, are accumulated over all orbits. The elimination of the bias and orbital parameters is accomplished by an application of Gaussian elimination that takes account of the special structure of the normal equations. Generally, the elimination of X, fro: >m M T P Q ~ x l~ 3 . u. u 2 J (10-4) 56 produces reduced normals for X 9 of the form (Q - P T M _1 P) X 2 = U 2 - P T M~ 1 U 1 . (10-5) If desired, X, can be recovered from a "back solution" X ] _ = M 1 (U 1 - P X 2 ) . (10-6) If M is "block diagonal", that is, of the form M = m 11 m 22 m 33 m nn (10-7) with a corresponding partitioning of P and U, , P = m m m lc 2c 3c m nc U. u- u, u n the reduced normals become; n Q - E m. m. . m. j = 1 DC ]] ]c n X, T -1 U - l m. . m. . u. 2 . =1 11 :: ] (10-8) Clearly, each term of the summation can be accumulated as soon as m. • , m. ,and u. are available; the terms can be accumulated in any order. 57 T -1 In equation (10-3), F Z F, the submatrix of the normal ■*• T -1 equation associated with the bias parameters Ab, and G Z G, •*■ T I the submatrix of the normal equations associated with the orbital parameters Av, has a block diagonal form which is retained after the elimination of Ab. The scheme described above is applied in succession to eliminate first the bias parameters and then the orbital parameters. The contribution to the reduced normals arising from the elimination of the jth bias parameter is computed during the formation of the partial normals arising from the jth pass. After all observations involving the 1st orbit have been processed, the accumulated reduced normal equations (reduced by the elimination of all bias parameters) are solved for pre- liminary values of corrections to the orbital parameters of the 1st orbit (orbital elements, air drag and radiation pressure) and surface parameters (station coordinates and densities). The normal equations are then reduced a second time to eliminate the orbital parameters. These reduced normal equations, along with the solution for the densities and station coordinates, are stored on tape for introduction during the formation of normals from the second orbit according to formula (10-8). During the processing of the observations from the oth orbit, the reduced partial normals from this orbit (bias parameters eliminated) are accumulated, added to the fully reduced normals carried forward from the (o-l)th orbit, and solved for the 58 oth preliminary values of orbital (for the oth orbit) and surface parameters. Then the oth orbital parameters are eliminated to produce reduced normals to be carried forward, along with the oth preliminary values, to the (0+l)st solution. The final solution for station coordinates and density values is obtained after all orbits are processed. In fact, the accumulation involved in the formation of reduced normals, as above, is already an intrinsic part of most Gaussian elimination programs, such as the Gauss-Jordan used here. Retrieval of the reduced normals from such a program requires only minor modifications. The results of the adjustment for the density values x 5 are now taken to compute normalized coefficients C and S r nm nm in order to compare these values with existing satellite results. If C and S denote the harmonic coefficients nmu nmu introduced in eq (5-2) to define the normal potential U, we obtain [Koch 1968] n 104 C = C + - £ y \ \ n — nm nmu fn , . s, „ n n , A Z {* r P (sine})) cos mA dE (10-9) (2n+l)kM a £=1 AE p nm -, 104 S = S + I Xn I! r n P (sine})) sin mX dE (10-10) nmu (2n+1)kM a n ^H ^ The integral over the surface elements AR in ((10-9) and (10-10)) is solved numerically by dividing AE, into 9 subdivisions. The origin of the earth-fixed coordinate system and of the orbital system--as mentioned earlier--is the center of 59 mass of the earth. Hence the harmonic coefficients C, n , C, -. , and S _ , must equal zero. Furthermore Co-i and S^-, should be small in comparison to the rest of the harmonic coefficients since during the orbit computations the z-axis effectively coincides with the rotational axis of the earth. To insure this, constraints in the form of observation equations with small variances (see Koch and Pope 1969) are set up ac- cording to (eq (10-9) and (10-10)) and their contribution is added to the normal equations. Five such constraint equations are used, setting each of C, n , C, ,-j S-, , , C„, , and S~-. equal to zero. As an option, C can also be constrained. r ' oo In addition to these constraints the longitude of one station was held fixed in the solution to prevent a singularity which would arise if all the station longitudes and the right ascension of each orbit were all unknown. We get the follow- ing condition equation j- y + Ay . y arc tan - — t — tt~ - arc tan 2 - x + Ax x or by linearization &■ - 2$& ) = o . (lo-ii) 1 ♦ C^) 2 V x x 2 / X v For some of the stations, which have been moved during the observation time of different orbits, the surveyed distances between the old and new location are introduced as constraints. 60 The differences in the coordinates (Au,Av,Aw) are weighted based on a mean square root error according to the obtained accuracy and observation equations with these differences are added. — _ — — i u s m u s n V s m + V s n w s _ m w s L- n_ Au mn Av mn Aw L mn. (10-12) Acknowledgment The author wishes to thank the director of the NOS Geodetic Research and Development Laboratory, Dr. Hellmut H. Schmid, and the members of this organization for their generous help, especially Mr. Bernard H. Chovitz, Mr. Lewis J. Gulick, Dr. Karl R. Koch, Mr. Foster Morrison, and Mr. Allen J. Pope. 61 REFERENCES Anderle, R. J. , "Geodetic Parameter Set NWL-5E-6 Based on Doppler Satellite Observations," NWL Report No. 19 78 , Dahlgren, Va. , 1965. Bean, B. R., and Dutton, E. J., "Radio Meteorology," National Bureau of Standards Monograph 92 , U.S. Govern - ment Printing Office, Washington, D.C., 19 66. Brouwer, D., and Cleraence, G. M. , Methods of Celestial Mechanics , Academic Press, New York and London, 19 61. Danby , J.M.A., Fundamentals of Celestial Mechanics , The Macmillan Co., New York, 19 62. Explanatory Supplement to the Astronomic Ephemeris and American Almanac, Her Majesty's Stationery Office, London, 19 61. Fricke, W. , et al . , "Report to the Executive Committee of the Working Group on the System of Astronomical Constants," Bulletin Geodesique IMo. 75, p. 59-67, Paris, 1965 Gaposchkin, E. M. , and Lambeck, K. , "1969 Smithsonian Standard Earth (II)," Smithsonian Astrophysical Observa - tory Special Report 315 , Cambridge, Mass., 1970. Gulick, L. J., "A Comparison of Methods for Computing Gravitational Potential Derivatives," ESSA Technical Report CSGS M- , Rockville, Md., 1970. Heiskanen, W. A., and Moritz, H. , Physical Geodesy , W. H. Freeman, San Francisco, 1967. Henrici, P., Discrete Variable Methods in Ordinary Differential Equations , John Wiley and Sons , New York, 1962. Henrici, P., Elements of Numerical Analysis , John Wiley and Sons, New York, 19 64. Hopfield, H. S., "The Effect of Tropospheric Refraction on the Doppler Shift of a Satellite Signal," Journal of Geophysical Research 68, pp. 5157-68, 1963. 62 Hopfield, H. S., "Two-Quartic Tropospheric Refractivity Profile for Correcting Satellite Data," Journal of Geophysical Research 74, pp. 4487-99, 1969. Jacchia, L. G., "New Static Models of the Thermosphere and Exosphere With Empirical Temperature Profiles," Smithsonian Astrophysical Observatory Special Report 313 , Cambridge, Mass., 19 70. Kaula, W. M. , Lee, H. K. , Taylor, P. T., and Lee, H. S., "Orbital Perturbations From Terrestrial Gravity Data," Final Report, Contract AFC 601 ) -4171 , 213 pp., Institute of Geophysics and Planetary Physics , University of California, Los Angeles, 1966. Koch, K. R., "Alternate Representation of the Earth's Gravitational Field for Satellite Geodesy," Bollettino Geofisica Teorica ed Applicata 10, pp. 318-325, Trieste, 1968. Koch, K. R., and Morrison, F. , "A Simple Layer Model of the Geopotential From a Combination of Satellite and Gravity Data," Journal of Geophysical Research 75, pp. 1483-92, 1970. Koch, K. R., and Pope, A. J., "Least Squares Adjustment with Zero Variances," Zeitschrift fur Vermessungswesen 9 4 pp. 390-397, Stuttgart 1969. Koch, K. R., and Witte , B., The Earth's Gravity Field Represented by a Simple Layer Potential from Doppler Tracking of Satellites, presented to the Fall AGU Meeting in San Francisco, Cal . , December 1970. Riley, J. D., Bennett, M. M. , and McCormick, E. , "Numerical Integration of Variational Equations," Mathematics of Computations 21, 12-17, 19 67. ' ~ Romanelli, M. J., "Runge-Kutta Methods for the Solution of Ordinary Differential Equations," Mathematical Methods for Digital Computers , Ralston, A., and Wilf, H. S., eds., John Wiley and Sons, New York, 19 60. Schmid, H. H. , "Application of Photogrammetry to Three- Dimensional Geodesy," EOS, Transactions, American Geophysical Union 50, pp. 4-12, 19 69. 63 Shapiro, I. I., "The Prediction of Satellite Orbits," Dynamics of Satellites , Roy, Maurice, ed., Springer Verlag, New York, 1963. Smithsonian Institution, "Geodetic Parameters for a 1966 Smithsonian Institution Standard Earth," Lundquist, C. A., and Veis , G., eds . , Smithsonian Astrophysical Observatory Special Report 200 , Cambridge, Mass., 19 66. "U.S. Standard Atmosphere Supplements, 19 66," prepared under sponsorship of Environmental Science Services Administration, National Aeronautics and Space Administra- tion and United States Air Force, Government Printing Office, Washington, D.C. Witte, B. U. , "Vergleich verschiedener tropospharischer Refraktionsmodelle fiir die Korrektion von Doppler- Messungen nach kiinstlichen Erdsatelliten , " (Comparison of Tropospheric Refraction Models for the Correction of Doppler Observation to Artificial Satellites), Allgemeine Vermessungs-Nachrichten , Karlsruhe 19 71, in press. v U. S. GOVERNMENT PRINTING OFFICE : 1971 O - 420-626 (Continued from inside front cover) C5GS 18. Submarine Physiography of the U.S. Continental Margins. G. F. Jordan, March 1962. Price $0.20. *C5GS 19. Analytic Absolute Orientation in Photogrammetry. G. C. Tewinkel, March 1962. *CSGS 20. The Earth as Viewed from a Satellite. Erwin Schmid, April 1962., *CSGS 21. Analytic Aerotriangulation. W. D. Harris, G. C. Tewinkel, C. A. Whitten, July 1962 (Corrected July 1963) *C5GS 22. Tidal Current Surveys by Photogrammetric Methods. Morton Keller, October 1963. *CSGS 23. Aerotriangulation Strip Adjustment. M. Keller and G. C. Tewinkel, August 1964. *C$GS 24. Satellite Triangulation in the Coast and Geodetic Survey, February 1965. CSGS 25. Aerotriangulation: Image Coordinate Refinement. M. Keller and G. C. Tewinkel, March 1965. Price $0.15. C5GS 26. Instrumented Telemetering Deep Sea Buoys. H. W. Straub, J. M. Arthaber, A. L. Copeland, and D. T. Theodore, June 1965. Price $0.25. *C£GS 27. Survey of the Boundary Between Arizona and California. Lansing G. Simmons, August 1965. *C§GS 28. Marine Geology of the Northeastern Gulf of Maine. R. J. Malloy and R. N. Harbison, February 1966. CSGS 29. Three-Photo Aerotriangulation. M. Keller and G. C. Tewinkel, February 1966. Price $0. 35. C5GS 30. Cable Length Determinations for Deep-Sea Oceanographic Operations. Capt. Robert C. Darling, June 1966. Price $0.10. C§GS 31. The Automatic Standard Magnetic Observatory. L. R. Alldredge and I. Saldukas, June 1966. Price $0.25. ESSA Technical Reports--C$GS C$GS 32. Space Resection in Photogrammetry. M. Keller and G. C. Tewinkel, September 1966. Price $0.15. CSGS 33. The Tsunami of March 28, 1964, as Recorded at Tide Stations. M. G. Spaeth and S. C. Berkman, July 1967. Price $0.50. CSGS 34. Aerotriangulation: Transformation of Surveying and Mapping Coordinate Systems, Lt. Cdr. Melvin J. Umbach, July 1967. Price $0.25. C5GS 35. Block Analytic Aerotriangulation. M. Keller and G. C. Tewinkel, November 1967. Price $0.55. CSGS 36. Geodetic and Grid Angles — State Coordinate Systems. Lansing G. Simmons, January 1968. Price $0.10. CSGS 37. Precise Echo Sounding in Deep Water. G. A. Maul, January 1969. Price $0.25. CSGS 38. Grid Values of Total Magnetic Intensity IGRF--1965. E. B. Fabiano and N. W. Peddie, April 1969. Price $0.60. CSGS 39. An Advantageous, Alternative Parameterization of Rotations for Analytical Photogrammetry. Allen Pope, April 1970. Price $0.30. C5GS 40. A Comparison of Methods of Computing Gravitational Potential Derivatives. L. J. Gulick, September 1970. Price $0.40. C£GS 41. A User's Guide to a Computer Program for Harmonic Analysis of Data at Tidal Frequencies. R. E. Dennis and E. E. Long, April 1971. PENN STATE UNIVERSITY LIBRARIES lilllllllllllll ADD00720no^2