Created
December 18, 2010 01:40
-
-
Save steos/746022 to your computer and use it in GitHub Desktop.
calculate orbital satellite position using Kepler's laws
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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