rfcp/backend/app/services/los_service.py

import numpy as np
from typing import Tuple, List
from app.services.terrain_service import terrain_service, TerrainService


class LineOfSightService:
    """
    Line-of-Sight calculations with terrain
    """

    EARTH_RADIUS = 6371000  # meters
    K_FACTOR = 4 / 3  # Standard atmospheric refraction

    def __init__(self, terrain: TerrainService = None):
        self.terrain = terrain or terrain_service

    async def check_line_of_sight(
        self,
        tx_lat: float, tx_lon: float, tx_height: float,
        rx_lat: float, rx_lon: float, rx_height: float = 1.5,
        num_samples: int = 50
    ) -> dict:
        """
        Check line-of-sight between transmitter and receiver

        Args:
            tx_lat, tx_lon: Transmitter coordinates
            tx_height: Transmitter antenna height above ground (meters)
            rx_lat, rx_lon: Receiver coordinates
            rx_height: Receiver height above ground (meters), default 1.5m (person)
            num_samples: Number of points to sample along path

        Returns:
            {
                "has_los": bool,
                "clearance": float,  # minimum clearance in meters (negative = blocked)
                "blocked_at": float | None,  # distance where blocked (meters)
                "profile": [...]  # elevation profile with LOS line
            }
        """
        # Get elevation profile
        profile = await self.terrain.get_elevation_profile(
            tx_lat, tx_lon, rx_lat, rx_lon, num_samples
        )

        if not profile:
            return {"has_los": True, "clearance": 0, "blocked_at": None, "profile": []}

        # Get endpoint elevations
        tx_ground = profile[0]["elevation"]
        rx_ground = profile[-1]["elevation"]

        tx_total = tx_ground + tx_height
        rx_total = rx_ground + rx_height

        total_distance = profile[-1]["distance"]

        min_clearance = float('inf')
        blocked_at = None

        # Check each point along path
        for point in profile:
            d = point["distance"]
            terrain_elev = point["elevation"]

            if total_distance == 0:
                los_height = tx_total
            else:
                # Linear interpolation of LOS line
                los_height = tx_total + (rx_total - tx_total) * (d / total_distance)

            # Earth curvature correction (with atmospheric refraction)
            # Effective Earth radius = K * actual radius
            effective_radius = self.K_FACTOR * self.EARTH_RADIUS
            curvature = (d * (total_distance - d)) / (2 * effective_radius)

            # LOS height after curvature correction
            los_height_corrected = los_height - curvature

            # Clearance at this point
            clearance = los_height_corrected - terrain_elev

            # Add to profile for visualization
            point["los_height"] = los_height_corrected
            point["clearance"] = clearance

            if clearance < min_clearance:
                min_clearance = clearance
                if clearance <= 0:
                    blocked_at = d

        has_los = min_clearance > 0

        return {
            "has_los": has_los,
            "clearance": min_clearance,
            "blocked_at": blocked_at,
            "profile": profile
        }

    async def calculate_fresnel_clearance(
        self,
        tx_lat: float, tx_lon: float, tx_height: float,
        rx_lat: float, rx_lon: float, rx_height: float,
        frequency_mhz: float,
        num_samples: int = 50
    ) -> dict:
        """
        Calculate Fresnel zone clearance

        60% clearance of 1st Fresnel zone = good signal

        Returns:
            {
                "clearance_percent": float,  # worst-case clearance as % of required
                "has_adequate_clearance": bool,  # >= 60%
                "worst_point_distance": float,
                "fresnel_profile": [...]
            }
        """
        profile = await self.terrain.get_elevation_profile(
            tx_lat, tx_lon, rx_lat, rx_lon, num_samples
        )

        if not profile:
            return {
                "clearance_percent": 100.0,
                "has_adequate_clearance": True,
                "worst_point_distance": 0,
                "fresnel_profile": []
            }

        tx_ground = profile[0]["elevation"]
        rx_ground = profile[-1]["elevation"]

        tx_total = tx_ground + tx_height
        rx_total = rx_ground + rx_height

        total_distance = profile[-1]["distance"]

        if total_distance <= 0:
            return {
                "clearance_percent": 100.0,
                "has_adequate_clearance": True,
                "worst_point_distance": 0,
                "fresnel_profile": profile
            }

        # Wavelength (lambda = c / f)
        wavelength = 300.0 / frequency_mhz  # meters

        worst_clearance_pct = 100.0
        worst_distance = 0.0

        for point in profile:
            d = point["distance"]
            terrain_elev = point["elevation"]

            if d <= 0 or d >= total_distance:
                continue  # Skip endpoints

            # LOS height at this point
            los_height = tx_total + (rx_total - tx_total) * (d / total_distance)

            # 1st Fresnel zone radius at this point
            d1 = d
            d2 = total_distance - d
            fresnel_radius = float(np.sqrt((wavelength * d1 * d2) / total_distance))

            # Required clearance (60% of 1st Fresnel zone)
            required_clearance = 0.6 * fresnel_radius

            # Actual clearance
            actual_clearance = los_height - terrain_elev

            # Clearance as percentage of required
            if required_clearance > 0:
                clearance_pct = (actual_clearance / required_clearance) * 100
            else:
                clearance_pct = 100.0

            # Add to profile
            point["fresnel_radius"] = fresnel_radius
            point["required_clearance"] = required_clearance
            point["clearance_percent"] = clearance_pct

            if clearance_pct < worst_clearance_pct:
                worst_clearance_pct = clearance_pct
                worst_distance = d

        return {
            "clearance_percent": float(worst_clearance_pct),
            "has_adequate_clearance": worst_clearance_pct >= 60.0,
            "worst_point_distance": float(worst_distance),
            "fresnel_profile": profile
        }


# Singleton instance
los_service = LineOfSightService()