<?php
/*
 * Part of the VastuIQ plugin: the Vastu scoring engine, the same rules VastuIQ.com uses.
 * License: GPLv2 or later.
 */
/**
 * VastuIQ Home Report Engine — sun path for a location
 *
 * NOAA Solar Calculator equations (Global Monitoring Laboratory, based on
 * Meeus, "Astronomical Algorithms"). Accurate to well under a degree for
 * the years we report on — far finer than a 22.5° Vastu zone.
 *
 * All bearings are TRUE bearings in degrees clockwise from true north.
 */

namespace VastuIQWP\Engine;

defined( 'ABSPATH' ) || exit;

final class SunPath
{
    public const TROPIC = 23.44;   // latitude of the tropics (Cancer / Capricorn)

    /** Sun position at a UTC timestamp: ['azimuth' => deg, 'elevation' => deg]. */
    public static function position(float $lat, float $lon, int $ts): array
    {
        [$dec, $eot] = self::declinationAndEot($ts);
        $minutes = ($ts % 86400) / 60.0;
        $tst = fmod($minutes + $eot + 4.0 * $lon, 1440.0);
        if ($tst < 0) $tst += 1440.0;
        $ha = $tst / 4.0 - 180.0;

        $latR = deg2rad($lat); $decR = deg2rad($dec); $haR = deg2rad($ha);
        $cosZ = sin($latR) * sin($decR) + cos($latR) * cos($decR) * cos($haR);
        $cosZ = max(-1.0, min(1.0, $cosZ));
        $zen = acos($cosZ);
        $elev = 90.0 - rad2deg($zen);

        $den = cos($latR) * sin($zen);
        if (abs($den) < 1e-9) {
            $az = $lat > 0 ? 180.0 : 0.0;
        } else {
            $x = (sin($latR) * cos($zen) - sin($decR)) / $den;
            $x = max(-1.0, min(1.0, $x));
            $a = rad2deg(acos($x));
            $az = $ha > 0 ? fmod($a + 180.0, 360.0) : fmod(540.0 - $a, 360.0);
        }
        if ($az >= 360.0) $az -= 360.0;
        return ['azimuth' => $az, 'elevation' => $elev + self::refraction($elev)];
    }

    /** Solar declination (deg) and equation of time (minutes) at a UTC timestamp. */
    public static function declinationAndEot(int $ts): array
    {
        $jd = $ts / 86400.0 + 2440587.5;
        $t = ($jd - 2451545.0) / 36525.0;
        $l0 = fmod(280.46646 + $t * (36000.76983 + $t * 0.0003032), 360.0);
        $m = 357.52911 + $t * (35999.05029 - 0.0001537 * $t);
        $e = 0.016708634 - $t * (0.000042037 + 0.0000001267 * $t);
        $mR = deg2rad($m);
        $c = sin($mR) * (1.914602 - $t * (0.004817 + 0.000014 * $t))
           + sin(2 * $mR) * (0.019993 - 0.000101 * $t)
           + sin(3 * $mR) * 0.000289;
        $omega = deg2rad(125.04 - 1934.136 * $t);
        $lambda = deg2rad($l0 + $c - 0.00569 - 0.00478 * sin($omega));
        $eps0 = 23 + (26 + (21.448 - $t * (46.815 + $t * (0.00059 - $t * 0.001813))) / 60) / 60;
        $eps = deg2rad($eps0 + 0.00256 * cos($omega));
        $dec = rad2deg(asin(sin($eps) * sin($lambda)));

        $y = tan($eps / 2) ** 2;
        $l0R = deg2rad($l0);
        $eot = 4 * rad2deg(
            $y * sin(2 * $l0R) - 2 * $e * sin($mR) + 4 * $e * $y * sin($mR) * cos(2 * $l0R)
            - 0.5 * $y * $y * sin(4 * $l0R) - 1.25 * $e * $e * sin(2 * $mR)
        );
        return [$dec, $eot];
    }

    /** Atmospheric refraction correction (deg), NOAA approximation. */
    private static function refraction(float $elev): float
    {
        if ($elev > 85.0) return 0.0;
        $te = tan(deg2rad($elev));
        if ($elev > 5.0) $r = 58.1 / $te - 0.07 / $te ** 3 + 0.000086 / $te ** 5;
        elseif ($elev > -0.575) $r = 1735.0 + $elev * (-518.2 + $elev * (103.4 + $elev * (-12.79 + $elev * 0.711)));
        else $r = -20.774 / $te;
        return $r / 3600.0;
    }

