Calculating viewing times of satellite from groundstation

[Ed]

Hi folks,
If I have the elements of a spacecraft 's (eliptical) orbit (a,e,i,etc) and the lat/long of a ground station, how can I calculate the times at which the spacecraft will be visibile from the groundstation (preferably in terms of time since perifocal passage)? I've been googling for hours on this one with no luck so far...

 

[Ed]

On a (related) side note, if I have
x - A sin(x) = B
where A, B are known constants, how do I solve for x?

 

[Jenab]

On a (related) side note, if I have
x - A sin(x) = B
where A, B are known constants, how do I solve for x?

Use a differential calculus procedure for finding roots by successive approximations, such as Newton's method (uses the function and the first derivative) or Danby's method (uses the function and the 1st, 2nd, and 3rd derivatives).

Note that you are seeking the value of a transcendental variable having some arguments inside and outside trigonometric functions. That means you must keep x in radians.

Newton's method.

define:

F(x) = B + A sin x - x

differentiate:

dF/dx = F1(x) = A cos x - 1

You want to find the roots of F(x); i.e., the values of x that make F(x)=0.

Choose an initial guess for x.

x = X0

Repeat while incrementing i,

X(i+1) = X(i) - F(X(i)) / F1(X(i))

Until |X(i+1) - X(i)| approaches zero.

x = X(i+1)

The converged value of X is assigned to x, and it is a root of F(x).

Beware! The root you find might not be the root you need. If F(x) has several roots, you could get the wrong one. It helps to know enough about F(x) so that you can pick a good initial guess for X0.

(I think I've seen Newton's method fail once or twice by trying to find a root at infinity.)

Danby's method.

define:

F(x) = B + A sin x - x

derivatives:

F1(x) = A cos x - 1
F2(x) = -A sin x
F3(x) = -A cos x

Find the roots of F(x).

x = X0

Repeat while incrementing i,

d1 = -F(X(i))/F1(X(i))
d2 = -F(X(i))/{ F1(X(i)) + d1 F2(X(i))/2 }
d3 = -F(X(i))/{ F1(X(i)) + d1 F2(X(i))/2 + d2^2 F3(X(i))/6 }
X(i+1) = X(i) + d3

Until |X(i+1) - X(i)| approaches zero.

x = X(i+1)

The converged value of X is assigned to x, and it is a root of F(x).

Bisection method.

define:

F(x) = B + A sin x - x

choose x0 and x1 such that

F(x0) F(x1) < 0

Repeat

Xmid = (x0 + x1)/2

Q0 = F(x0) F(Xmid)
Q1 = F(x1) F(Xmid)

if Q0<0 then x1 = Xmid
if Q1<0 then x0 = Xmid

Until |x1 - x0| approaches zero, or until F(Xmid)=0.

If F(Xmid)=0 then x=Xmid

else...

x = (x0 + x1)/2

This method is slow! But it doesn't require differentiation.

Peck-Peck-Peck method.

define:

F(x) = B + A sin x - x

choose an initial value, x0, and a small incremental value dx.

Repeat

x1 = x0 + dx

Until F(x1) F(x0) < 0

x ~ (X1 + x0)/2

This method is slow and imprecise. The smaller dx is, the better is the precision but the slower is the procedure. It might be used just to get an initial guess for a more exact procedure.

Jerry Abbott

 

[BobG]

On a (related) side note, if I have
x - A sin(x) = B
where A, B are known constants, how do I solve for x?

This looks suspiciously like a very generic Kepler's equation. x must be the Eccentric Anomaly and A must the eccentricity. B will be your True Anomaly.

Jenab detailed pretty well how you solve this equation. I'd just like to make one note about his comment. For Earth orbiting satellite orbits, use the mean anomaly as your initial guess and you'll have no problems. For your more exotic orbits (comets, interplanetary trajectories, etc), you might run into a situation where the Newton-Raphson method fails to solve the problem (in other words, the Newton method breaks down at very, very high eccentricities). If you use pi as your initial guess in place of the Mean Anomaly, even those orbits will converge (mean anomaly is much faster for the 'typical' satellite orbit).

 

[Jenab]

This looks suspiciously like a very generic Kepler's equation. x must be the Eccentric Anomaly and A must the eccentricity. B will be your True Anomaly.

Jenab detailed pretty well how you solve this equation. I'd just like to make one note about his comment. For Earth orbiting satellite orbits, use the mean anomaly as your initial guess and you'll have no problems. For your more exotic orbits (comets, interplanetary trajectories, etc), you might run into a situation where the Newton-Raphson method fails to solve the problem (in other words, the Newton method breaks down at very, very high eccentricities). If you use pi as your initial guess in place of the Mean Anomaly, even those orbits will converge (mean anomaly is much faster for the 'typical' satellite orbit).

