/***********************************************/
/*           Copyright (c) 2025 Belmu          */
/*             All Rights Reserved             */
/***********************************************/

/*
    [Credits]:
        Jessie - providing Minecraft date to milliseconds conversion equation (https://github.com/Jessie-LC)

    [References]:
        Museo Galileo. (2015). In-depth Ecliptic coordinates. https://catalogue.museogalileo.it/indepth/EclipticCoordinates.html#:~:text=The%20first%20of%20the%20two,celestial%20body%20and%20the%20ecliptic.
        Strous, L. (2021). Astronomy Answers Positions in the Sky. https://aa.quae.nl/en/reken/hemelpositie.html
        Strous, L. (2021). Astronomy Answers Position of the Sun. https://aa.quae.nl/en/reken/zonpositie.html
        Wikipedia. (2024). Equation of the center. https://en.wikipedia.org/wiki/Equation_of_the_center
        Wikipedia. (2024). Mean anomaly. https://en.wikipedia.org/wiki/Mean_anomaly
*/

const float ms_in_day = 8.64e7; // Milliseconds in a day

const float J2000 = 2451545.0;   // January 1st 2000 (used as a reference in astronomy) [Julians]
const float M0_E  = 357.5291;    // Earth's mean anomaly at epoch J2000 [deg]
const float M1_E  = 0.98560028;  // Rate of change of the Earth's mean anomaly in degrees per day [deg/day]
const float st0_E = 280.1470;    // Sidereal time of the Earth at longitude 0° at the instant defined by J2000. [deg]
const float st1_E = 360.9856235; // Rate of change of the Earth's sidereal time in degrees per day [deg/day]
const float P     = 102.9373;    // Ecliptic longitude of the Earth's perihelion (point in the orbit closest to the Sun) [deg]
const float obl   = 23.4393;     // Obliquity of Earth's equator (inclination of the equator with respect to the ecliptic) [deg]

const float obl_rad = obl * deg_to_rad;

float minecraftDateToScaledMillisecondsDate() {
    /* 
        24'000 ticks per Minecraft day
        50 milliseconds per tick
        72 Minecraft days per real day
        6000 ticks offset to line up with Minecraft noon
    */
    return worldDay * 24000.0 + (worldTime + 6000.0) * 50.0 * 72.0;
}

float minecraftDaysToJulianDays() {
    /*
        The -0.5 comes from the fact that Julian days start at noon.
        There are 72 Minecraft days in a real day, hence the division by 72.
    */
    return (minecraftDateToScaledMillisecondsDate() / ms_in_day) + J2000 - 0.5;
}

float julianDaysToGregorianDays(float julianDays) {
    return julianDays - J2000;
}

/*
    The mean anomaly is the angular distance from the pericenter which a body would have if it moved in a circular
    orbit, with constant speed, with the same orbital period as the body in its elliptical orbit.
*/
float meanAnomaly(float days, float M0, float M1) {
    return M0 + M1 * days;
}

/*
    Kepler's equation of the center, cannot be solved, only approximated.
    Angular difference between the actual position of a body in its elliptical orbit and the position it would
    occupy if its motion were uniform, in a circular orbit of the same period (hence why the mean anomaly is used).
*/
float keplerEquationApproximation(float M) {
    return 1.9148 * sin(M) + 0.02 * sin(2.0 * M) + 0.0003 * sin(3.0 * M);
}

/*
    The ecliptic longitude is the angle computed from west to east from the vernal equinox to the intersection between 
    the ecliptic meridian of the celestial body and the ecliptic.
*/
float eclipticLongitudeEarth(float M) {
    float C = keplerEquationApproximation(deg_to_rad * M);
    return deg_to_rad * (M + P + C) + PI;
}

// In the equatorial coordinate system, the right ascension is how far the body is from the vernal equinox. [deg]
float rightAscension(float longitude, float latitude) {
    return atan(sin(longitude) * cos(obl_rad) - tan(latitude) * sin(obl_rad), cos(longitude));
}

