A Cowell-based semi-analytical procedure for generating orbits of low eccentricity

A Cowell-based semi-analytical procedure for generating orbits of low eccentricity

Adv. Space Res. Vol. 6. No. 9, pp. 119—133. 1986 Printed in Great Britain. All rights reserved. 027~1177:86S0.~+ .50 Copyright © COSPAR A COWELL-BAS...

1MB Sizes 0 Downloads 33 Views

Adv. Space Res. Vol. 6. No. 9, pp. 119—133. 1986 Printed in Great Britain. All rights reserved.

027~1177:86S0.~+ .50 Copyright © COSPAR

A COWELL-BASED SEMIANALYTICAL PROCEDURE FOR GENERATING ORBITS OF LOW ECCENTRICITY R. H. Gooding Royal Aircraft Establishment, Farnborough, Hampshire, GUJ46TD, U.K.

ABSTRACT

A new approach to the semi—analytical modelling of Earth—satellite motion is proposed. It is based on the standard Cowell procedure in which the satellite’s equations of motions are integrated numerically, but it uses analytical formulae for the high—frequency perturbations associated with geopotential harmonics of high degree; this avoids the need for very short integration steps. As the integration procedure dominates in the new approach, it differs fundamentally from the author’s earlier method in which all the short—period perturbations are analytically represented and a minimal Euler integration operates on a set of mean orbital elements. Thus much less mathematical development is now necessary and it should be possible to produce operational software in a relatively short timescale. I NTRODUCT ION The author has previously /1/ proposed a mixed procedure for orbit modelling in which almost all the perturbations are represented analytically, via formulae for position and velocity in a system of polar coordinates defined by a mean orbital plane. The coordinate system was originally /1,2/ chosen to be cylindrical polar, but for a complete (e—untruncaced) second—order theory for perturbations due to .J 2 and J3 /3/ it was found that the associated spherical—polar system is better. Full details of this theory are still to be published /4/, but /3/ covers the essentials. The trouble with the analytically dominated approach is that an immense amount of work has still to be done. It is easy to cover the zonal harmonics, using recurrence relations, but harder to cover the tesserals, though general formulae for an orbit of low eccentricity have already been given /1,2/. The present coverage of lunisolar perturbations is quite inadequate, however, and other perturbation sources have yet to be considered at all, so the steps that remain (to achieve an accuracy of, say, 1 m in position) are still as summarized in /5/. Much more rapid progress might be made if the existing approach could be reversed, so that sources of high—frequency variation being treated analytically. The object almost all 1ythetheperturbations would be covered by numerical integration, using Cowell’s here is to method, on avoid the inefficiency of the very short step—size that becomes necessary, in a pure Cowell procedure, when the accuracy requirements are such that the geopotential field must be taken to high degree. Thus it has been found /6/ that the length of a step (a fixed step—size being appropriate for the low—e orbits we are assuming) should not be greater than about 10% of the period of the perturbation of highest frequency present; for a 180 x 180 geopotential field this would lead to 1800 steps per orbital revolution, though 40 would suffice for ltn accuracy in an unperturbed orbit. A cubic law of inefficiency is indicated here, since if the field is increased from (say) 36 x 36 to 180 x 180, one factor of 5 arises from the increased (maximum) degree, another from the increased order and a third from the reduced step—size; thus (and as a rough guide only) the result is an increase in program run time by a factor of more than 100. The premise of the present paper is that there is a r6le for a Cowell—based semi—analytical orbit generator, for the reasons indicated. An obvious application would be to an orbit such as that proposed for ERS—1, since it has been found by Wakker et al /7/ that a field of degree up to 50 would be needed to provide a radial accuracy of 20 cm for this satellite. Since high—frequency perturbations are eliminated from the Cowell procedure of the new generator, we may regard this procedure as generating ‘filtered’ position and velocity, 119

120

r

R. H. Gooding

and

r

say, where the true position and velocity are related to these by r

=

r+Sr

and

~

Here 6~ and d~ are computed by the analytical component of the generator and we can suppose that cSr is normally quite small (less than 1 1cm) — we cannot regard 6? particularly small, precisely because 6? is the derivative of a quantity varying at high frequency, but this is of little consequence. The significance of the smallness of 6,~ is that a first—order analysis is adequate, so that we can neglect O(Jp.m) effects associated with any analytical perturbations due to the harmonic ~im . We also truncate the analysis on the hypothesis that O(eJim) terms are negligible, making the assumption (correct for the proposed ERS—1 orbit) that e is of the same order of magnitude (10 ) as Though J 2 effects will always be fully covered by acceleration terms in the Cowell component, the variation (due to J2 ) of quantities required in the analytical formulae will be neglected (or only partially covered at best), so for any ~tm represented by the analytical component there will be neglected 0(J2JLm) terms as well as the O(eJim) terms. In the next section we introduce the spherical—polar system, based on a smoothed orbital plane that is instantaneously defined by the filtered position and velocity. All the relevant formulae relating perturbations, expressed in this system, to the (rectangular) components of S~ and ~ are given, after which the subsequent section presents the actual perturbation expressions that arise from a given ~tm Here, and in the following section which describes a pilot computer study, the particular difficulties associated with orbital resonance are indicated. The conclusion of this introductory paper is that the new approach is promising, but that further study is needed, in particular of the treatment of resonance. THE SMOOTHED ORBITAL PLANE AND SPHERICAL—POLAR COORDINATES

