I'm creating an orbit simulator using Java. I have elliptical orbits working great, and now I'm trying to do hyperbolic orbits. The problem I'm having is that when the ship enters the orbit at periapsis, it moves very slowly at first, and quickly accelerates until it reaches infinite velocity at infinite altitude. Obviously, it should start at its fastest and go to zero, which is the opposite of what it's doing. The orbit path is rendering properly and the ship is following it, but the problem seems to be that Kepler's equation is growing exponentially instead of converging on a value of the hyperbolic anomaly, and I don't know why. I'm following my references precisely, and the equivalent code for e < 1 works perfectly. Here is the function that calculates the ship's position at time t, with n being the mean motion.
if(e > 1) {
M = n * t; //current mean anomaly
//iterate Kepler's equation to converge on a value of the current hyperbolic anomaly
for(int i = 0; i < KEPLER; i++) {H = (e * sinh(H)) - M;}
nu = 2 * atan((tanh(H / 2)) * sqrt((e + 1) / (e - 1))); //current true anomaly
//position variables
r = (a * (1.0 - (e * e))) / (1.0 + (e * cos(nu))); //current orbital radius
phi = atan((e * sin(nu)) / (1.0 + (e * cos(nu)))); //current flight angle
v = sqrt(GM * ((2.0 / r) - (1.0 / a))); //current orbital velocity
}
I've quadruple checked every equation, verified the values of the variables, and gone over each step with a fine toothed comb. I think the problem is that H = (e * sinh(H)) - M is growing exponentially with every iteration instead of converging, but I don't think that's what it's supposed to do. I wrote this simple test program to watch it work, and this is what it gave me.
public class Main {
public static void main(String[] args) throws java.io.IOException {
double E = 0, H = 0, E0 = 0, H0 = 0;
int k = 0, k1 = 11;
double M = 0.04;
double e1 = 0.1;
double e2 = 2.0;
while(k < k1) {
E = M + (e1 * Math.sin(E));
H = (e2 * Math.sinh(H)) - M;
System.out.printf("k = %d, E = %.6f, H = %.6f\n", k, E, H);
System.out.printf(" dE = %.6f, dH = %.6f\n", E - E0, H - H0);
E0 = E;
H0 = H;
k++;
try {
Thread.sleep(100);
} catch (InterruptedException e) {
e.printStackTrace();
}
}
}
}
Output:
k = 0, E = 0.040000, H = -0.040000
dE = 0.040000, dH = -0.040000
k = 1, E = 0.043999, H = -0.120021
dE = 0.003999, dH = -0.080021
k = 2, E = 0.044398, H = -0.280619
dE = 0.000400, dH = -0.160598
k = 3, E = 0.044438, H = -0.608634
dE = 0.000040, dH = -0.328014
k = 4, E = 0.044442, H = -1.333825
dE = 0.000004, dH = -0.725191
k = 5, E = 0.044443, H = -3.572066
dE = 0.000000, dH = -2.238241
k = 6, E = 0.044443, H = -35.601966
dE = 0.000000, dH = -32.029899
k = 7, E = 0.044443, H = -2895592024320601.000000
dE = 0.000000, dH = -2895592024320565.500000
k = 8, E = 0.044443, H = -Infinity
dE = 0.000000, dH = -Infinity
k = 9, E = 0.044443, H = -Infinity
dE = 0.000000, dH = NaN
k = 10, E = 0.044443, H = -Infinity
dE = 0.000000, dH = NaN
I'm trying to teach this to myself, so I'm hoping it's something obvious like using the wrong equation or skipping a step.
Trying to find a solution of X = f(X) by just iterating it as an assignment works only if d(f)/dX < 1 ... just use a different numerical method, like for example https://en.wikipedia.org/wiki/Newton%27s_method