Kepler's Equation to calculate transfer time

[brownnr]

I'm working out of a book called "Analytical Mechanics of Space Systems," by Junkins and Schaub.

In that book, they give a modified shooting method method for solving Lambert's two-point boundary value problem (connecting two points in space with a closed orbit). One of the problems that I face is that I am getting the right initial velocity that connects the two points, but AFTER I have that I am not getting the correct transfer time. Essentially, I have eccentricity e, true anomaly f1 (corresponding to the first point), and f2 (true anomaly of the second point), length of the semi-major axis a, and the gravitational parameter mu.

I convert the true anomaly f1, f2 to eccentric anomaly E1, E2 using the following code:

// Find Eccentric Anomaly given eccentricity and true anomaly
double convertftoE( double e, double f){

double E = 0.0;

E = 2*atan(sqrt((1-e)/(1+e))*tan(f/2));

return E;
}

Then I convert eccentric anomalies E1, E2 to mean anomalies M1, M2 using:

// Calculates the mean anomaly given e (eccentricity)
// and E (eccentric anomaly)
double calcMeanAnomaly(double e, double E) {

double M = 0.0;

M = E - e*sin(E);

return M;

}

Then I calculate the transfer time dt_hat using Kepler's equation:

dt_hat = sqrt( (a*a*a) / mu ) * ( M2 - M1 );

I integrate using the initial time over dt_hat with a 4th order Runge Kutta method, but it always ends up short. If I let it go longer, the orbit goes right overtop of the point that I am trying to reach.

I need to have a "nice" equation or method for finding the proper transfer time, because I use that dt_hat to move through a number of neighboring solutions to the sought after transfer time Ts through a linear parameter 0 <= alpha <=1: dt_hat*(1 - alpha) + alpha*Ts, where alpha slowly increments from 0 to 1 changing the velocity and orbit each time via a numerically calculated Jacobian acting as a sensitivity matrix to get the satellite from point A to point B in the required amount of time Ts.

The way I have it set now, I using the results from the RK4 method as a starting dt. But this should be something I can calculate prior to using a numerical method such as RK4.

Can anyone spot a mistake in what I'm doing in calculating the transfer time? Am I misinterpreting the proper use of these equations? Thanks very much.

-Nelson

 

[D H]

The mean anomaly at some later time is easily calculable as mean anomaly is a linear function of time. Use Kepler's equation backwards to find the eccentric anomaly, and then find the true anomaly.

What you are missing is the inverse of Kepler's equation, that is, determine the eccentric anomaly given the mean anomaly. One easy way to do this is via Newton's method.

Regarding the use of RK4: You either need to ensure you have an integral number of steps or do something special on the final step to make sure you hit $t_0+\Delta \hat t$ exactly. Even if you do this, you will still have some problems because an RK4 integrator does not conserve angular momentum or orbital energy. You can overcome these problems to some extent by making your time step sufficiently small. On the other hand, if you make the time step too small you will run into two problems. First, it will take a long time to perform the integration. This will impinge on the number of interations you can make in your Lambert targeting algorithm. Second, and even worse, you will find the error will first decrease and the increase as you make the time step ever smaller.

Since you are working in a domain with an exact solution (the two body problem), I recommend against using numerical integration.

 

[brownnr]

Some modifications.

I apologize for having put the original post in the other forum. There wasn't any malicious intent, I posted it first in general astronomy and thought that the question was more suited here after returning to the main page. I don't come to the forums all that often.

I would completely agree with your assessment of the numeric integration. But for someone who has no experience with the analytic solution, I was using it to "fill in the gaps" and check my answers until I could get it working. I keep filling in more pieces of the analytic solution each time I progress with this and it greatly speeds up the convergence rates each time I do so.

In the context of the problem, you don't start from an observers point of view, you start with two vectors w.r.t. the main body. So you have the true anomaly, and a method of calculating the true anomaly difference (angle between the vectors). To calculate the transfer time you need the mean anomaly (and hence the eccentric anomaly).

I made a slight modification in finding the eccentric anomaly:

// Find Eccentric Anomaly given eccentricity and true anomaly
double convertftoE( double e, double f){

double E = 0.0;

E = atan(sqrt((1-e)/(1+e))*tan(f/2));

if(E < 0.0) E+=3/2*PI;

E *= 2.0;

return E;
}

As atan returns a range of [-Pi/2, Pi/2], I had to offset by two quadrants to get values between [0, 2*Pi]. It seems to work nicely now.

Thank you for your input, and sorry again about the redundent post.