The author has previously /1,2,3/ demonstrated the advantages of working with analytical perturbations in terms of polar coordinates, defined by a mean orbital plane, rather than sticking with perturbations in standa~dorbital elements (a, e, i, 17, w and M). The mean orbital plane, and the position of the ‘mean satellite’ lying within it, can only be specified implicitly, in terms of the mean orbital elements that are fundamental to the theory. Thus if (X, Y, Z) are rectangular coordinates relative to the mean orbital plane, with (x, y, z) the normal equatorial—based coordinates, we have

(x

y

z)T

R3(—17) R1(—i) (X

Y

Z)

,

(1)

where Rj( ) describes a coordinate rotation about the j’th axis and T denotes matrix transposition. Also /3/ the (spherical) polar coordinates (r, b, w) are related to (X, Y, Z) by X

rcosbcosw

y

rcosbsinw

Z

rsinb

,

(2)

so that w is the in—plane quasi—longitudinal angle and b is the out—of—plane quasi— latitudinal angle. Since the mean satellite lies in the mean orbital plane, b is identically zero; also ~‘ 17, where u ( = W + v, v being true anomaly) is argument of latitude. Since b ( b + Sb) remains small, no polar singularity can arise in the polar coordinate system. In the present work the concept of the mean orbital plane becomes that of a ‘smoothed orbital plane’ which is implicitly defined by the Cowell integration and relative to which the analytical perturbations are formed. Thus the smoothed orbital plane is the plane of the osculating orbit for the_filtered position and velocity, and a complete set of smoothed elements (~,~, i, 17, ~ and M) can be defined on this basis. It follows from the way in which smoothed elements are introduced that they are actually much less ‘smooth’ than mean elements and that they are much less important in the new theory than mean elements are in the other theory. It may be useful to register their values at regular intervals (epochs of orbit determination), but only on the basis that they are equivalent to filtered position and velocity components. Further, their precise meaning will vary from satellite to satellite.

Semi-Analytical Orbit Generation

121

Though we follow the older approach, in computing the appropriate d~ and 6? in the (r, b, w) coordinate system, we do not apply them as directly as before (Sr added to ? etc, followed by the application of equation (1) and then equation (2), with full details still to be published /4/), because the integration process yields (5~, ~, ~) and not ~‘, b,~). The formulae for (dx, dy, dz) and (6?, di’, 6?) from (dr, Sb, Sw) and,,~ (6?, db,_S~)will now be developed, therefore. For convenience we Set C Cos i I ~sin i ,_C cos 17 and S = sintl , and also write W for the matrix product R~ 3(—17) R1(i) . Then C

~S

—~S

S

~C

0

1

—~C

From the analogue of equation (1) we have (since

(dx

dy

Sz)T

=

(3)

.

i

and

14 (SX

ISY

12

are fixed)

ISZ)T

(4)

where, following equation (2),

dx

cosbcos~ Sr

SY

cos b sin ~

SZ

sin b

+ ~—sinbcos~db —

sin

+

~—cos~sin~

sin ~

Sw

.

(5)

cos b cos

cos b

0

Thus expressions for dx, Sy and Sz are obtained by pre—multiplying the three terms of equation (5) by W and setting b = 0 and ~ = as we proceed, where, by inversion of equations (1) and (3), -

~r _i rj cos U

r~ C

=

sin’17

—~

0 (x

S

y

z)

(6)

.

I

~

~)‘r

