import math

from .geo_bounding_box_types import BoundingBox
from .math_atan import atan
from .math_sin_cos import sin_cos

# The latitude/longitude box that contains every point within a distance of a
# centre, by the method of J. P. Matuschek, "Finding Points Within a Distance
# of a Latitude/Longitude Using Bounding Coordinates"
# (http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates).
#
# Sine, cosine and arctangent come from math.sin-cos and math.atan instead of
# math.sin and friends, because the platform maths library may differ in the
# last bit between Python, JavaScript and Rust. The shared helpers use only
# + - * /, and the rest here only + - * / sqrt floor ceil, all correctly
# rounded, in the same order in all three.

# IUGG mean radius R1 of the GRS 80 ellipsoid (Moritz, Journal of Geodesy 74 (2000)).
EARTH_RADIUS_METRES = 6371008.8
# Ten millionths of a degree, about a centimetre.
SCALE = 10000000.0
# Within a millionth of a grid step of a 7-decimal value counts as that value,
# so a centre of 51.5074 stays 51.5074 rather than becoming 51.5073999.
SNAP = 0.000001

PI = 3.141592653589793
HALF_PI = PI / 2
DEGREES = PI / 180


def _asin_positive(x: float) -> float:
    # asin for 0 <= x, on math.atan: asin x = atan(x / sqrt(1 - x^2)), with
    # (1 - x)(1 + x) to keep precision near 1. There is no shared asin; x >= 1
    # (reachable only by rounding, just short of a pole) is a quarter turn.
    if x >= 1:
        return HALF_PI
    return atan(x / math.sqrt((1 - x) * (1 + x)))


# Outward rounding: a lower bound goes down and an upper bound goes up, so
# rounding can only grow the box, never cut the circle.
def _round_down(x: float) -> float:
    return float(math.floor(x * SCALE + SNAP)) / SCALE + 0.0


def _round_up(x: float) -> float:
    return float(math.ceil(x * SCALE - SNAP)) / SCALE + 0.0


def _check_number(value: object, what: str) -> None:
    if isinstance(value, bool) or not isinstance(value, (int, float)) or not math.isfinite(value):
        raise TypeError("%s must be a finite number" % what)


def bounding_box(lat: float, lng: float, distance_metres: float) -> BoundingBox:
    """The box containing every point within distance_metres of (lat, lng)."""
    _check_number(lat, "latitude")
    _check_number(lng, "longitude")
    _check_number(distance_metres, "distance")
    if lat < -90 or lat > 90:
        raise ValueError("latitude must be between -90 and 90 degrees")
    if lng < -180 or lng > 180:
        raise ValueError("longitude must be between -180 and 180 degrees")
    if distance_metres < 0:
        raise ValueError("distance must be 0 or more metres")
    lat, lng, distance_metres = float(lat), float(lng), float(distance_metres)

    r = distance_metres / EARTH_RADIUS_METRES
    r_degrees = r / DEGREES
    min_lat = lat - r_degrees
    max_lat = lat + r_degrees

    if min_lat > -90 and max_lat < 90:
        sin_r = sin_cos(r).sin
        cos_lat = sin_cos(lat * DEGREES).cos
        d_lng = _asin_positive(sin_r / cos_lat) / DEGREES
        min_lng = lng - d_lng
        max_lng = lng + d_lng
        # Past the antimeridian the bound comes round the other side, leaving
        # min_lng > max_lng: the box is the two strips either side of 180.
        if min_lng < -180:
            min_lng += 360
        if max_lng > 180:
            max_lng -= 360
    else:
        # The circle reaches a pole, so it covers every longitude.
        if min_lat < -90:
            min_lat = -90.0
        if max_lat > 90:
            max_lat = 90.0
        min_lng = -180.0
        max_lng = 180.0

    return BoundingBox(
        min_lat=_round_down(min_lat),
        min_lng=_round_down(min_lng),
        max_lat=_round_up(max_lat),
        max_lng=_round_up(max_lng),
    )