RadioPropagationApi/api/app/core/coverage.py

from __future__ import annotations

from typing import Any

from pyproj import Geod

from app.core.antenna import AntennaPattern, gain
from app.core.geo import GeoPoint, sample_path
from app.core.propagation import fspl, itm_loss, p1812_field
from app.models.antenna import AntennaPatternSpec
from app.models.coverage import CoverageRequest
from app.services.vegetation import vegetation_loss_along

_GEOD = Geod(ellps="WGS84")


def _antenna_pattern(spec: AntennaPatternSpec) -> AntennaPattern:
    return AntennaPattern(
        pattern=spec.pattern,
        azimuth_deg=spec.azimuth_deg,
        tilt_deg=spec.tilt_deg,
        gain_dbi=spec.gain_dbi,
        beamwidth_h=spec.beamwidth_h,
        beamwidth_v=spec.beamwidth_v,
        front_to_back_db=spec.front_to_back_db,
        sidelobe_floor_db=spec.sidelobe_floor_db,
        pattern_file=spec.pattern_file,
    )


def _path_loss_db(
    request: CoverageRequest,
    tx: GeoPoint,
    rx: GeoPoint,
    distance_km: float,
    azimuth_deg: float,
    elevation_profile: list[float] | None = None,
) -> float:
    freq_mhz = request.tx.frequency_mhz
    if request.model == "fspl":
        return fspl(freq_mhz, distance_km)
    if request.model == "itm":
        return itm_loss(
            tx,
            rx,
            tx_height_agl=request.tx.height_agl,
            rx_height_agl=request.rx.height_agl,
            freq_mhz=freq_mhz,
            elevation_profile_m=elevation_profile,
        )
    if request.model == "p1812":
        return p1812_field(
            tx,
            rx,
            tx_height_agl=request.tx.height_agl,
            rx_height_agl=request.rx.height_agl,
            freq_mhz=freq_mhz,
            environment=request.environment,
            elevation_profile_m=elevation_profile,
        )
    raise ValueError(f"Unsupported coverage model: {request.model}")


def _rx_power_dbm(
    request: CoverageRequest,
    distance_km: float,
    azimuth_deg: float,
    elevation_profile: list[float] | None = None,
) -> float:
    tx = GeoPoint(lat=request.tx.lat, lon=request.tx.lon)
    rx_point = _GEOD.fwd(request.tx.lon, request.tx.lat, azimuth_deg, distance_km * 1000)
    rx = GeoPoint(lat=rx_point[1], lon=rx_point[0])
    pattern = _antenna_pattern(request.antenna)
    antenna_gain = gain(pattern, azimuth_deg, 0.0)
    path_loss = _path_loss_db(request, tx, rx, distance_km, azimuth_deg, elevation_profile)
    vegetation_db = 0.0
    if request.include_vegetation:
        points = sample_path(tx, rx, max(2, min(64, int(distance_km * 1000 / 250) + 2)))
        vegetation_db = vegetation_loss_along(
            points,
            freq_hz=request.tx.frequency_mhz * 1_000_000,
            include_vegetation=True,
        )
    eirp = request.tx.power_dbm + antenna_gain
    return eirp - path_loss - vegetation_db + request.rx.gain_dbi


def _contour_points(
    request: CoverageRequest,
    level_dbm: float,
    dem_path: str | None = None,
) -> list[list[float]]:
    coords: list[list[float]] = []
    azimuth = 0.0
    while azimuth < 360.0:
        last_good: list[float] | None = None
        distance_m = request.range_step_m
        while distance_m <= request.radius_m:
            distance_km = distance_m / 1000.0
            elevation_profile = None
            if request.model in {"itm", "p1812"} and dem_path is not None:
                tx = GeoPoint(lat=request.tx.lat, lon=request.tx.lon)
                rx_lon, rx_lat, _ = _GEOD.fwd(
                    request.tx.lon, request.tx.lat, azimuth, distance_m
                )
                rx = GeoPoint(lat=rx_lat, lon=rx_lon)
                try:
                    from app.core.dem import elevations_along

                    points = sample_path(tx, rx, 64)
                    elevation_profile = elevations_along(points, dem_path=dem_path).tolist()
                except Exception:
                    elevation_profile = None
            rx_power = _rx_power_dbm(request, distance_km, azimuth, elevation_profile)
            if rx_power >= level_dbm:
                lon, lat, _ = _GEOD.fwd(request.tx.lon, request.tx.lat, azimuth, distance_m)
                last_good = [lon, lat]
                distance_m += request.range_step_m
                continue
            break
        if last_good is not None:
            coords.append(last_good)
        azimuth += request.azimuth_step_deg
    if coords and coords[0] != coords[-1]:
        coords.append(coords[0])
    return coords


def compute_coverage(
    request: CoverageRequest,
    *,
    dem_path: str | None = None,
) -> dict[str, Any]:
    features: list[dict[str, Any]] = []
    for level in request.levels_dbm:
        ring = _contour_points(request, level, dem_path=dem_path)
        if len(ring) < 4:
            continue
        features.append(
            {
                "type": "Feature",
                "properties": {
                    "level_dbm": level,
                    "model": request.model,
                    "frequency_mhz": request.tx.frequency_mhz,
                },
                "geometry": {"type": "Polygon", "coordinates": [ring]},
            }
        )

    if request.format == "geojson":
        return {
            "type": "FeatureCollection",
            "features": features,
            "properties": {
                "model": request.model,
                "radius_m": request.radius_m,
                "levels_dbm": request.levels_dbm,
            },
        }

    return {
        "format": request.format,
        "feature_count": len(features),
        "message": "Raster coverage export is not implemented; use format=geojson",
    }