The first term of equation (5) gives, trivially, (dr/~)(~ ~ (2 + =2 + -=2)½ . The second term is just ~I(O o l)’~ Sb

,

ie

where, of course, ~ ~ Sb

=

~(Is

Finally, the third term may be written ‘~

—~S cos

Sw

which, by equation (6), reduces to

C



(—~ — IC~

~

sin ~

— ~

1C~+ IS~)T Sw

Thus the developed formula from equation (4) is dx

=

(dr/~)

I

+ ~

IS

Sb

+

-~-IC~

Sw

.

(7)

-~ I

ICI+~S7

It will be noted that C and S only occur with the multiplying factor 1 , this being an inevitable consequence of the fact that dx, dy and 6? must still be determinate when I is zero so that 12 is indeterminate (equatorial orbit). For Si, ~ order,

JASR6~9-I

and

51

,

the starting point is to observe that equation (1) gives, to first

(x

y

z)

=

W (X

Y

Z)

,

(8)

122

R. H. Gooding

since this result would be exact (b~the whole concept of osculation) if W were defined in terms of i and 12 rather than i and 12 . Thus the equation (to first order) corresponding to equation (2) becomes (paralleling equation (5) of course),

x

I

cos b cos w

Y

cos b sin w

Z

sib

+

rb



sin b cog w

+

r~’— cog b sin w



sin b sin w

cos b cos w

cosb

0

.

(9)

We now get (dx SY SZ)T as the sum of nine terms, three from each of the three terms of equation (9). Thus the first term of equation (9) gives

6?

cosbcosl

+

~Sb

—sin~cog1

cos b sin 1

+

~dw—cosbsin1

sin b sin I



5mb

cos b cos

cosb

I

0

the second term gives (~db+bSr)—sin~cos1

—I~Sbcos~cosI

—sinbsinl

+~bSw

cosbslnl

sinbsinl

cosb

—sinbcosl

sib

0

and the third term gives

(~dsn

+ 1

6?)



cos

sin 1

‘~

cos b cos

+

~

Sb

‘I

sin b sin 1 —

sin

cos

0



~

Sw cos

cos I

cos

sin 1

I

0

0

Of these nine terms, the fifth, sixth and eighth vanish, both

b and b

being zero;

the first and last terms combine and reduce, like the first term of equation (5), to —

1

l)~

Sw)(~ ~

;

terms combine and reduce, like the

the second and fourth

second term of equation (5), to (~Sb + I Sb)(IS —C ~)T finally, the third and seventh terms combine and reduce, like the third term of equation (5), to

[s~

+ (~Sw +

1

6r)/~](—If — 55c d~

together constitute (

IsI



6?)T

,

(11

+ yy + zz)/r

,

and

w

~~)T

.

The three reduced combinations

1

corresponding to equation (7), ~

-

course by

I~ +

by

being given of

-=2

n(a/r) (1



e

EXPLICIT ANALYTICAL PERTURBATIONS IN SPHERICAL—POLAR COORDINATES

General formulae for perturbations due to the geopotential It was shown previously /2/, via the integration of Lagrange’s planetary equations, how general formulae for perturbations due to the geopotential field could be expressed in cylindrical—polar coordinates (r’, u’, c’), defined such that, corresponding to equation (2) here,

(X

Y

Z)T



(r’ cos u’

r’ sin u’

c~)T

We summarize the development of /2/, giving the results now in spherical—polar coordinates, with expressions for velocity components (omitted before) included. We take the general term of the Earth’s gravitational potential as notation,

Ut

=

E.(!)

P~(sinB)(C~ cos mA + St

ULm

sin mA)

where, in the usual

.

(10)

Semi-Analytical Orbit Generation

For convenience

we introduce

‘~fm and

Aim

123

the usual polar equivalents of the harmonic

,

coefficients, defined by (Ci,

5Lm~

(cos mX~, sin mAt)



(11)

and also the normalized equivalents such that (Cm, 5Lm’ ~tm~ where for

m ~ 0

Nj(Ct,



the standard value of

Nm

(12)

is given by

2(2t + 1)(L



atm’ ~tm~

m)!/(t



+

m)!

;

(13)

for m — 0 the factor 2 must be omitted from equation (13), but it Is alternatively possible to identify ~i o with the familiar unnormalized ‘zonal harmonic’ ~3t by taking Ni,O — —1 . Using equadons (11) and (12) we can rewrite equation (10) as

Ut



(~)t ~tm Ni

P7(sin B) cog m(X



~‘tm~

(14)

To express ULm as a function of the orbital elements (with u retained), as required for the integration of the planetary equations, we invoke the expansion

P~(sin B) exp

~Ftmp(i)



exp

2p)u + m(Q





v)}

(15)

of Allan /8/, where v is the sidereal angle of the Greenwich meridian, and where the inclination function Fimp is real when 9. and m have the same parity and otherwise is (pure) imaginary. However, it is more symmetrical and convenient /2/, following /9/, to replace p , as an index, by k , where k

£



2p



(16)

.

Thus we write

F~m(i)

Ftmp(i)



(17)

where equation (15) shows that the index k runs from —t to +t ; only the F~m for which k has the same parity as 2. are actually needed, but the others are useful in recurrence relations /9/. Further, we may normalize

F~m , whilst at the same time making it always real, by defining

m(i)



(/ 1)k~

Ni

F~(i)

.

(18)

As an aid to compactness in the perturbation formulae, we change the notation yet again, writing 9

~

(19)

~tm1’~

and

Si/Si



Finally, and to the same end, we introduce

-r

i and —

(20)

.

x

½t—u

by

(21)

124

R. H. Gooding

and 1~It

x Then it follows posed into

from equation (15) etc that

~ ~m

U~m to

(22)

.

Uim

,

as given by equation (14), may be decom-

where

the planetary

O(eJtm) to get

~

~.(~)L

-

To integrate

12 + V + Aim





cos

(k’r

+ mx)

(23)

.

equations, some terms of which contain

0(Jtm) results;

1 + (9. + 1) e cog v + O(e



e1

factors,

we need

but this raises no difficulty, since 2

)

Had we required results to 0(eJtm) or beyond, however, it would have been necessary to deal with this expansion more thoroughly, the standard procedure /8,10/ being to introduce a further summation index (q), together with the Hansen—based eccentricity functions, G~(e) or (equivalently) G~q(e) . As we are able to avoid this complication here, the integration can now proceed /2/ on the basis that, in particular, cos (icr + mx) dt can be set to sin (ky + mx) /(Ic~ + ii~) . This is valid on two assumptions: first, that the quantities I and x are changing at a uniform rate, which is a legitimate hypothesis in a first—order analysis, no distinction being necessary between, for example, ‘r and -r in perturbation formulae; and secondly, that it is legitimate to express results in terms of indefinite integrals rather than definite integrals — there will be further consideration

f

of this

point

later.

On the basis of the first assumption we could take I ( —Cs) to be —~ , where ~ is the smoothed mean motion, but it is natural to allow for the main secular variation of s , ie for

17

as calculated

for

J 2 alone; so we write

I

—(1 +~)



Similarly, we do not just take

~

~



(say)

(24)

.

to be ~3 , the Earth’s angular velocity n~12

nE

,

but write

(25)

.

It is convenient to have a symbol for the ratio of these two angular rates, so (still following

/2/)

we write

s

(n%





17)/n

;

(26)

clearly s is the ratio of the. Earth’s rate of rotation (relative to the satellite’s orbital plane) to the satellite’s rate of (nodal) revolution, and as such it is the fundamental parameter of resonance. As a final piece of shorthand we define a quantity c for each pair of indices, k and m , by writing c

k—ms



(27)

.

For brevity we omit the formulae for perturbations in the elements, given by equations 84—91 of /2/. Then the perturbations in the spherical coordinates are given by linear combinations of these, viz a(Se cos v + e S4, sin v)

Sr



Sa

Sb



Si sin

Sw



2(Se sin v



u



512 sin

(28)

,

I cos u

(29)

and —

e 6$ cos v) + SL

.

(30)

125

Semi-Analytical Orbit Generation

Here

5$

and SL

are non—singular perturbations 64,

Sw

defined by +

SQ cos i

(31)

and SL

Substitution

SM



of the element—perturbation Sr



ai’[{2k

Sb



(1/2c)

[{~‘ {~‘





c2)] cos (kI + mx)

cosec i (k cos I

+

cosec I (k cos i



(32)

.

e



formulae into equations (28)—(3O) gives

c(9. + 1)}/c(l



2) 5$

~f(1

+

m)} cos {(k



(33)

,



1)y +

m)} cog {(k + 1)y +



mx}

mx}

(34)

and Sw



~[{3k

2c(i + 1) + kc2}/c2(1



c2)] sin (k~’ + mx)



(35)

.

There are obvious problems if the coefficient of any term becomes large in spite of the small factors F and F’ . No difficulty arises with the cosec i factors in equation (34), since F(k cos I — m) turns Out to be always 0(sin i) . There is a real difficulty when c is close to 0 or ±1 , however, which can happen in the following ways. If m — 0 , then c = k and so k — 0 or k — ±1 ; we then have secular (k — 0) or long—period (k — ±1) perturbations, and here we shall find that the difficulty is easy to cope with. If m 0 on the other hand, the zero denominator is associated with resonance, which we shall find harder to deal with. We also require the formulae for

6?

,

‘Sb and S~r .

Ict +.n~(

Noting that

—cn



we differentiate equations (33)—(35), getting 6?



a n

~[(2k

(n/2c)[(c

Sb

c(t + 1)}/(1





1){~’ +

c2)] sin (Icy + mx)



cosec I (k cos i



m )} sin {(k

(c + i){F’ —F cosec I (k cos i



(36)

,





1)1 +

mxl

m)} sin {(k+ 1)y +m~}] (37)

and —



n~[{3k



2c(t + 1) + kc2}/c(1



c2)] cos (Icy + mx)

(38)

.

It is worth remarking that when positions only (not velocities) are to be generated, then equations (36)—(38) will only be required in the preliminary operation of filtering the position and velocity at epoch; if the filtered values are already available these equations will not be required at all. Special treatment possible for the zonal harmonics First, we observe that for

m



0

the inclination functions are symmetrical, in that

F~ 0(i)

Thus we can dispense with negative values of k’O.

F~0(i)

k

by doubling the expression for

U~ ~

when

This then fits in with the author’s original analysis of satellite motion in an axi— symmetric field /11/, equation 20 of which was (with minimal change of notation) 1ci A~(i)/2~~ kI} cos Icy —



.(~)i

J2. {u,,~P~(O)sin

,

(39)

126

R. H. Gooding

where

Uk



1

if

k

normalized so that



0

but 2 if

,

A~(i) — 1

if

1 ~ k ~ £

P~(O) — 0

when

£ and k

and

is even, but

9.

consistent with the fact that non—zero since

,

A~(i) is an inclination function cos I

if

U~ only arise when

have opposite parity.

£

is odd.

9. and k

Equation (39) is

have the same parity,

The relation between the

A and F

inclination functions may be seen from equations (18), (19), (23) and (39), with u.~ U~~ and U~identified and Ni 0 Set to —1 ; thus F~~(i) — —(fi)~F~ 0(i) and the conn~c— ting relation between the

F and A

F~0(i)

functions is





{P~(0) gj0kj /2k k’~

4(i)

(40)

.

It is now postulated that we want to avoid the inefficiency associated with the step—size that would be demanded if the Cowell component of the generator were required to cope with all the hi~h—frequency cos Icy terms in equation (39), say those with _k ~ 2. for some specific £ . Then the Cowell component can cover all U9. ~ with £ ~ & , but it does not follow that for all £ > £ the U2. ~ must be dealt with entirely analytically. What we would like is for the U~~ with (& )) Ic > 9. to be dealt with analytically, but for the residual U~0

(with

k ~

~)

to be covered by the integration.

useful for secular effects (k



0) and long—period effects (Ic

if the combination of the residual The desired splitting equation (14),

of

U9. ~

U~,o has the

This will be particularly —

1), but it is only possible

form of a legitimate potential.

turns out to be quite straightforward.





•_(R)t J~ P~(sin B)

We recall that, from

(41)

and note that when U9.,0 is used unsplit in a Cowell integration the appropriate acceleration terms are generated, as in /12/, from

~

taking Now

r and B

~u’~

(aU/ax, aU/ay, 3U/3z)



in equation (41) as functions of

Pi(sin 8) in equation (41) is a polynomial in sin B



,

(42)

x, y and z sin B

sin I sin u

,

where (43)

.

Hence equation (41) can be re—expressed as a polynomial in sin u , and then further re—expressed (since sin u — cog y) as a series in cos ky . This is merely the route that leads to_equation (39), of course, but it shows us how, after the series has been truncated at cos £y , the route can be reversed to produce a new polynomial in sin B . The only snag is the appearance of sin I , in the coefficients of this polynomial, which cannot be eliminated. But this is a minor snag, since each coefficient is just a polynomial in sin I , the essential point of this being that only positive powers of sin I are involved. Thus the Cowell procedure must carry a knowledge of i , but in a form to which there is no objection so long as a fixed value is used. The details of a general implementation of the zonal split for arbitrary 9. have not yet been worked out, but the use of recurrence relations would almost certainly be involved. The details for a J3 split are given in the present paper, in the description of the pilot computer study. General treatment for tesseral and sectorial harmonics From equation

(23), with

I

given by equation (21) and

f

x

by equation (22), we see that

I

the frequency associated with ~ is kCs — m(n% — 12) ‘ ~ en,,, by equations (24)—(27). We now assume, as in the previous section, that the highest frequency tolerable in the Cowell procedure is £m.~ . In view of the simplicity of the zonal split, we naturally look for a similar split for the tesseral harmonics (taking ‘tesseral’ to cover also ‘sectorial’, the harmonics with m — £). Since c — Ic — ma and ma > 0 , c j is maximum when Ic —t c then being negative. For a given value of a , therefore, (s being defined by the orbit)

Semi-Analytical Orbit Generation

127

we should like to split for each ~9.m with 9. > £ — as , except that no Cowell contribution (and hence no split) would be possible if at the same time 2. < ma — 9. ; but, since m ~ 9. , the latter would only be possible if a > 1 , which would imply a super—synchronous orbit and hence be of little interest. Thus the generalization of the situation for the zonal harmonics would in practice be that for & ~ 2. — ma we would operate with U unsplit in the Cowell procedure but that otherwise we would always want to split. The foregoing appears academic, however, as the desired split does not seem possible.

What

we should have liked to do involves the retention of a subset of the U~ , given by equation (23) as developed from equation (14), tracing the combination o~these back to an expression like equation (14) but with P~(sinB) appropriately replaced. This would involve the elimination of x as well as y , for which equation (43) and two associated relations would be available; unfortunately the elimination appears to be impossible. If the effects of the tesseral harmonics cannot be split, then we are forced to the analytical treatment of the whole of U9.m whenever 9. > & — ma . In some ways this is a simplification, but it means that, unless a self—defeatingly large value of £ is selected, a number of near—resonant terms must be treated analytically that we would prefer not to. The resonance problem We have seen that a difficulty arises with small denominators in equations (33)—(38), Ic with c close to 0 or ±1 . The source of the difficulty was our assumption that it is legitimate to integrate, in particular, cos (ky + m~) , assuming constant ? and ~ , into the indefinite integral (Ic? + ~~)l sin (Icy + mx) , which is unbounded as c approaches zero. Over a given time interval (from t 0 to t) It is quite a different matter if the indefinite integral is replaced by the appropriate definite integral, in fact by {sin (Icy + mx) — sin (ky0 + mxo)}/(lci + ni), since then the integral is bounded by

IttoI

There is a double rationale for the use of indefinite integrals in normal circumstances: first, that it obviously leads to much simpler expressions; and secondly, that it avoids the effective bias that otherwise attaches to the use, away from epoch, of a particular set of (smoothed) elements. But this second point only applies when the rate of change of ky + is great enough for this variable to complete a number of ‘revolutions’ over the period of time with which the epoch is associated. When ky + mx is confined to a small ‘arc’ of a single revolution, on the other hand, the ‘epoch bias’ is positively desirable as it reflects the true situation.

mx

Thus it is very easy, in principle, to allow for particular resonances, when Lagrange’s planetary equations are integrated, just by replacing the appropriate indefinite integrals by definite integrals. The corresponding adjustment to equations (33)—(38) is not quite so straightforward, however, because v and u in equations (28)—(30) have a rBle which is quite different from the r8le of y and x in the integrals. As an example, the way to allow for small c in equation (33) is not just by the subtraction of cos (IIO + mxo) from cos (ky + mx) . From the general formulae for perturbations in the elements /2/, we find that it is only in 6? ,• with coefficient 2aFk/c , that this subtraction is2),necessary; which (by not contribute close to zero. In hypothesis Se and 6$ , here) which isalso to Sr , the overall denominator is (1 — c

I

What we now have is an extra element of arbitrariness associated with the proposed method of orbit generation. In addition to assigning a value to £ , when running a particular computer program, we would have to lay down, by decree, the duplets (k, m) that were to be regarded as resonant in the analytical component of the program. It is normal to extract the greatest common divisor from k and m , writing /13/ (k, m) — y(a, B) ; so that it would only be necessary to specify the appropriate duplets (a, B) for all relevant (Ic, m) to be picked up. Instead of decreeing the appropriate duplets manually, it would be possible for the computer program to locate them automatically, on the basis of the proximity of s to a/B , but the location would have to be done once and for all at epoch (to avoid a change of decree) — even then there would be the danger that the decree would change (with s) from one epoch to another. It might seem that the whole resonance problem has now been disposed of, apart from the technicalities of a general implementation. This is not so, however, as we shall see in the account of the pilot study. Application to orbits of the ERS—1 type It is worth considering the implications of the foregoing for an orbit of the ERS—1 type, for which (as already remarked) a field to 50 x 50 might be necessary. Approximate values of the three main orbital parameters, as proposed, are 7153 km (a), 0.001 (e) and 98.5~(i),

128

R. H. Gooding

corresponding to a deep resonance (frozen orbit) with (B, a) value of a ( — a/B) just under 0.07.



(43, 3).

This implies a

Suppose we set 9. — 10 , so that 100 steps per orbit should suffice in the Cowell integration. If we follow Wakker et al /7/ in taking c < s as the criterion for a resonance decree, then the decree will operate for five (o, Ic) duplets, viz (14,1), (15,1), (28,2), (29,2) and (43,3), the last duplet being for the deep resonance. These resonance duplets involve a total of 64 ~im harmonics, or 32 for a field truncated at 36 x 36.

I

Wakker et al /7/ tabulate the amplitudes of the main resonance effects, using the rule of thumb of Kaula /10/ for the magnitude of each

Cia and ~ 9.m , viz 1O512.2 for both. On this basis they quote the most important harmonic, for shallow resonance, as J15,14 (with k — 1), which gives amplitudes of 3.4 m, 17.6 m and 162.0 a for Sr , a Sb and a Sw respectively, though J1414 (with k 01 gives a much larger amplitude (17.2 a) for Sr (Note that accurate values of C15,14 and S15,14 are now available /14/.) With 100 steps per orbital revolution, rather than the 500 that a 50 x 50 field would demand, a considerable saving in computer time should be possible, but it is difficult to put a figure to the saving. According to the cube law, the saving in the integration would be by the factor 125, but this is reduced by the requirements of the analytical component of the computer program. There are two major considerations here: first, the efficiency with which the analytical terms can be lumped together in an appropriate way, to offset the fact that we now have three indices (L, a and k) to be covered, not just the first two; and secondly, the number of times per orbit that, on average, it is necessary to compute these terms at all. A PILOT COMPUTER STUDY Introductory remarks The raison d’être for the Cowell—based approach to semi—analytical orbit modelling is that much less time should be required to produce operational software than with the author’s earlier approach (for which a great deal of theoretical development is still required). A fair amount of work will still be necessary for the full implementation, however, so it seemed sensible to make a preliminary assessment by means of a pilot study. The primary objective of this study was to assess the viability of the approach, with achievable accuracy the main consideration. Since all perturbations other than those due to the geopotential field would be covered by the integration anyway, the study was confined to the geopotential. (Special attention will have to be paid to the modelling of the acceleration due to solar radiation pressure, to avoid loss of accuracy in the integration as the satellite enters and leaves eclipse.) As the basic principle was to divide the modelled geopotential between the Cowell integration and the analytical representation, the division being specified by an 2.—based frequency above which all terms had to be treated analytically, it was decided that the minimally meaningful study would be for a field containing ~ ~2 2 and , with 2- set to 2. On this basis the study was made in three phases: in the first, J2,2 was set to zero so that the splitting of a zonal harmonic (J3) could be assessed with no tesseral harmonics present; in the second phase, J22 was added, all its effects being non—resonant and covered entirely analytically; in the final phase, part of the J22 perturbation was forced to act resonantly so that the magnitude of the resonance problem could be assessed. Only one orbit has been studied, viz a circular orbit with a — 6678 km (period 1½ h) and i — 41 . This bears no resemblance to the ERS—1 orbit, but the study has in any case been too primitive for this to matter. The results, in each phase of the study, were evaluated by comparison with results from pure integration, the starting point for both procedures being an appropriate set of position and velocity components for the given orbit. At regular intervals over a period of 6 h (covering four revolutions of the satellite) the differences between the two sets of results were output, chosen fortheoutput the following: iSh )( to = Sr , where the overall 2 — x2 +the y2 quantities + z2), to indicate error being in height; dp ( — indicate error in position; SV ( = ) to indicate the overall error in velocity; and the set r of quantities iSa, Se, di, 511 sin i, e Sw and 5(N + w), to indicate the interpretation of the errors in~terms of orbital elements.

I

I

I

I

Semi-Analytical Orbit Generation Effects of the split

J

The potential due to

J3

3 is given by

=

and this yields

129

.~ ~3

~.

2B)

(-~)~

sin

sin B (3



5

sin i (4



5 sin2i) sin u

the two components given by

u~

=

-~ ~3

=



~

~

and U3 3

With

9.

=

2

,

5 J ufR\ —‘ — i

8

sin 3 i sin 3u

3r\r/

this constitutes a ‘split’, with

program and U~ treated analytically. program /12/, we merely rewrite it as

3

retained in the Cowell component of the

To treat

ii iR~

(4—5

U~ in the same way as for the comparison

2— sin 1 0)sInB

the resulting acceleration is then 3

3 3 5 sin~i0)IR\ 1—1 (r9.—5sinBr)

~(4



8r where

~

-

is a unit vector along the z—axis.

The required changes to the representation of

U 3

for pure integration /12/ are then Immediate.

The analytical expressions to cover

U~ follow from equations (33)—(38), with doubled to allow for the convention that U~= U2. ~ + U~~ (k 0) . Writing 3 2. ‘ C (5/96) J3 (R/a) sIn i , for convenience, we have Sr

=

Sb and

Sw

=

3Ga sin I sin 3u

F and F’

(44)

,

—1212 cos 9. (2 cos 2u + 3)

(45)

—20

(46)

sin I cos 3u

,

with also (5?, Sb,

d’r)

=

3Gn(3a sin i cos 3u, 16 cos i sin 2u, 2 sin I sin 3u).

The results of the computation (with this split) showed that the maximal positional error, was less than I a over the whole 6h period. The height error, Sh , was dominated by an effect of the same period as the orbit itself, of amplitude growing at the rate of about 0.1 m per orbit, so that the value of dh did not exceed ½ a over 6 h ; this applied when the satellite was ‘started’ (epoch) on the equator, the amplitude growth being smaller with any other starting point. The errors are large in comparison with those reported /1/ for the almost pure analytical model (5 cm accuracy was quoted for

J3

modelling).

for most purposes when the value of

9.

is increased from 2 to at least 10.

However, they stay be small enough

Thus, allowing

for a reduction of effect by (3/11)2 , from Kaula’s rule, an accuracy of around 30 cm/d should be achievable. Since the observed effects were almost entirely due to O(J2J3) terms, as seen by their disappearance on performing runs with J2 set to zero, there is also the possibility that a simple enhancement of the analytical terms can be found that removes the bulk of these effects.

130

R. H. Gooding

Non—resonant effects of

~

2 2

Because it was not found possible to split the -~2 2 potential, it was necessary to cover J2,2 completely in the analytical part of the orbit generator, J2 being covered still entirely by integration and J3 being split. Most runs in this phase of the study were with all three harmonics present, but for a few test runs either J2 or J3 (or even both) was set to zero. To express the required formulae, we use 8 , defined as x — ½ir ( — V — 12 + A22), and u rather than x and y , in particularizing equations (33)—(38). Then the required expressions for 6?, ‘Sb and dw take the form

and

dr



R2 cos 2(8 + u) + it0 cos 28 + R2 cos 2(8

Sb



B1 sin (28 + u) + B1 sin (28

Sw



W2 sin 2(0 + u) + W0 sin 28 + W2 sin 2(8



u)

(47)

,

(48)

u)





u)

,

(49)

where for the it and W coefficients the suffix is simply the appropriate Ic from equations (33) and (35); but equation (34) contains a pair of terms, of which the first coefficient vanishes when k — —2 and the second when k — 2 , so B_1 combines the coefficient of the second term for Ic — —2 with that of the first term for k = 0 , whilst B1 combines the coefficient of the second term for k — 0 with that of the first coefficient for k — 2 . (This is not coincidental and extends to arbitrary ~L,m ~) To express the coefficients in equations (47)—(49) we write, for convenience, H

2

~ ~ 3 2,2l~a) fit\~

/T~fR\ “64 ~2,2~a)

Then R—2



aH (1 1)(2s — cos (3s++ 3)1) (a + + ~)2 1)(2s 2i

2

6aHsln



aH + cos — ~)2 (3s — — 3)1) (s (1 — 1)(2s l)(2s H sin s(s i (1 +— 1) cos 1)

(53)

B1

=

H sin s(s 1(1 + cos — 1) 1)

(54)

—2 2 H (1 I) + (4sl)(2s + 2s+ + 2(s —+cos 1)2(2s 3) 1)

(55)

W

5?, Sb and S4

,

(52)



=

6H 21 s(4s 2sin — 1)









(56)

—2 2 H (1 + cos I) (4s 2(s

To obtain

(51)

,



B—1

—2

and

(50)

-

1)2(2s





l)(2s

2s —

+

)

(57)

3)

we simply differentiate equations (47)—(49), with

8

=

sn~ and

The errors in the computed comparisons were now found to be large, and totally dominated by an along—track effect. As well as a secular build—up at about 10 m/d , which in orbit determinations would be subsumed in a fitted value for the mean motion, the along—track error contained a sinusoidal component of period about half a day (the period of 28 in fact) and amplitude about 5 m. All the error was shown (by runs with J 2 — 0) to be O(J2J2,2) and was assumed to be of the same origin as in the previous study /1/. It was decided to introduce the same additional terms as before, therefore.

Semi-Analytical Orbit Generation

131

As just intimated, it is a simple matter to eliminate a constant secular along—track variation empirically. When this was done, the maximum height error over 6 h was less than I m, and the maximum error in overall position less than 3 m. On extending to 24h coverage, the (maximum) height error increased to only about 1.5 m, but the overall positional error rose to more than 20 m. Even with an allowance for an increase In 2. to a realistic value, It would seem that the error in position cannot be held to negligible proportions (even without resonance considerations), but that the error in height, which is often all—important, probably can be. Simulation of resonance It was desired to have a first assessment of the new approach for a condition of near resonance, and in particular for c , given by equation (27), close to zero. With the only tesseral harmonic present, this would have meant the selection of a new orbit at roughly geosynchronous height. However, the perturbations would then be greatly reduced, unless the value of J 2,2 was correspondingly magnified, and if any geophysical constant was to be tampered with it seemed much to be preferred that it should be ( = V) , the Earth’s rotation rate. The test program was modified to permit different values of to be input, therefore, and runs were performed mainly with a magnification factor of 15 (very shallow resonance) but occasionally 15.9 (deep resonance). For k = 2 , which is the only value relevant to J22 resonance, the associated values of s increased to about 0.944 and 1.00064 respectively, with c equal to 0.11 and —0.0013. The necessity for the introduction, by ‘resonance decree’, of certain definite—integral expressions, to replace the corresponding indefinite integrals, has been explained. The requirement, in practice, was to avoid the occurrence of very large terms associated with R2, B1 and W2 , in equations (52), (54) and (57). On evaluating the appropriate definite integrals for Sa, Si, SQ and SL , and then substituting in equations (28)—(30)~ it is found that the required correction terms are as follows, writing Q = H(1_+ cos I)/(s — 1) for convenience. For Sr , we just need to add in the term 2aQ(1 + cos I) cos 2(80 — uj) that is essentially constant. For Sb , the term to be added is the oscillating — Q sin i sin (280 — 2uO + u). Lastly, for ‘Sw , the term to be added varies linearly with time, being —

3Q (1 + coa

2(s

1) sin 2(eo



00) +

nt cos 2(eo



u0)

There are_corresponding corrections to Sb and Sw , of course, represented by —Qm,, sin i cos (280 — 200 + u) and —3Qn.,, (1 + cos i) cos 2(8~ — , the correction 6? being zero.

uo)

to

To make the two terms of each definite integral compatible, it was necessary to supersede the normal computation of a — 1 (by subtraction) by the computation of a composite quantity, as explained for the general case. Even then, however, the results obtained were most disappointing. Unexpectedly, the error attributable to semi—major axis exceeded that attributable to all the other elements, approaching ½ km after four orbits and with even ‘Sh reaching 100 m in consequence. Just how different the situation is from the non—resonant one will be seen from the observation that the effect was found to be O(J~,2) , the O(J2J2,2) component being negligible by comparison. How can J~,2terms exceed J2J2,2 terms when one is expecting the opposite effect with a ratio of up to three orders of magnitude? A clue was offered by the fact that the along— track error appeared to have a quartic behaviour. Consider the behaviour of equation (23) to give

a, n and M

due to

U~2

F

=

H (1 + cos 1)2

and, from Lagrange’s first planetary equation, a where

O=2(u—O)

near resonance.

U~,2

—4~nasinU

We have,

for

R. H. Gooding

132

If we take 0 as constant (deep resonance), then (using equation for it gives

—3na/2a) the integral of the



2tsin0

Sri



6~n

Then a further integration yields a contribution to

SM

n



M

given by

3F(nt)2 sin 0

Now the perturbation in M goes directly into L and SL goes directly into dw via equation (30). But equation (30) is only first—order and so, in particular, neglects a term in (‘SM)2.

With

‘SM

quadratic in

t

,

this term is necessarily quartic.

of the order 3.7 x i06, for the given orbit, and with order of magnitude of (SM)2, when multiplied by exceeding 3m

a

,

n,,,t

Bit



Further, with

~

after four orbits, the

produces an effect (for maximum sin 0)

This is still two orders of magnitude off the observed behaviour, but at first sight it had seemed that the discrepancy was mitch greater still. It remains to be confirmed that the source of the error lies in the first—order nature of equations (28)—(30). If it Is confirmed, then a cure might be found from the inclusion of the main second—order terms, in (SM)2 in particular. It is at any rate worth remarking that the errors became no larger when deep resonance, rather than shallow resonance, was simulated. Two other ways of tackling the large resonance—induced errors were considered.

The first is

related to the observation that (SM)2 errors would not be introduced for SM perturbations applied directly to M itself. This would be a slight retreat from the argument for the simplicity of perturbations in sphericel—polar coordinates, but it would not necessarily matter for just a few resonance—related terms. To confirm that there should be a real advantage, we note that if the equation for sequent term in (in

F2t

.

~

is taken to

This will integrate to a term in

O(J~2~ there will be a con-

F2t2

(in

Sn) and hence a term in

SM), so one expects the error term to reduce from quartic to cubic.

The computer test—program has been modified along the lines suggested, but with a most disappointing outcome, since the results became much worse — the errors are now 0(J 2,2) rather than 0(J~,2) . Despite intense analysis, the reason for this has not yet been discovered, though it seems not just to be due to programming error. The other way to tackle the error that was referred to is by ‘rectification’ of the orbit. The meaning of this term is the same as in the integration of orbits by Encke’s method; indeed the principle is the same, since as in Encke’s method we are computing the evolution of a small deviation from an independently generated reference orbit. (The difference between Encke’s method and the approach here is that Encke’s reference orbit is Keplerian, and hence generated analytically, with the deviation produced by integration, whilst here it is the reference orbit that is generated by integration.) The problem with both methods of orbit generation lies in the loss of accuracy as the deviation vector grows in magnitude, and this can in principle be rectified by, from time to time, consolidating to position and velocity at an updated epoch, with a fresh start to the whole generating procedure. The objections to rectification are twofold. First, there is likely to be a serious loss of efficiency, since the restarting of a Cowell procedure is not a trivial matter. Secondly, the evidence so far is that the reduction in error is disappointing. (The test program has been modified to permit a rectification half—way through the period to be covered in the orbit generation, but the effects have only been partially analysed.) CONCLUSIONS A new approach to orbit generation has been outlined, the aim of which is to improve the efficiency of a Cowell—based integration procedure. Only the perturbations due to the geo— potential, and associated with a frequency above a given pre—eet level, are treated analytically, outside the Cowell procedure, but this should be enough to permit a reasonable size of step in the integration (not more than 100 steps per orbital revolution, say). The approach is at this stage only applicable to low—eccentricity orbits. The results of a pilot study indicate that, for such an orbit, there should be no difficulty in covering the higher—degree zonal harmonics as required, since the effect of these harmonics can be split

Semi-Analytical Orbit Generation

133

(between the numerical and analytical components of the orbit generation) in a natural way. The tesseral harmonics are not so easy to handle, since the splitting procedure breaks down, but serious difficulty only arises with harmonic terms that are associated with resonance. For resonance—associated terms that cannot be dealt with in the integration component of the program, the natural way to proceed is by computing terms that effectively convert certain indefinite integrals into definite integrals. This works well for short periods of time, but then along—track errors proportional to the fourth power of elapsed time build up. Three methods of reducing these errors have been considered, but further work is necessary to establish whether orbital resonance can be satisfactorily covered by the new approach. REFERENCES 1.

R.H. Gooding, Adv. Space Rca., !~83—93 (1981)

2.

R.H. Gooding, Phil. Trans. R. Soc., A, 299, 425—474 (1981)

3.

R.H. Gooding, Acta Astronautica, 10, 309—317 (1983)

4.

R.H. 000ding, RAE Technical Report (to be published)

5.

R.H. Gooding, Adv. Space Res., 5, 175—183 (1985)

6.

R.H. Gooding, RAE Technical Memorandum Space 357 (1986)

7.

K.F. Wakker, B.A.C. ~brosius, L. Asrdoom, Precise orbit determination for ERS—I (Delft University of Technology). Report under ESOC contract 5227/82/D/IM(SC) (1983)

8.

R.R. Allan, Proc. R. Soc., A, 288, 60—68 (1965)

9.

R.H. Gooding, Celest. Mech., 4, 91—98 (1971)

10.

W.M. Kaula, Theory of satellite geodesy, Blaisdell, Waltham (Mass), 1966

11.

R.H. Gooding, RAE Technical Report 66018 (1966)

12.

R.H. Merson, A.W. Odell, RAE Technical Report 75093 (1975)

13.

R.R. Allan, Planet. Space Sci.,

14.

D.C. King—Hele, D.M.C. Walker, Planet. Space Sd., 34, 183—195 (1986)

15, 53—76 (1967)