<?php
/*
 * Part of the VastuIQ plugin: the Vastu scoring engine, the same rules VastuIQ.com uses.
 * License: GPLv2 or later.
 */
/**
 * VastuIQ Home Report Engine — magnetic declination (true north vs compass north)
 *
 * A direct PHP port of NOAA's World Magnetic Model (WMM) legacy C code,
 * using the official WMM-2025 coefficients (data/wmm-2025.json, public
 * domain, NOAA NCEI / British Geological Survey), valid 2025.0–2030.0.
 * Ported from the pygeomag reference implementation and checked against it.
 *
 * Declination D = angle from TRUE north to MAGNETIC north, east positive.
 *   true bearing = magnetic bearing + D
 */

namespace VastuIQWP\Engine;

defined( 'ABSPATH' ) || exit;

final class MagneticModel
{
    private const MAXORD = 12;
    private float $epoch;
    private string $model;
    private array $c = [];
    private array $cd = [];
    private array $snorm = [];
    private array $k = [];
    private array $fn = [];
    private array $fm = [];

    public function __construct(?string $cofFile = null)
    {
        $cofFile = $cofFile ?? dirname(__DIR__) . '/data/wmm-2025.json';
        $wmm = json_decode(is_readable($cofFile) ? (string)file_get_contents($cofFile) : '', true); // phpcs:ignore WordPress.WP.AlternativeFunctions.file_get_contents_file_get_contents -- reads a file bundled with the plugin.
        $lines = is_array($wmm) && isset($wmm['header'], $wmm['rows']) ? array_merge([(string)$wmm['header']], array_map(fn($r) => implode(' ', (array)$r), (array)$wmm['rows']), ['9999']) : [];
        if (!$lines) {
            throw new \RuntimeException('Magnetic model file missing'); // phpcs:ignore WordPress.Security.EscapeOutput.ExceptionNotEscaped -- message is never printed.
        }
        $head = preg_split('/\s+/', trim(array_shift($lines)));
        if (count($head) !== 3) {
            throw new \RuntimeException('Invalid magnetic model header'); // phpcs:ignore WordPress.Security.EscapeOutput.ExceptionNotEscaped -- message is never printed.
        }
        $this->epoch = (float)$head[0];
        $this->model = $head[1];

        $size = self::MAXORD + 1;
        $c = $cd = $k = array_fill(0, $size, array_fill(0, $size, 0.0));
        $snorm = array_fill(0, $size * $size, 0.0);
        $fn = $fm = array_fill(0, $size, 0.0);

        foreach ($lines as $line) {
            if (strncmp($line, '9999', 4) === 0) break;
            $v = preg_split('/\s+/', trim($line));
            if (count($v) !== 6) throw new \RuntimeException('Corrupt magnetic model record'); // phpcs:ignore WordPress.Security.EscapeOutput.ExceptionNotEscaped -- message is never printed.
            [$n, $m] = [(int)$v[0], (int)$v[1]];
            if ($m > self::MAXORD) break;
            $c[$m][$n] = (float)$v[2];
            $cd[$m][$n] = (float)$v[4];
            if ($m !== 0) {
                $c[$n][$m - 1] = (float)$v[3];
                $cd[$n][$m - 1] = (float)$v[5];
            }
        }

        // Schmidt semi-normalisation (same steps as the NOAA C code)
        $snorm[0] = 1.0;
        for ($n = 1; $n <= self::MAXORD; $n++) {
            $snorm[$n] = $snorm[$n - 1] * (2 * $n - 1) / $n;
            $j = 2;
            for ($m = 0; $m <= $n; $m++) {
                $k[$m][$n] = (($n - 1) * ($n - 1) - $m * $m) / ((2 * $n - 1) * (2 * $n - 3));
                if ($m > 0) {
                    $flnmj = (($n - $m + 1) * $j) / ($n + $m);
                    $snorm[$n + $m * $size] = $snorm[$n + ($m - 1) * $size] * sqrt($flnmj);
                    $j = 1;
                    $c[$n][$m - 1] *= $snorm[$n + $m * $size];
                    $cd[$n][$m - 1] *= $snorm[$n + $m * $size];
                }
                $c[$m][$n] *= $snorm[$n + $m * $size];
                $cd[$m][$n] *= $snorm[$n + $m * $size];
            }
            $fn[$n] = (float)($n + 1);
            $fm[$n] = (float)$n;
        }
        $k[1][1] = 0.0;

        [$this->c, $this->cd, $this->snorm, $this->k, $this->fn, $this->fm] = [$c, $cd, $snorm, $k, $fn, $fm];
    }

    public function modelName(): string { return $this->model; }
    public function validFrom(): float { return $this->epoch; }
    public function validTo(): float { return $this->epoch + 5.0; }

    /** Decimal year for a timestamp, e.g. 2026.79. */
    public static function decimalYear(?int $ts = null): float
    {
        $ts = $ts ?? time();
        $y = (int)gmdate('Y', $ts);
        $start = gmmktime(0, 0, 0, 1, 1, $y);
        $end = gmmktime(0, 0, 0, 1, 1, $y + 1);
        return $y + ($ts - $start) / ($end - $start);
    }