B is the Mean Anomaly. The usual way to write it is

M = u - e sin u

I've seen Newton's method fail to converge Kepler's equation for high-eccentricity orbits.

I'd written a graphical ephemeris code that showed the planets, asteroids, and anything interesting that I had elements for, plotted on a screen. The program was set up so you could arrow-key the program ahead one day at a time. I tried it for the last passage of Halley's comet and noticed that sometimes the comet would disappear from the inner solar system and reappear somewhere out in interstellar space. The problem turned out to be exactly that kind of convergence failure.

I used to fix those problems with reverse interpolation.

Since I knew the correct value of M and wanted the correct value of u, I would set up a peck-peck-peck search on trial u values from 0 to 2 pi radians, with increments of a milliradian, until I found the value of u that returned values for M that bracketed the correct value of M.

I had retained in memory one prior trial point, which gave me

u0 M0
u1 M1
u2 M2

(The known correct value of M is between M1 and M2. The unknown correct value of u must be between u1 and u2.)

Then, using the M's as the independent variable and the u's as the dependent variable, I found the 2nd degree Lagrange interpolating polynomial that included those three points.

u(M) = a M^2 + b M + c

Putting in the correct value of M gave me something tolerably close to the proper value of u.

Jerry Abbott

 

[Jenab]

Now that I've been thinking about it, one might improve the reverse interpolation one step further. Start, as before, with a [0, 2 pi) peck-peck-peck hunt for values of the eccentric anomaly, u, that return values for the mean anomaly, M, which bracket the known correct value of M. Retain the point prior to the pair of bracket points, but continue to gather one point further.

u0 M0
u1 M1
u2 M2
u3 M3

The known correct value of M is between M1 and M2. The unknown correct value of u must be between u1 and u2.

Using the M's as the independent variable and the u's as the dependent variable, find the 2nd degree Lagrange interpolating polynomial that includes points 0, 1, and 2; and the additional 2nd degree Lagrange interpolating polynomial that includes the points 1, 2, and 3.

u1(M) = a1 M^2 + b1 M + c1
u2(M) = a2 M^2 + b2 M + c2

Put the known, correct value of M into each polynomial, and assign the average of the results to the eccentric anomaly.

u = { u1(M) + u2(M) }/2

I would do this in preference to using the same four points to construct a 3rd degree Lagrange interpolating polynomial. You never know how cubics are going to swish around.

Jerry Abbott

 

[BobG]

B is the Mean Anomaly. The usual way to write it is

M = u - e sin u

I've seen Newton's method fail to converge Kepler's equation for high-eccentricity orbits.

Thanks. That will probably work much better, especially since he won't even know the true anomaly, yet.

Try using pi radians for those problem orbits a couple times. That came from a 1998 paper by Charles and Tatum, Celestial Mechanics and Dynamical Astronomy. I tried it for a few examples just out of curiosity and it seemed to work no matter how close the eccentricity came to one. It has a minor drawback for satellite orbits in that low eccentricity orbits will converge slower than if your initial guess was M. But if you work with high eccentricity orbits a lot and use a computer program or spreadsheet, it should be a lot easier than interpolation. If you're using a calculator, it might be a push, since you have to do each of those iterations manually, regardless.

 

[Ed]

Hi folks,
thanks for those replies, that's very useful.

Any ideas to the original problem? I'm currently working with assuming I have a fully defined elliptical orbit about a central body with a groundstation at a given lat & long. At a given time, I can calculate the eccentric anomaly, and hence the position of the satellite .. the issue now becomes finding the angle described between the normal of the orbitted body's surface at the ground station and the line between station and satellite.

Its been a long time since I did this kind of mathematics, and I don't think I tried a problem of this great a complexity. :-/ Am I right in thinking by best bet is to start by converting the positions of the satellite and station into cartesian co-ords and working from there?

 

[Jenab]

Hi folks,
thanks for those replies, that's very useful.

Any ideas to the original problem? I'm currently working with assuming I have a fully defined elliptical orbit about a central body with a groundstation at a given lat & long. At a given time, I can calculate the eccentric anomaly, and hence the position of the satellite .. the issue now becomes finding the angle described between the normal of the orbitted body's surface at the ground station and the line between station and satellite.

