package de.arxsilex.maptools.converter.core

import kotlin.math.PI
import kotlin.math.asin
import kotlin.math.cos
import kotlin.math.acos
import kotlin.math.max
import kotlin.math.min
import kotlin.math.sin
import kotlin.math.sqrt

data class GeoRectangleMetrics(
    val centerLatitude: Double,
    val centerLongitude: Double,
    val widthKilometers: Double,
    val heightKilometers: Double,
    val areaSquareKilometers: Double,
)

data class GeoRegionMetrics(
    val centerLatitude: Double,
    val centerLongitude: Double,
    val widthKilometers: Double,
    val heightKilometers: Double,
    val areaSquareKilometers: Double,
)

fun GeoBounds.rectangleMetrics(): GeoRectangleMetrics {
    val centerLatitude = (south + north) / 2.0
    val centerLongitude = (west + east) / 2.0
    return GeoRectangleMetrics(
        centerLatitude = centerLatitude,
        centerLongitude = centerLongitude,
        widthKilometers = greatCircleKilometers(centerLatitude, west, centerLatitude, east),
        heightKilometers = greatCircleKilometers(south, centerLongitude, north, centerLongitude),
        areaSquareKilometers = sphericalRectangleAreaSquareKilometers(west, south, east, north),
    )
}

fun GeoSelectionRegion.metrics(): GeoRegionMetrics = when (this) {
    is GeoRectangleRegion -> bounds.rectangleMetrics().let {
        GeoRegionMetrics(it.centerLatitude, it.centerLongitude, it.widthKilometers, it.heightKilometers, it.areaSquareKilometers)
    }
    is GeoCircleRegion -> GeoRegionMetrics(
        centerLatitude = centerLatitude,
        centerLongitude = centerLongitude,
        widthKilometers = radiusKilometers * 2.0,
        heightKilometers = radiusKilometers * 2.0,
        areaSquareKilometers = 2.0 * PI * EARTH_MEAN_RADIUS_KM * EARTH_MEAN_RADIUS_KM *
            (1.0 - cos(radiusKilometers / EARTH_MEAN_RADIUS_KM)),
    )
}

/**
 * Union area for rectangles and circles. Rectangle-only selections use the exact method below;
 * mixed selections use deterministic spherical strip integration and merge every strip before
 * adding it, so overlaps are never counted twice.
 */
fun Iterable<GeoSelectionRegion>.unionAreaSquareKilometers(clipTo: GeoBounds? = null): Double {
    val regions = toList()
    if (regions.isEmpty()) return 0.0
    if (regions.all { it is GeoRectangleRegion } && clipTo == null) {
        return regions.map { it.bounds }.unionAreaSquareKilometers()
    }
    val rawWest = maxOf(regions.minOf { it.bounds.west }, clipTo?.west ?: -180.0)
    val rawEast = minOf(regions.maxOf { it.bounds.east }, clipTo?.east ?: 180.0)
    if (rawEast <= rawWest) return 0.0
    val slices = 2_048
    val longitudeStep = (rawEast - rawWest) / slices
    var total = 0.0
    repeat(slices) { index ->
        val west = rawWest + longitudeStep * index
        val east = west + longitudeStep
        val longitude = (west + east) / 2.0
        val intervals = regions.mapNotNull { region ->
            region.latitudeIntervalAt(longitude)?.let { (south, north) ->
                maxOf(south, clipTo?.south ?: -90.0) to minOf(north, clipTo?.north ?: 90.0)
            }?.takeIf { it.second > it.first }
        }.sortedBy { it.first }
        if (intervals.isEmpty()) return@repeat
        var south = intervals.first().first
        var north = intervals.first().second
        intervals.drop(1).forEach { (nextSouth, nextNorth) ->
            if (nextSouth > north) {
                total += sphericalRectangleAreaSquareKilometers(west, south, east, north)
                south = nextSouth
                north = nextNorth
            } else {
                north = maxOf(north, nextNorth)
            }
        }
        total += sphericalRectangleAreaSquareKilometers(west, south, east, north)
    }
    return total
}

private fun GeoSelectionRegion.latitudeIntervalAt(longitude: Double): Pair<Double, Double>? = when (this) {
    is GeoRectangleRegion -> if (longitude in bounds.west..bounds.east) bounds.south to bounds.north else null
    is GeoCircleRegion -> {
        if (longitude !in bounds.west..bounds.east) return null
        val horizontal = greatCircleKilometers(centerLatitude, centerLongitude, centerLatitude, longitude)
        if (horizontal > radiusKilometers) return null
        val latitudeDelta = sqrt(max(0.0, radiusKilometers * radiusKilometers - horizontal * horizontal)) /
            EARTH_MEAN_RADIUS_KM * 180.0 / PI
        (centerLatitude - latitudeDelta).coerceAtLeast(-90.0) to
            (centerLatitude + latitudeDelta).coerceAtMost(90.0)
    }
}

/** Exact union area for axis-aligned geographic rectangles on the mean-earth sphere. */
fun Iterable<GeoBounds>.unionAreaSquareKilometers(): Double {
    val rectangles = toList()
    if (rectangles.isEmpty()) return 0.0
    val longitudeEdges = rectangles.flatMap { listOf(it.west, it.east) }.distinct().sorted()
    var result = 0.0
    for (index in 0 until longitudeEdges.lastIndex) {
        val west = longitudeEdges[index]
        val east = longitudeEdges[index + 1]
        if (east <= west) continue
        val latitudeIntervals = rectangles.asSequence()
            .filter { it.west < east && it.east > west }
            .map { it.south to it.north }
            .sortedBy(Pair<Double, Double>::first)
            .toList()
        if (latitudeIntervals.isEmpty()) continue
        var south = latitudeIntervals.first().first
        var north = latitudeIntervals.first().second
        latitudeIntervals.drop(1).forEach { (nextSouth, nextNorth) ->
            if (nextSouth > north) {
                result += sphericalRectangleAreaSquareKilometers(west, south, east, north)
                south = nextSouth
                north = nextNorth
            } else {
                north = maxOf(north, nextNorth)
            }
        }
        result += sphericalRectangleAreaSquareKilometers(west, south, east, north)
    }
    return result
}

fun greatCircleKilometers(
    latitudeA: Double,
    longitudeA: Double,
    latitudeB: Double,
    longitudeB: Double,
): Double {
    val latitudeDelta = radians(latitudeB - latitudeA)
    val longitudeDelta = radians(longitudeB - longitudeA)
    val a = sin(latitudeDelta / 2.0) * sin(latitudeDelta / 2.0) +
        cos(radians(latitudeA)) * cos(radians(latitudeB)) *
        sin(longitudeDelta / 2.0) * sin(longitudeDelta / 2.0)
    return 2.0 * EARTH_MEAN_RADIUS_KM * asin(sqrt(a.coerceIn(0.0, 1.0)))
}

private fun sphericalRectangleAreaSquareKilometers(
    west: Double,
    south: Double,
    east: Double,
    north: Double,
): Double = EARTH_MEAN_RADIUS_KM * EARTH_MEAN_RADIUS_KM *
    radians(east - west) * (sin(radians(north)) - sin(radians(south)))

private fun radians(degrees: Double): Double = degrees * PI / 180.0