// In the equatorial coordinate system, the declination is how far the body is from the celestial equator. [deg]
float declination(float longitude, float latitude) {
    return asin(sin(latitude) * cos(obl_rad) + cos(latitude) * sin(obl_rad) * sin(longitude));
}

float azimuth(float hourAngle, float latitude, float declination) {
    return atan(sin(hourAngle), cos(hourAngle) * sin(latitude) - tan(declination) * cos(latitude));
}

float altitude(float hourAngle, float latitude, float declination) {
    return asin(sin(latitude) * sin(declination) + cos(latitude) * cos(declination) * cos(hourAngle));
}

// Sidereal time measures how the rotation of a celestial body aligns with the stars, rather than its system's star.
float siderealTime(float days, float longitudeWest, float st0, float st1) {
    return st0 + st1 * days - longitudeWest;
}

vec3 calculateSunPosition(float latitude, float longitude) {
    float days = minecraftDaysToJulianDays();
          days = julianDaysToGregorianDays(days);

    float meanAnomalyDeg    = meanAnomaly(days, M0_E, M1_E);
    float eclipticLongitude = eclipticLongitudeEarth(meanAnomalyDeg);

    // As we are using the ecliptic plane as a reference, the latitude is always 0°.
    vec2 equatorialCoords;
    equatorialCoords.x = rightAscension(eclipticLongitude, 0.0);
    equatorialCoords.y = declination   (eclipticLongitude, 0.0);

    // Using the West longitude instead of the East longitude is a choice made for facility.
    float longitudeWest   = -longitude;
    float siderealTimeRad = deg_to_rad * siderealTime(days, longitudeWest, st0_E, st1_E);

    // Hour angle, indicates how long ago the body passed through the celestial meridian in sidereal time.
    float hourAngle   = siderealTimeRad - equatorialCoords.x;
    float latitudeRad = deg_to_rad * latitude;

    // The sun is always positioned on the ecliptic plane, so the Z component is zero.
    vec3 horizontalCoordinates;
    horizontalCoordinates.x = azimuth (hourAngle, latitudeRad, equatorialCoords.y);
    horizontalCoordinates.y = altitude(hourAngle, latitudeRad, equatorialCoords.y);
    horizontalCoordinates.z = 0.0;

    return horizontalCoordinates;
}

vec3 calculateMoonPosition(float latitude, float longitude) {
    float days = minecraftDaysToJulianDays();
          days = julianDaysToGregorianDays(days);

    /*
        Very simplified equations, because the Moon's orbit actually rotates with a period of around 18 years,
        making everything much more complicated. Those three values are in radians.
    */
    float meanAnomaly           = deg_to_rad * (134.963 + 13.064993 * days);
    float meanEclipticLongitude = deg_to_rad * (218.316 + 13.176396 * days);
    float meanDistance          = deg_to_rad * (93.272  + 13.229350 * days);

    float eclipticLongitude = meanEclipticLongitude + deg_to_rad * 6.289 * sin(meanAnomaly);
    float eclipticLatitude  = deg_to_rad * 5.128 * sin(meanDistance);
    float distanceKm        = 385001.0 - 20905.0 * cos(meanAnomaly);

    vec3 equatorialCoords;
    equatorialCoords.x = rightAscension(eclipticLongitude, eclipticLatitude);
    equatorialCoords.y = declination   (eclipticLongitude, eclipticLatitude);
    equatorialCoords.z = distanceKm;

    // Using the West longitude instead of the East longitude is a choice made for facility.
    float longitudeWest   = -longitude;
    float siderealTimeRad = deg_to_rad * siderealTime(days, longitudeWest, st0_E, st1_E);

    // Hour angle, indicates how long ago the body passed through the celestial meridian in sidereal time.
    float hourAngle   = siderealTimeRad - equatorialCoords.x;
    float latitudeRad = deg_to_rad * latitude;

    vec3 horizontalCoordinates;
    horizontalCoordinates.x = azimuth (hourAngle, latitudeRad, equatorialCoords.y);
    horizontalCoordinates.y = altitude(hourAngle, latitudeRad, equatorialCoords.y);
    horizontalCoordinates.z = equatorialCoords.z;

    return horizontalCoordinates;
}