Its been a long time since I did this kind of mathematics, and I don't think I tried a problem of this great a complexity. :-/ Am I right in thinking by best bet is to start by converting the positions of the satellite and station into cartesian co-ords and working from there?

The satellite's orbital elements will give you the geocentric position of the satellite with respect to the stars, the inertial (non-rotating) frame of reference.

The Earth rotates. You treat your position on the Earth like an orbit of a different kind, in which your position is a function of your latitude, your longitude, and the time of observation.

Your longitude and the time of observation will give you your local sidereal time. Your local sidereal time is the same as your current right ascension. Your latitude is the same as your declination. Your distance from the center of the Earth is approximately one Earth radius. You resolve these spherical vector components to the equivalent rectangular vector Xyou, Yyou, Zyou.

You use the satellite's elements to predict its geocentric RA & DEC & distance. Then resolve the spherical vector into rectangular components Xsat, Ysat, Zsat.

Both of these vectors are in geocentric celestial coordinates.

You subtract them to get the vector from you to the satellite.

rsat = Rsat - Ryou

You want to know whether the angle, subtended at you, between the center of the Earth and the satellite, is greater than or less than 90 degrees. If it is greater than 90 degrees, the satellite will be above your local horizon and be visible. Otherwise it will be below your horizon and be hidden.

Use the cosine law.

Jerry Abbott

 

[Ed]

OK, thanks again, I'm getting there slowly one step at a time! I'd find this impossible without your help to be sure.

just a quick confirmation of my maths really.. please let me know if I've gone wrong..
I have the orbit in terms of the arguments of ascending node and perigee (n and w)which I can convert into "normal" polar co-ords for right ascension (a) and declination (d) with the inclination (i) as follows:
d = w sin i
a = n + w cos i

which can be substituted into
x = r cos a sin d
y = r sin a sin d
z = r cos d

That gives me the co-ords of the perigee, not the satellite, which is located at angle theta from the line of apsides...

so ... this is where my brain really melts...

RA of the satellite = a + theta cos i
dec of satellite = d + theta sin i

giving
xsat = r cos asat sin dsat
ysat = r sin asat sin dsat
zsat = r cos dsat

or
xsat = r cos (n + w.cos i + theta.cos i) . sin (w.sin i + theta.sin i)
etc

Or have I gone completely and utterly wrong as I suspect?

 

[Jenab]

OK, thanks again, I'm getting there slowly one step at a time! I'd find this impossible without your help to be sure.

just a quick confirmation of my maths really.. please let me know if I've gone wrong..
I have the orbit in terms of the arguments of ascending node and perigee (n and w)which I can convert into "normal" polar co-ords for right ascension (a) and declination (d) with the inclination (i) as follows:
d = w sin i
a = n + w cos i

I'm not sure what you're calculating (a,d) for here. I notice that there is no time dependence on them. The right ascension and declination of an Earth-orbiting satellite would be time dependent.

which can be substituted into
x = r cos a sin d
y = r sin a sin d
z = r cos d

Nope. Declination is a "latitude" kind of angle, not a direction-cosine kind of angle. It is measured as a deviation from the XY plane, which in this case is Earth's equator. It is not measured as a deviation from the Z azis. Those equations should be:

x = r cos a cos d
y = r sin a cos d
z = r sin d

That gives me the co-ords of the perigee, not the satellite, which is located at angle theta from the line of apsides...

Now you tell me.

so ... this is where my brain really melts...

RA of the satellite = a + theta cos i
dec of satellite = d + theta sin i

I don't do it that way. I'm not sure whether that is valid or not.

My way of solving for a position vector involves figuring the current mean anomaly, solving for the eccentric anomaly, using the eccentric anomaly to get the position in orbit canonical coordinates, and rotating it into the coordinate system of standard reference. Maybe someone else around here does it your way.

giving
xsat = r cos asat sin dsat
ysat = r sin asat sin dsat
zsat = r cos dsat

Nope.

xsat = r cos asat cos dsat
ysat = r sin asat cos dsat
zsat = r sin dsat

or
xsat = r cos (n + w.cos i + theta.cos i) . sin (w.sin i + theta.sin i)
etc

I don't know about this.

Jerry Abbott

 

[Jenab]

Lagrange Interpolating Polynomials

Thanks for suggesting pi as an initial guess for Newton's method in solving Kepler's equation for high eccentricity elliptical orbits, Bob.

Since I mentioned Lagrange interpolating polynomials in connection with another way of solving that equation, I'd better describe that curve-fitting technique more fully.

Lagrange interpolating polynomial of degree two.

You're given three points along a segment of some sort of curve whose 2nd derivative doesn't change sign across the interval you want to approximate.

