I am trying to model a two dimensional orbit for a two body Kepler problem but have gotten stuck when introducing the time variable.
For a satellite with known semi major axis (a), eccentricity (e), and true anomaly (theta), I have:
r = a*(1-e**2)/(1+e*cos(theta))
How can I calculate theta as a function of time using the Wiki suggested Universal variable formulation method? I have no idea how to implement (am using Python but any algorithm advice much appreciated!)
Alternatively, how do I calculate r as a function of time?
Note: All other orbital elements and masses are available. Also I am trying to come up with a general solution for elliptical, hyperbolic and parabolic orbits.
Cheers!
Answer
I calculated true anomaly as function of time, for planetary motion , in c# , in this way:
Compute mean anomaly (time: current time , G: newton grav.connstant, M: planet mass or the sum of the two orbiting objects , a: semi major axis)
//M = nt
double n = Math.Sqrt((G * (M)) / (a * a * a));
double Mt = n * time;Compute the eccentric anomaly E by solving Kepler's equation:
//For orbits with ε > 0.8, an initial value of E0 = π should be used.
if (eccentr>0.8)
E = NumApprox(150, Math.PI,Mt, 10E-15);
else
E = NumApprox(150, Mt, Mt, 10E-15);true anomaly (angle)
true_anom = 2.0 * Math.Atan2(Math.Sqrt(1.0 + eccentr) * Math.Sin(E / 2.0), Math.Sqrt(1.0 - eccentr) * Math.Cos(E / 2.0));
distance from planet
d = a * ((1.0 - eccentr * eccentr) / (1.0 + eccentr * Math.Cos(true_anom)));
Finaly Numerical approximation of inverse problem:
private double NumApprox(int intr, double prev,double Mt, double err)
{
double ret = prev;
double retprev = prev;
for (int i=0 ; i retprev = ret;
ret = ret - (ret - eccentr * Math.Sin(ret) - Mt) / (1.0 - eccentr * Math.Cos(ret));
if ( Math.Abs(ret - retprev) < err)
break;
}
return ret;
}
EDIT: calculate Position and Velocity
What we done is half of the work in getting Cartesain Orbit Elements from kepler Orbit Elements where :
//some kepler Orbit Elements:
public double d;
public double true_anom;
public double eccentr;
public double a;
public double E;
public double w=0; //small omega ω : Argument of periapsis (in rad)
//Cartesain Orbit Elements:
//Position Vector
private double x;
private double y;
private double z; //2d :not used
//Velocity Vector
private double vx;
private double vy;
private double vz; //2d :not used
as final step we can calculate position and velocity vectors as :
//Position
x = d * Math.Cos(true_anom);
y = d * Math.Sin(true_anom);
//apply ω
double xx = x * Math.Cos(w) - y * Math.Sin(w);
double yy = x * Math.Sin(w) + y * Math.Cos(w);
x = xx;
y = yy;
//Velocity
double v = Math.Sqrt(G * M * a) / d;
vx = -v * Math.Sin(E);
vy = -v * Math.Sqrt(1.0-eccentr*eccentr) * Math.Cos(E);
Newton's method for numerical approximation
The Kerbal Space Program (KSP) Physics Documantation (pdf)
No comments:
Post a Comment