calculate orbital satellite position using Kepler's laws

Kepler.java

import static java.lang.Math.*;

public class Kepler
{
	public static class Point2d
	{
		public double x;
		public double y;
		
		public Point2d(double x, double y)
		{
			this.x = x;
			this.y = y;
		}

		public Point2d toCartesian()
		{
			return new Point2d((x * cos(y)), 
				(x * sin(y)));
		}	
	}
	
	/**
	 * gravitational constant
	 */
	public static final double G = 6.67428E-11;

	/**
	 * 2 * pi
	 */
	public static final double TWO_PI = PI * 2;

	/**
	 * pi squared
	 */
	public static final double PI2 = pow(PI, 2);

	/**
	 * solves kepler (find position of satellite for given point in time)
	 * 
	 * @param majorAxis the semi-major axis of the satellites orbit
	 * @param e the numeric eccentricity of the satellites orbit
	 * @param T the total time of circulation
	 * @param time time since satellite reached periapsis
	 * @param iterations number of iterations to calculate eccentric anomaly
	 * 
	 * @return the cartesian coordinates of the planet in the orbital plane
	 */
	public static Point2d getOrbitalPosition(double majorAxis, double e, 
			double T, double time, int iterations)
	{
		double mean = getMeanAnomaly(T, time);
		double ecc = getEccentricAnomaly(e, mean, iterations);
		double phi = getTrueAnomaly(e, ecc);
		double r = getRadiusLength(majorAxis, e, phi);
		
		return new Point2d(r, phi).toCartesian();
	}

	/**
	 * returns the circulation time according to
	 * keplers third law
	 * 
	 * @param majorAxis the semi-major axis of the orbit
	 * @param M the mass of the central body
	 * @param m the satellites mass
	 * 
	 * @return the total time of circulation
	 */
	public static double getCirculationTime(double majorAxis, double M, double m)
	{
		return sqrt((pow(majorAxis, 3) * 4 * PI2) / 
				(G * (M + m)) );
	}

	/**
	 * iteratively calculates the eccentric anomaly of an orbit for a given
	 * mean anomaly and eccentricity (converges very fast - at about 4 iterations)
	 */
	public static double getEccentricAnomaly(double eccentricity, double meanAnomaly, int iterations)
	{
		double tempResult = meanAnomaly;
		double numerator;

		for (int i = 0; i < iterations; i++)
		{
			numerator = meanAnomaly + eccentricity * sin(tempResult)
					- eccentricity * tempResult * cos(tempResult);
			tempResult = numerator / (1 - eccentricity * cos(tempResult));
		}

		return tempResult;
	}

	/**
	 * calculates the mean anomaly of an orbit for the specified point in
	 * time
	 * 
	 * @param circulationTime total 
	 * @param time time since periapsis
	 */
	public static double getMeanAnomaly(double circulationTime, double time)
	{
		return TWO_PI / circulationTime * time;
	}

	/**
	 * calculates the true anomaly for a given orbits eccentricity and
	 * eccentric anomaly
	 */
	public static double getTrueAnomaly(double eccentricity, double eccentricAnomaly)
	{
		return 2 * atan(sqrt((1 + eccentricity) / (1 - eccentricity))
				* tan(eccentricAnomaly / 2));
	}

	/**
	 * calculates the length of the radius vector
	 */
	public static double getRadiusLength(double majorAxis, double eccentricity, double trueAnomaly)
	{
		double nominator = majorAxis * (1 - pow(eccentricity, 2));
		double denominator = 1 + eccentricity * cos(trueAnomaly);

		return nominator / denominator;
	}
}