Skip to content

Instantly share code, notes, and snippets.

@steos
Created December 18, 2010 01:40
Show Gist options
  • Save steos/746022 to your computer and use it in GitHub Desktop.
Save steos/746022 to your computer and use it in GitHub Desktop.
calculate orbital satellite position using Kepler's laws
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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment