private double[] doCalc(final double[] POS) {
        // Ecliptic to equatorial coordinates
        double t2  = this.t / 100.0;
        double tmp = t2 * (27.87 + t2 * (5.79 + t2 * 2.45));
        tmp = t2 * (-249.67 + t2 * (-39.05 + t2 * (7.12 + tmp)));
        tmp = t2 * (-1.55 + t2 * (1999.25 + t2 * (-51.38 + tmp)));
        tmp = (t2 * (-4680.93 + tmp)) / 3600.0;
        double angle = (23.4392911111111 + tmp) * DEG_TO_RAD; // obliquity

        // Add nutation in obliquity
        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  = .002558 * Math.cos(M1) - .00015339 * Math.cos(M2);
        angle += d * DEG_TO_RAD;

        POS[0] *= DEG_TO_RAD;
        POS[1] *= DEG_TO_RAD;
        double cl = Math.cos(POS[1]);
        double x  = POS[2] * Math.cos(POS[0]) * cl;
        double y  = POS[2] * Math.sin(POS[0]) * cl;
        double z  = POS[2] * Math.sin(POS[1]);
        tmp = y * Math.cos(angle) - z * Math.sin(angle);
        z = y * Math.sin(angle) + z * Math.cos(angle);
        y = tmp;

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

        // Obtain topocentric rectangular coordinates
        // Set radiusAU = 0 for geocentric calculations
        // (rise/set/transit will have no sense in this case)
        double radiusAU     = EARTH_RADIUS / AU;
        double correction[] = new double[]{radiusAU * Math.cos(longitude) * Math.cos(lst), radiusAU * Math.cos(longitude) * Math.sin(lst), radiusAU * Math.sin(longitude)};
        double xtopo        = x - correction[0];
        double ytopo        = y - correction[1];
        double ztopo        = z - correction[2];

        // Obtain topocentric equatorial coordinates
        double ra  = 0.0;
        double dec = PI_OVER_TWO;
        if (ztopo < 0.0) dec = -dec;
        if (ytopo != 0.0 || xtopo != 0.0) {
            ra = Math.atan2(ytopo, xtopo);
            dec = Math.atan2(ztopo / Math.sqrt(xtopo * xtopo + ytopo * ytopo), 1.0);
        }
        double dist = Math.sqrt(xtopo * xtopo + ytopo * ytopo + ztopo * ztopo);

        // Hour angle
        double angh = lst - ra;

        // Obtain azimuth and geometric alt
        double sinlat = Math.sin(longitude);
        double coslat = Math.cos(longitude);
        double sindec = Math.sin(dec), cosdec = Math.cos(dec);
        double h      = sinlat * sindec + coslat * cosdec * Math.cos(angh);
        double alt    = Math.asin(h);
        double azy    = Math.sin(angh);
        double azx    = Math.cos(angh) * sinlat - sindec * coslat / cosdec;
        double azi    = Math.PI + Math.atan2(azy, azx); // 0 = north

        // Get apparent elevation
        if (alt > -3 * DEG_TO_RAD) {
            double r    = 0.016667 * DEG_TO_RAD * Math.abs(Math.tan(PI_OVER_TWO - (alt * RAD_TO_DEG + 7.31 / (alt * RAD_TO_DEG + 4.4)) * DEG_TO_RAD));
            double refr = r * (0.28 * 1010 / (10 + 273.0)); // Assuming pressure of 1010 mb and T = 10 C
            alt = Math.min(alt + refr, PI_OVER_TWO); // This is not accurate, but acceptable
        }

        switch (twilight) {
            case GOLDEN_HOUR:
                tmp = 6 * DEG_TO_RAD;
                break;
            case HORIZON_34arcmin:
                // Rise, set, transit times, taking into account Sun/Moon angular radius (pos[3]).
                // The 34' factor is the standard refraction at horizon.
                // Removing angular radius will do calculations for the center of the disk instead
                // of the upper limb.
                tmp = -(34.0 / 60.0) * DEG_TO_RAD - POS[3];
                break;
            case BLUE_HOUR:
                tmp = -4 * DEG_TO_RAD;
                break;
            case TWILIGHT_CIVIL:
                tmp = -6 * DEG_TO_RAD;
                break;
            case TWILIGHT_NAUTICAL:
                tmp = -12 * DEG_TO_RAD;
                break;
            case TWILIGHT_ASTRONOMICAL:
                tmp = -18 * DEG_TO_RAD;
                break;
        }

        // Compute cosine of hour angle
        tmp = (Math.sin(tmp) - Math.sin(longitude) * Math.sin(dec)) / (Math.cos(longitude) * Math.cos(dec));
        double celestialHoursToEarthTime = RAD_TO_DAY / SIDEREAL_DAY_LENGTH;

        // Make calculations for the meridian
        double transit_time1 = celestialHoursToEarthTime * normalizeRadians(ra - lst);
        double transit_time2 = celestialHoursToEarthTime * (normalizeRadians(ra - lst) - TWO_PI);
        double transit_alt   = Math.asin(Math.sin(dec) * Math.sin(longitude) + Math.cos(dec) * Math.cos(longitude));
        if (transit_alt > -3 * DEG_TO_RAD) {
            double r    = 0.016667 * DEG_TO_RAD * Math.abs(Math.tan(PI_OVER_TWO - (transit_alt * RAD_TO_DEG + 7.31 / (transit_alt * RAD_TO_DEG + 4.4)) * DEG_TO_RAD));
            double refr = r * (0.28 * 1010 / (10 + 273.0)); // Assuming pressure of 1010 mb and T = 10 C
            transit_alt = Math.min(transit_alt + refr, PI_OVER_TWO); // This is not accurate, but acceptable
        }

        // Obtain the current event in time
        double transit_time  = transit_time1;
        double jdToday       = Math.floor(jd_UT - 0.5) + 0.5;
        double transitToday2 = Math.floor(jd_UT + transit_time2 - 0.5) + 0.5;
        // Obtain the transit time. Preference should be given to the closest event
        // in time to the current calculation time
        if (jdToday == transitToday2 && Math.abs(transit_time2) < Math.abs(transit_time1)) transit_time = transit_time2;
        double transit = jd_UT + transit_time;

        // Make calculations for rise and set
        double rise = -1, set = -1;
        if (Math.abs(tmp) <= 1.0) {
            double ang_hor    = Math.abs(Math.acos(tmp));
            double rise_time1 = celestialHoursToEarthTime * normalizeRadians(ra - ang_hor - lst);
            double set_time1  = celestialHoursToEarthTime * normalizeRadians(ra + ang_hor - lst);
            double rise_time2 = celestialHoursToEarthTime * (normalizeRadians(ra - ang_hor - lst) - TWO_PI);
            double set_time2  = celestialHoursToEarthTime * (normalizeRadians(ra + ang_hor - lst) - TWO_PI);

            // Obtain the current events in time. Preference should be given to the closest event
            // in time to the current calculation time (so that iteration in other method will converge)
            double rise_time  = rise_time1;
            double riseToday2 = Math.floor(jd_UT + rise_time2 - 0.5) + 0.5;
            if (jdToday == riseToday2 && Math.abs(rise_time2) < Math.abs(rise_time1)) rise_time = rise_time2;

            double set_time  = set_time1;
            double setToday2 = Math.floor(jd_UT + set_time2 - 0.5) + 0.5;
            if (jdToday == setToday2 && Math.abs(set_time2) < Math.abs(set_time1)) set_time = set_time2;
            rise = jd_UT + rise_time;
            set = jd_UT + set_time;
        }

        return new double[]{azi, alt, rise, set, transit, transit_alt, ra, dec, dist};
    }
Calculates everything for the Sun and the Moon.
src/main/java/eu/hansolo/tilesfx/tools/SunMoonCalculator.java