public double[] getMoonDiskOrientationAngles() {
        double outS[]    = doCalc(getSun());
        double moonPos[] = getMoon();
        double outM[]    = doCalc(moonPos);
        double moonLon   = moonPos[0], moonLat = moonPos[1], moonRA = outM[6], moonDEC = outM[7];
        double sunRA     = outS[6], sunDEC = outS[7];

        // Moon's argument of latitude
        double F = (93.2720993 + 483202.0175273 * t - 0.0034029 * t * t - t * t * t / 3526000.0 + t * t * t * t / 863310000.0) * DEG_TO_RAD;
        // Moon's inclination
        double I = 1.54242 * DEG_TO_RAD;
        // Moon's mean ascending node longitude
        double omega = (125.0445550 - 1934.1361849 * t + 0.0020762 * t * t + t * t * t / 467410.0 - t * t * t * t / 18999000.0) * DEG_TO_RAD;
        // Obliquity of ecliptic (approx, better formulae up)
        double eps = 23.43929 * DEG_TO_RAD;

        // Obtain optical librations lp and bp
        double W     = moonLon - omega;
        double sinA  = Math.sin(W) * Math.cos(moonLat) * Math.cos(I) - Math.sin(moonLat) * Math.sin(I);
        double cosA  = Math.cos(W) * Math.cos(moonLat);
        double A     = Math.atan2(sinA, cosA);
        double lp    = normalizeRadians(A - F);
        double sinbp = -Math.sin(W) * Math.cos(moonLat) * Math.sin(I) - Math.sin(moonLat) * Math.cos(I);
        double bp    = Math.asin(sinbp);

        // Obtain position angle of axis p
        double x    = Math.sin(I) * Math.sin(omega);
        double y    = Math.sin(I) * Math.cos(omega) * Math.cos(eps) - Math.cos(I) * Math.sin(eps);
        double w    = Math.atan2(x, y);
        double sinp = Math.sqrt(x * x + y * y) * Math.cos(moonRA - w) / Math.cos(bp);
        double p    = Math.asin(sinp);

        // Compute bright limb angle bl
        double bl = (Math.PI + Math.atan2(Math.cos(sunDEC) * Math.sin(moonRA - sunRA),
                                          Math.cos(sunDEC) * Math.sin(moonDEC) * Math.cos(moonRA - sunRA) - Math.sin(sunDEC) * Math.cos(moonDEC)));

        // Paralactic angle par (first obtain local apparent sidereal time)
        double jd0   = Math.floor(jd_UT - 0.5) + 0.5;
        double T0    = (jd0 - J2000) / JULIAN_DAYS_PER_CENTURY;
        double secs  = (jd_UT - jd0) * SECONDS_PER_DAY;
        double gmst  = (((((-6.2e-6 * T0) + 9.3104e-2) * T0) + 8640184.812866) * T0) + 24110.54841;
        double msday = 1.0 + (((((-1.86e-5 * T0) + 0.186208) * T0) + 8640184.812866) / (SECONDS_PER_DAY * JULIAN_DAYS_PER_CENTURY));
        gmst = (gmst + msday * secs) * (15.0 / 3600.0) * DEG_TO_RAD;
        double lst = gmst + latitude;

        y = Math.sin(lst - moonRA);
        x = Math.tan(longitude) * Math.cos(moonDEC) - Math.sin(moonDEC) * Math.cos(lst - moonRA);
        double par = x != 0 ? Math.atan2(y, x) : (y / Math.abs(y)) * PI_OVER_TWO;
        return new double[]{ lp, bp, p, bl, par };
    }
Calculates everything for the Sun and the Moon.
src/main/java/eu/hansolo/tilesfx/tools/SunMoonCalculator.java