Tuesday, January 5, 2016

physics - Orbital mechanics: orbit as a function of time. Universal variable formulation?


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:




  1. 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;


  2. 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);


  3. 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));





  4. 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;
}

enter image description here



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);

EDIT: references enter image description here


Orbital elements


Kepler's equation


True anomaly



Eccentricity vector


Argument of periapsis


Newton's method for numerical approximation


The Kerbal Space Program (KSP) Physics Documantation (pdf)


No comments:

Post a Comment

Simple past, Present perfect Past perfect

Can you tell me which form of the following sentences is the correct one please? Imagine two friends discussing the gym... I was in a good s...