    /**
     * Declination in degrees (east positive) at latitude/longitude (degrees),
     * altitude in km above the WGS84 ellipsoid, and decimal year.
     */
    public function declination(float $lat, float $lon, float $altKm = 0.0, ?float $year = null): float
    {
        $year = $year ?? self::decimalYear();
        // Outside the model's life span, clamp to it (the error grows slowly; reports show the year used).
        $year = max($this->epoch, min($this->epoch + 4.999, $year));
        $dt = $year - $this->epoch;

        $size = self::MAXORD + 1;
        $p = $this->snorm;
        $tc = $dp = array_fill(0, $size, array_fill(0, $size, 0.0));
        $sp = $cp = $pp = array_fill(0, $size, 0.0);
        $sp[0] = 0.0; $cp[0] = $pp[0] = 1.0; $dp[0][0] = 0.0;

        $a = 6378.137; $b = 6356.7523142; $re = 6371.2;
        $a2 = $a * $a; $b2 = $b * $b; $c2 = $a2 - $b2;
        $a4 = $a2 * $a2; $b4 = $b2 * $b2; $c4 = $a4 - $b4;

        $rlon = deg2rad($lon); $rlat = deg2rad($lat);
        $srlon = sin($rlon); $srlat = sin($rlat);
        $crlon = cos($rlon); $crlat = cos($rlat);
        $srlat2 = $srlat * $srlat; $crlat2 = $crlat * $crlat;
        $sp[1] = $srlon; $cp[1] = $crlon;

        // geodetic -> spherical
        $q = sqrt($a2 - $c2 * $srlat2);
        $q1 = $altKm * $q;
        $q2 = (($q1 + $a2) / ($q1 + $b2)) ** 2;
        $ct = $srlat / sqrt($q2 * $crlat2 + $srlat2);
        $st = sqrt(1.0 - $ct * $ct);
        $r2 = $altKm * $altKm + 2.0 * $q1 + ($a4 - $c4 * $srlat2) / ($q * $q);
        $r = sqrt($r2);
        $d = sqrt($a2 * $crlat2 + $b2 * $srlat2);
        $ca = ($altKm + $d) / $r;
        $sa = $c2 * $crlat * $srlat / ($r * $d);

        for ($m = 2; $m <= self::MAXORD; $m++) {
            $sp[$m] = $sp[1] * $cp[$m - 1] + $cp[1] * $sp[$m - 1];
            $cp[$m] = $cp[1] * $cp[$m - 1] - $sp[1] * $sp[$m - 1];
        }

        $aor = $re / $r;
        $ar = $aor * $aor;
        $br = $bt = $bp = $bpp = 0.0;

        for ($n = 1; $n <= self::MAXORD; $n++) {
            $ar *= $aor;
            for ($m = 0; $m <= $n; $m++) {
                if ($n === $m) {
                    $p[$n + $m * $size] = $st * $p[$n - 1 + ($m - 1) * $size];
                    $dp[$m][$n] = $st * $dp[$m - 1][$n - 1] + $ct * $p[$n - 1 + ($m - 1) * $size];
                } elseif ($n === 1 && $m === 0) {
                    $p[$n + $m * $size] = $ct * $p[$n - 1 + $m * $size];
                    $dp[$m][$n] = $ct * $dp[$m][$n - 1] - $st * $p[$n - 1 + $m * $size];
                } elseif ($n > 1 && $n !== $m) {
                    if ($m > $n - 2) {
                        $p[$n - 2 + $m * $size] = 0.0;
                        $dp[$m][$n - 2] = 0.0;
                    }
                    $p[$n + $m * $size] = $ct * $p[$n - 1 + $m * $size] - $this->k[$m][$n] * $p[$n - 2 + $m * $size];
                    $dp[$m][$n] = $ct * $dp[$m][$n - 1] - $st * $p[$n - 1 + $m * $size] - $this->k[$m][$n] * $dp[$m][$n - 2];
                }

                $tc[$m][$n] = $this->c[$m][$n] + $dt * $this->cd[$m][$n];
                if ($m !== 0) {
                    $tc[$n][$m - 1] = $this->c[$n][$m - 1] + $dt * $this->cd[$n][$m - 1];
                }

                $par = $ar * $p[$n + $m * $size];
                if ($m === 0) {
                    $temp1 = $tc[$m][$n] * $cp[$m];
                    $temp2 = $tc[$m][$n] * $sp[$m];
                } else {
                    $temp1 = $tc[$m][$n] * $cp[$m] + $tc[$n][$m - 1] * $sp[$m];
                    $temp2 = $tc[$m][$n] * $sp[$m] - $tc[$n][$m - 1] * $cp[$m];
                }
                $bt -= $ar * $temp1 * $dp[$m][$n];
                $bp += $this->fm[$m] * $temp2 * $par;
                $br += $this->fn[$n] * $temp1 * $par;

                if ($st == 0.0 && $m === 1) {
                    $pp[$n] = ($n === 1) ? $pp[$n - 1] : $ct * $pp[$n - 1] - $this->k[$m][$n] * $pp[$n - 2];
                    $bpp += $this->fm[$m] * $temp2 * ($ar * $pp[$n]);
                }
            }
        }

        $bp = ($st == 0.0) ? $bpp : $bp / $st;
        $bx = -$bt * $ca - $br * $sa;
        $by = $bp;
        return rad2deg(atan2($by, $bx));
    }
}