x1, y1
x2, y2
x3, y3

The parabola that fits those three points exactly, and which should approximate the intervening points, is found from:

y(x) = y1 { (x - x2) (x - x3) } / { (x1 - x2) (x1 - x3) }
+ y2 { (x - x1) (x - x3) } / { (x2 - x1) (x2 - x3) }
+ y3 { (x - x1) (x - x2) } / { (x3 - x1) (x3 - x2) }

Which reduces to

y(x) = A x^2 + B x + C

where

A = y1 / {(x1-x2)(x1-x3)} + y2 / {(x2-x1)(x2-x3)} + y3 / {(x3-x1)(x3-x2)}

B = - y1 (x2+x3) / {(x1-x2)(x1-x3)} - y2 (x1+x3) / {(x2-x1)(x2-x3)} - y3 (x1+x2) / {(x3-x1)(x3-x2)}

C = y1 x2 x3 / {(x1-x2)(x1-x3)} + y2 x1 x3 / {(x2-x1)(x2-x3)} + y3 x1 x2 / {(x3-x1)(x3-x2)}

The nice thing about curve-fitting orbits with 2nd degree Lagrange interpolating polynomials is that they fit small arcs of the orbit closely, they are easy to differentiate, and both they and their first derivatives are smooth.

dy/dx = 2Ax + B

Among other things, this means that the derivative of a 2nd degree Lagrange interpolating polynomial can be used to approximate the components of an orbital velocity vector at any point between the first and third points of position, assuming they were well-chosen in relation to the orbit's curvature.

By contrast, estimating the velocity vector at point 2 by dividing the vector different r3-r1 by the time difference t3-t1 requires that the intermediate point occur exactly halfway (in time) between the two endpoints.

That constraint that is avoided when the velocity vector is approximated with the derivative of a 2nd degree Lagrange interpolating polynomial.

Jerry Abbott

 

[Ed]

I'm not sure what you're calculating (a,d) for here. I notice that there is no time dependence on them. The right ascension and declination of an Earth-orbiting satellite would be time dependent.

I'm calculating a,d for the perigee, which isn't time dependent. I realized this after a little while. Hence my edit of the post to add theta later on... it was still a work-in-progress when I initially posted, sorry.

Nope. Declination is a "latitude" kind of angle, not a direction-cosine kind of angle.

