// Orbital parameters
    public double[] semimajorAxis = new double[2]; //[au, au/century]
    public double[] eccentricity = new double[2]; //[ , /century]
    public double[] inclination = new double[2]; //inclination [degrees, degrees/century]
    public double[] meanLongitude = new double[2]; //inclination [degrees, degrees/century]
    public double[] longitudeOfPerihelion = new double[2]; //inclination [degrees, degrees/century]
    public double[] longitudeOfTheAscendingNode = new double[2]; //inclination [degrees, degrees/century]

void ComputePositionAtDate(){
        double tMillisFromJ2000 = DateTime.Now.ToUniversalTime().Subtract(new DateTime(2000, 1, 1, 12, 0, 0, DateTimeKind.Utc)).TotalMilliseconds;
        double tCenturiesFromJ2000 = tMillisFromJ2000 / (1000 * 60 * 60 * 24 * 365.25 * 100);
        double a = semimajorAxis [0] + semimajorAxis [1] * tCenturiesFromJ2000;
        double e = eccentricity [0] + eccentricity [1] * tCenturiesFromJ2000;
        double i = inclination [0] + inclination [1] * tCenturiesFromJ2000;
        double L = Math.IEEERemainder(Math.IEEERemainder(meanLongitude[1] * tCenturiesFromJ2000, 360) + meanLongitude[0], 360);

        Debug.Log (this.name + ": " + L);

        double p = longitudeOfPerihelion [0] + longitudeOfPerihelion [1] * tCenturiesFromJ2000;
        double W = longitudeOfTheAscendingNode [0] + longitudeOfTheAscendingNode [1] * tCenturiesFromJ2000;

        double M = L - p;
        double w = p - W;

        double E = M;
        while (true) {
            double dE = (E - e * Mathd.Sin (E) - M) / (1 - e * Mathd.Cos (E));
            E -= dE;
            if (Mathd.Abs (dE) < 1e-6)
                break;
        }

        double P = a * (Mathd.Cos (E) - e);
        double Q = a * Mathd.Sin (E) * Math.Sqrt (1 - Mathd.Pow (e, 2d));

        // rotate by argument of periapsis
        double x = Mathd.Cos(w) * P - Mathd.Sin(w) * Q;
        double y = Mathd.Sin (w) * P + Mathd.Cos (w) * Q;
        // rotate by inclination
        double z = Mathd.Sin(i) * x;
        x = Mathd.Cos (i * x);
        // rotate by longitude of ascending node
        double xTemp = x;
        x = Mathd.Cos (W) * xTemp - Mathd.Sin (W) * y;
        y = Mathd.Sin (W) * xTemp + Mathd.Cos (W) * y;

        position = new Vector3d (x * 149597870700d, z * 149597870700d, y * 149597870700d);
    }