private double[] getMoon() {
        // MOON PARAMETERS (Formulae from "Calendrical Calculations")
        double phase = normalizeRadians((297.8502042 + 445267.1115168 * t - 0.00163 * t * t + t * t * t / 538841 - t * t * t * t / 65194000) * DEG_TO_RAD);

        // Anomalistic phase
        double anomaly = (134.9634114 + 477198.8676313 * t + .008997 * t * t + t * t * t / 69699 - t * t * t * t / 14712000);
        anomaly = anomaly * DEG_TO_RAD;

        // Degrees from ascending node
        double node = (93.2720993 + 483202.0175273 * t - 0.0034029 * t * t - t * t * t / 3526000 + t * t * t * t / 863310000);
        node = node * DEG_TO_RAD;

        double E = 1.0 - (.002495 + 7.52E-06 * (t + 1.0)) * (t + 1.0);

        // Now longitude, with the three main correcting terms of evection,
        // variation, and equation of year, plus other terms (error<0.01 deg)
        // P. Duffet's MOON program taken as reference
        double l = (218.31664563 + 481267.8811958 * t - .00146639 * t * t + t * t * t / 540135.03 - t * t * t * t / 65193770.4);
        l += 6.28875 * Math.sin(anomaly) + 1.274018 * Math.sin(2 * phase - anomaly) + .658309 * Math.sin(2 * phase);
        l += 0.213616 * Math.sin(2 * anomaly) - E * .185596 * Math.sin(sanomaly) - 0.114336 * Math.sin(2 * node);
        l += .058793 * Math.sin(2 * phase - 2 * anomaly) + .057212 * E * Math.sin(2 * phase - anomaly - sanomaly) + .05332 * Math.sin(2 * phase + anomaly);
        l += .045874 * E * Math.sin(2 * phase - sanomaly) + .041024 * E * Math.sin(anomaly - sanomaly) - .034718 * Math.sin(phase) - E * .030465 * Math.sin(sanomaly + anomaly);
        l += .015326 * Math.sin(2 * (phase - node)) - .012528 * Math.sin(2 * node + anomaly) - .01098 * Math.sin(2 * node - anomaly) + .010674 * Math.sin(4 * phase - anomaly);
        l += .010034 * Math.sin(3 * anomaly) + .008548 * Math.sin(4 * phase - 2 * anomaly);
        l += -E * .00791 * Math.sin(sanomaly - anomaly + 2 * phase) - E * .006783 * Math.sin(2 * phase + sanomaly) + .005162 * Math.sin(anomaly - phase) +
             E * .005 * Math.sin(sanomaly + phase);
        l += .003862 * Math.sin(4 * phase) + E * .004049 * Math.sin(anomaly - sanomaly + 2 * phase) + .003996 * Math.sin(2 * (anomaly + phase)) +
             .003665 * Math.sin(2 * phase - 3 * anomaly);
        l += E * 2.695E-3 * Math.sin(2 * anomaly - sanomaly) + 2.602E-3 * Math.sin(anomaly - 2 * (node + phase));
        l += E * 2.396E-3 * Math.sin(2 * (phase - anomaly) - sanomaly) - 2.349E-3 * Math.sin(anomaly + phase);
        l += E * E * 2.249E-3 * Math.sin(2 * (phase - sanomaly)) - E * 2.125E-3 * Math.sin(2 * anomaly + sanomaly);
        l += -E * E * 2.079E-3 * Math.sin(2 * sanomaly) + E * E * 2.059E-3 * Math.sin(2 * (phase - sanomaly) - anomaly);
        l += -1.773E-3 * Math.sin(anomaly + 2 * (phase - node)) - 1.595E-3 * Math.sin(2 * (node + phase));
        l += E * 1.22E-3 * Math.sin(4 * phase - sanomaly - anomaly) - 1.11E-3 * Math.sin(2 * (anomaly + node));
        double longitude = l;

        // Let's add nutation here also
        double M1 = (124.90 - 1934.134 * t + 0.002063 * t * t) * DEG_TO_RAD;
        double M2 = (201.11 + 72001.5377 * t + 0.00057 * t * t) * DEG_TO_RAD;
        double d  = -.0047785 * Math.sin(M1) - .0003667 * Math.sin(M2);
        longitude += d;

        // Get accurate Moon age
        double Psin = 29.530588853;
        moonAge = normalizeRadians((longitude - slongitude) * DEG_TO_RAD) * Psin / TWO_PI;

        // Now Moon parallax
        double parallax = .950724 + .051818 * Math.cos(anomaly) + .009531 * Math.cos(2 * phase - anomaly);
        parallax += .007843 * Math.cos(2 * phase) + .002824 * Math.cos(2 * anomaly);
        parallax += 0.000857 * Math.cos(2 * phase + anomaly) + E * .000533 * Math.cos(2 * phase - sanomaly);
        parallax += E * .000401 * Math.cos(2 * phase - anomaly - sanomaly) + E * .00032 * Math.cos(anomaly - sanomaly) - .000271 * Math.cos(phase);
        parallax += -E * .000264 * Math.cos(sanomaly + anomaly) - .000198 * Math.cos(2 * node - anomaly);
        parallax += 1.73E-4 * Math.cos(3 * anomaly) + 1.67E-4 * Math.cos(4 * phase - anomaly);

        // So Moon distance in Earth radii is, more or less,
        double distance = 1.0 / Math.sin(parallax * DEG_TO_RAD);

        // Ecliptic latitude with nodal phase (error<0.01 deg)
        l = 5.128189 * Math.sin(node) + 0.280606 * Math.sin(node + anomaly) + 0.277693 * Math.sin(anomaly - node);
        l += .173238 * Math.sin(2 * phase - node) + .055413 * Math.sin(2 * phase + node - anomaly);
        l += .046272 * Math.sin(2 * phase - node - anomaly) + .032573 * Math.sin(2 * phase + node);
        l += .017198 * Math.sin(2 * anomaly + node) + .009267 * Math.sin(2 * phase + anomaly - node);
        l += .008823 * Math.sin(2 * anomaly - node) + E * .008247 * Math.sin(2 * phase - sanomaly - node) + .004323 * Math.sin(2 * (phase - anomaly) - node);
        l += .0042 * Math.sin(2 * phase + node + anomaly) + E * .003372 * Math.sin(node - sanomaly - 2 * phase);
        l += E * 2.472E-3 * Math.sin(2 * phase + node - sanomaly - anomaly);
        l += E * 2.222E-3 * Math.sin(2 * phase + node - sanomaly);
        l += E * 2.072E-3 * Math.sin(2 * phase - node - sanomaly - anomaly);
        double latitude = l;

        return new double[]{longitude, latitude, distance * EARTH_RADIUS / AU, Math.atan(1737.4 / (distance * EARTH_RADIUS))};
    }
Calculates everything for the Sun and the Moon.
src/main/java/eu/hansolo/tilesfx/tools/SunMoonCalculator.java