Writing a matlab program to solve the Kepler Problem (Universal Vars)

[Julez55]

TL;DR Summary

I'm writing a matlab program to solve the Kepler problem for all orbit cases. It works successfully for Elliptic, Circular, and Hyperbolic. The case of hyperbolic orbits where a blows up is giving me trouble.

Introduction: I am writing a Matlab program that solves the Kepler problem for any Elliptical, Hyperbolic, or Circular orbit. Given an initial position vector, velocity vector, find the position vector and velocity vector after some given time. I am not doing this for any class, I have been doing this to work on my computer science abilities. The textbook I am using is "Fundamentals of Astrodynamics" by Bate, Mueller, and White. The specific program I am working on is a suggested project in the back of the book.

The Problem: My program successfully solves the Kepler problem for all elliptical, hyperbolic, and circular orbits. However, I'm having a hard time solving the case of Parabolic orbits where the denominator of a goes to zero. I tried to resolve this by calculating the reciprocal of a, which then, in this case, is zero. This leads to the initial guess of x_n to be zero, and the initial guess of z_n to be zero as well, which causes the iteration to blow up. I also added a case inside the iteration which computes the values of C_n and S_n using their series representations when the reciprocal of a, called alpha in my program is equal to zero. I attached my Matlab program below.

If anybody thinks they can help me with this I would love to have further conversation, I'd really like to figure this one out!

 

[Dr_Nate]

It would really help if you wrote out here the equations of interest. Having a brief look at your code, I believe you are trying to calculate the eccentricity. My guess is that you are using the equation that only works with $0<e<1$, which is for an ellipse.