That'll teach me not to check stuff. I was staring blankly at my diagram yesterday in the 32° heat and opted to take the easy route and used google. Would you believe the equations I used were taken from a "Model answer" to an exercise set at the university of Alberta... http://www.chem.uAlberta.ca/~abrown/fortran/answers/ex4.html and it was only after seeing them I realized my old maths teachers would've strangled me for not realising it was simple trig. The circles put me off, heh. (or that's my excuse)
Maybe I should send whoever abrown is a polite email and point it out...

I don't do it that way. I'm not sure whether that is valid or not.

My way of solving for a position vector involves figuring the current mean anomaly, solving for the eccentric anomaly, using the eccentric anomaly to get the position in orbit canonical coordinates, and rotating it into the coordinate system of standard reference. Maybe someone else around here does it your way.

hmm. I'm going to be working from a,e,i,n,w,Tp (time since perifocal passage) ... so yeah I can follow you through getting M from Tp and then E from M... but E is of no direct use to me as it is not at the origin of my co-ordinate system? I can convert E to theta for the 2D problem of the elliptical orbit as a function of e, but this is independent of orbit orientation about the focus and so it still leaves me with the issue of projecting this onto the cartesian co-ord system..

Am I right in thinking you're approximating E=theta for a near-circular orbit? I'm actually using an eccentricity of approx 0.5, so the diff between E and theta is quite significant.

Actually... the more I look at my final equation (after the sin-cos mixup) the more I feel confident with it.. it's essentially just adding theta to w, as they are both in the orbital plane, and of course w is a constant. You could rewrite the above as

dsat = (Wsat).sin i
asat = n + (Wsat).sin i
where Wsat = w + theta (angle between rising node and satellite in the orbital plane)

 

[Ed]

For the satellite, this gives:

xsat = r.cos(n+Wsat.cos i).cos(Wsat.sin i)
ysat = r.sin(n+Wsat.cos i).cos(Wsat.sin i)
zsat = r.sin(Wsat.sin i)

which seems to make sense - if you consider the simplified cases when i=0 or 90°
if i=0
xsat = r.cos(n+Wsat)
ysat = r.sin(n+Wsat)
zsat = 0

and if i=90°
xsat = r.cos(n)cos(Wsat)
ysat = r.sin(n)cos(Wsat)
zsat = r.sin(Wsat)

so.. now to see if I can get my head around this initial calculation of E-e.sinE!
Thanks again guys, much appreciated.

 

[BobG]

1. If you divide the time since perigee by the orbital period and multiply by 2 \pi, you get the Mean Anomaly.

M=E-e sin E
You find E using an iterative method. Newton’s method is simplest and most common (followed by Danby’s method or an accelerated Newton’s method).

2. Convert Eccentric Anomaly to True Anomaly

tan\left(\frac{\nu}{2}\right) = \left(\frac{1+e}{1-e}\right)^{\frac{1}{2}} \ast tan \left(\frac{E}{2}\right)
Half angles are used to push your problems to one of the points where you know M = E = \nu. In other words, if the satellite's at apogee all your anomalies are 180 degrees, so don't bother using the equation.

3. To find the satellite radius:

r=\frac{a(1+e^2)}{1+e cos \nu}

4. To find the satellite position in perifocal coordinates (principal direction toward perigee, origin at center of Earth, orbital pole perpendicular to orbit plane)

r_\omega = (r cos \nu)_{e1} +(r sin \nu)_{e2}+0_{e3}

5. Now you need to rotate the satellite vector from perifocal coordinates to geocentric coordinates. Rotate clockwise about the Z axis by the argument of perigee. Clockwise about the X axis by the inclination. Clockwise about the Z axis by the right ascension of ascending node.

If you’re uncomfortable with matrix multiplication, you can combine steps 4 and 5 using three tensor equations (which is what I think you’re trying to do).

x=r(cos\Theta cos\Omega - sin\Theta sin\Omega sin\iota)
y=r(cos\Theta sin\Omega + sin\Theta cos\Omega cos\iota)
z=r(sin\Theta sin\iota)

\Theta = \omega + \nu

The equations you’re using may be equivalent based on some trig identities - I’m just having trouble following your terminology.

 

[BobG]

I do have one other comment on this. The orbital elements you're using are only valid for the time that observation was made. Because of atmospheric drag (for low orbiting satellites), the Earth's triaxiality (non-spherical gravitational field - especially that bulge around the equator), and the Sun and Moon pulling on the satellite (if a high altitude satellite where the Earth's gravity is weak), you need to adjust your orbital elements for the time you want to observe the satellite.

For example, on sun-synchronous satellites, the right ascension of ascending node drifts to the east about .986 degrees per day. Makes sense - 360 degrees in a circle and 365 days in a year.

This link explains how to adjust your element sets for the Earth's oblateness.
http://scienceworld.wolfram.com/physics/OrbitalPerturbation.html
As imposing as it looks, the bottom line is at the bottom. In other words, the two main equations you need to worry about are for:

\frac{d\omega}{dt}
and
\frac{d\Omega}{dt}

The other perturbations are small enough you can compensate for them by getting new element sets at least once a week (the ones available on the internet are already old, so they don't have a very long shelf life).

Atmospheric drag will also be a problem for low orbiting satellites. There is no easy way to compensate for that, since atmospheric drag varies in response to solar activity. Just accept the fact that your predicts might be off by up to 30 seconds or more during solar maximum. Fortunately we just passed the solar maximum, so your predicts will get better and better until the next one around 2010 or 2011.

That gives you time to learn all about it. I once made a spreadsheet that very, VERY roughly modeled atmospheric drag, but even good models don't predict rise/fade times very good for very long during the solar maximum.

 

[Jenab]

I'm going to be working from a,e,i,n,w,Tp (time since perifocal passage) ... so yeah I can follow you through getting M from Tp and then E from M... but E is of no direct use to me as it is not at the origin of my co-ordinate system?

You can get x,y,z from the eccentric anomaly, u, (and the elements) directly, without getting the true anomaly first.

x''' = a (cos u - e)

y''' = a sin u (1-e^2)^0.5

z''' = 0

Rotation by the argument of the perihelion, w.

x'' = x''' cos w - y''' sin w

y'' = x''' sin w + y''' cos w

z'' = z''' = 0

Rotation by the inclination, i.

x' = x''

y' = y'' cos i

z' = y'' sin i

Rotation by the longitude of ascending node, L.

x = x' cos L - y' sin L

y = x' sin L + y' cos L

z = z'

The unprimed position vector, [x, y, z], is the geocentric celestial rectangular position vector to the satellite at the time for which the input eccentric anomaly (u) is correct, provided that your elements are referred to the geocentric celestial coordinate system.

Jerry Abbott