    /**
     * Sunrise, solar noon and sunset for a UTC calendar date.
     * Returns null fields during polar day/night (not expected for our markets).
     */
    public static function day(float $lat, float $lon, int $year, int $month, int $dayOfMonth): array
    {
        $midnight = gmmktime(0, 0, 0, $month, $dayOfMonth, $year);
        // solar noon (UTC minutes), refined once
        $noon = 720.0 - 4.0 * $lon;
        for ($i = 0; $i < 2; $i++) {
            [, $eot] = self::declinationAndEot((int)round($midnight + $noon * 60));
            $noon = 720.0 - 4.0 * $lon - $eot;
        }
        $noonTs = (int)round($midnight + $noon * 60);
        $noonPos = self::position($lat, $lon, $noonTs);

        $events = []; $polar = null;
        foreach (['sunrise' => -1, 'sunset' => 1] as $name => $sign) {
            $tMin = $noon;
            $ok = true;
            for ($i = 0; $i < 3; $i++) {
                [$dec, $eot] = self::declinationAndEot((int)round($midnight + $tMin * 60));
                $latR = deg2rad($lat); $decR = deg2rad($dec);
                $cosH = cos(deg2rad(90.833)) / (cos($latR) * cos($decR)) - tan($latR) * tan($decR);
                if ($cosH < -1 || $cosH > 1) { $ok = false; $polar = $cosH < -1 ? 'day' : 'night'; break; }   // midnight sun / polar night (far north of Sweden)
                $h = rad2deg(acos($cosH));
                $tMin = 720.0 - 4.0 * ($lon - $sign * $h) - $eot;
            }
            $events[$name] = $ok ? (int)round($midnight + $tMin * 60) : null;
        }

        $rise = $events['sunrise']; $set = $events['sunset'];
        return [
            'date' => sprintf('%04d-%02d-%02d', $year, $month, $dayOfMonth),
            'sunrise_azimuth' => $rise ? self::position($lat, $lon, $rise)['azimuth'] : null,
            'sunset_azimuth' => $set ? self::position($lat, $lon, $set)['azimuth'] : null,
            'noon_elevation' => $noonPos['elevation'],
            'noon_azimuth' => $noonPos['azimuth'],
            'noon_side' => ($noonPos['azimuth'] > 90 && $noonPos['azimuth'] < 270) ? 'S' : 'N',
            'day_hours' => ($rise && $set) ? ($set - $rise) / 3600.0 : ($polar === 'day' ? 24.0 : ($polar === 'night' ? 0.0 : null)),
        ];
    }

    /**
     * Everything the report shows about the sun for one location:
     * the two solstices and the March equinox, plus plain-language facts.
     */
    public static function summary(float $lat, float $lon, ?int $year = null): array
    {
        $year = $year ?? (int)gmdate('Y');
        $jun = self::day($lat, $lon, $year, 6, 21);
        $dec = self::day($lat, $lon, $year, 12, 21);
        $mar = self::day($lat, $lon, $year, 3, 20);

        $south = $lat < 0;
        $tropical = abs($lat) < self::TROPIC;
        $winter = $south ? $jun : $dec;     // local winter solstice
        $summer = $south ? $dec : $jun;
        $equatorSide = $south ? 'N' : 'S';  // side the sun is on at midday for most of the year
        $poleSide = $south ? 'S' : 'N';

        // Days in the year when the midday sun is on the POLE side (only in the tropics)
        $poleSideDays = 0;
        if ($tropical) {
            for ($d = 0; $d < 365; $d += 2) {
                $ts = gmmktime(12, 0, 0, 1, 1, $year) + $d * 86400;
                [$decl] = self::declinationAndEot($ts);
                if (($south && $decl < $lat) || (!$south && $decl > $lat)) $poleSideDays += 2;
            }
        }

        return [
            'latitude' => $lat,
            'longitude' => $lon,
            'hemisphere' => $south ? 'southern' : 'northern',
            'tropical' => $tropical,
            'equator_side' => $equatorSide,
            'pole_side' => $poleSide,
            'winter_solstice' => $winter,
            'summer_solstice' => $summer,
            'equinox' => $mar,
            'pole_side_noon_days' => $poleSideDays,
            'sunrise_range' => self::range([$jun['sunrise_azimuth'], $dec['sunrise_azimuth'], $mar['sunrise_azimuth']]),
            'sunset_range' => self::range([$jun['sunset_azimuth'], $dec['sunset_azimuth'], $mar['sunset_azimuth']]),
        ];
    }

    /** Smallest and largest of the azimuths that exist (none on a polar day or night). */
    private static function range(array $az): array
    {
        $az = array_values(array_filter($az, fn($v) => $v !== null));
        return $az ? [min($az), max($az)] : [90.0, 90.0];
    }
}