@mytec: iter1.3 ready for test

backend/app/services/coverage_service.py

@@ -0,0 +1,331 @@
+import numpy as np
+import asyncio
+from typing import List, Optional, Tuple
+from pydantic import BaseModel
+from app.services.terrain_service import terrain_service, TerrainService
+from app.services.los_service import los_service
+from app.services.buildings_service import buildings_service, Building
+
+
+class CoveragePoint(BaseModel):
+    lat: float
+    lon: float
+    rsrp: float  # dBm
+    distance: float  # meters from site
+    has_los: bool
+    terrain_loss: float  # dB
+    building_loss: float  # dB
+
+
+class CoverageSettings(BaseModel):
+    radius: float = 10000  # meters
+    resolution: float = 200  # meters
+    min_signal: float = -120  # dBm threshold
+    use_terrain: bool = True
+    use_buildings: bool = True
+
+
+class SiteParams(BaseModel):
+    lat: float
+    lon: float
+    height: float = 30  # antenna height meters
+    power: float = 43  # dBm (20W)
+    gain: float = 15  # dBi
+    frequency: float = 1800  # MHz
+    azimuth: Optional[float] = None  # degrees, None = omni
+    beamwidth: Optional[float] = 65  # degrees
+
+
+class CoverageService:
+    """
+    RF Coverage calculation with terrain and buildings
+    """
+
+    EARTH_RADIUS = 6371000
+
+    def __init__(self):
+        self.terrain = terrain_service
+        self.buildings = buildings_service
+        self.los = los_service
+
+    async def calculate_coverage(
+        self,
+        site: SiteParams,
+        settings: CoverageSettings
+    ) -> List[CoveragePoint]:
+        """
+        Calculate coverage grid for a single site
+
+        Returns list of CoveragePoint with RSRP values
+        """
+        points = []
+
+        # Generate grid
+        grid = self._generate_grid(
+            site.lat, site.lon,
+            settings.radius,
+            settings.resolution
+        )
+
+        # Fetch buildings for coverage area (if enabled)
+        buildings = []
+        if settings.use_buildings:
+            # Calculate bbox with margin
+            lat_delta = settings.radius / 111000  # ~111km per degree
+            lon_delta = settings.radius / (111000 * np.cos(np.radians(site.lat)))
+
+            buildings = await self.buildings.fetch_buildings(
+                site.lat - lat_delta, site.lon - lon_delta,
+                site.lat + lat_delta, site.lon + lon_delta
+            )
+
+        # Calculate coverage for each point
+        for lat, lon in grid:
+            point = await self._calculate_point(
+                site, lat, lon,
+                settings, buildings
+            )
+
+            if point.rsrp >= settings.min_signal:
+                points.append(point)
+
+        return points
+
+    async def calculate_multi_site_coverage(
+        self,
+        sites: List[SiteParams],
+        settings: CoverageSettings
+    ) -> List[CoveragePoint]:
+        """
+        Calculate combined coverage from multiple sites
+        Best server (strongest signal) wins at each point
+        """
+        if not sites:
+            return []
+
+        # Get all individual coverages
+        all_coverages = await asyncio.gather(*[
+            self.calculate_coverage(site, settings)
+            for site in sites
+        ])
+
+        # Combine by best signal
+        point_map: dict[Tuple[float, float], CoveragePoint] = {}
+
+        for coverage in all_coverages:
+            for point in coverage:
+                key = (round(point.lat, 6), round(point.lon, 6))
+
+                if key not in point_map or point.rsrp > point_map[key].rsrp:
+                    point_map[key] = point
+
+        return list(point_map.values())
+
+    def _generate_grid(
+        self,
+        center_lat: float, center_lon: float,
+        radius: float, resolution: float
+    ) -> List[Tuple[float, float]]:
+        """Generate coverage grid points"""
+        points = []
+
+        # Convert resolution to degrees
+        lat_step = resolution / 111000
+        lon_step = resolution / (111000 * np.cos(np.radians(center_lat)))
+
+        # Calculate grid bounds
+        lat_delta = radius / 111000
+        lon_delta = radius / (111000 * np.cos(np.radians(center_lat)))
+
+        lat = center_lat - lat_delta
+        while lat <= center_lat + lat_delta:
+            lon = center_lon - lon_delta
+            while lon <= center_lon + lon_delta:
+                # Check if within radius (circular, not square)
+                dist = TerrainService.haversine_distance(center_lat, center_lon, lat, lon)
+                if dist <= radius:
+                    points.append((lat, lon))
+                lon += lon_step
+            lat += lat_step
+
+        return points
+
+    async def _calculate_point(
+        self,
+        site: SiteParams,
+        lat: float, lon: float,
+        settings: CoverageSettings,
+        buildings: List[Building]
+    ) -> CoveragePoint:
+        """Calculate RSRP at a single point"""
+
+        # Distance
+        distance = TerrainService.haversine_distance(site.lat, site.lon, lat, lon)
+
+        if distance < 1:
+            distance = 1  # Avoid division by zero
+
+        # Base path loss (Okumura-Hata for urban)
+        path_loss = self._okumura_hata(
+            distance, site.frequency, site.height, 1.5  # 1.5m receiver height
+        )
+
+        # Antenna pattern loss (if directional)
+        antenna_loss = 0.0
+        if site.azimuth is not None and site.beamwidth:
+            antenna_loss = self._antenna_pattern_loss(
+                site.lat, site.lon, lat, lon,
+                site.azimuth, site.beamwidth
+            )
+
+        # Terrain loss (LoS check)
+        terrain_loss = 0.0
+        has_los = True
+
+        if settings.use_terrain:
+            los_result = await self.los.check_line_of_sight(
+                site.lat, site.lon, site.height,
+                lat, lon, 1.5  # receiver at 1.5m
+            )
+            has_los = los_result["has_los"]
+
+            if not has_los:
+                # Add diffraction loss based on clearance
+                clearance = los_result["clearance"]
+                terrain_loss = self._diffraction_loss(clearance, site.frequency)
+
+        # Building loss
+        building_loss = 0.0
+
+        if settings.use_buildings and buildings:
+            for building in buildings:
+                intersection = self.buildings.line_intersects_building(
+                    site.lat, site.lon, site.height + await self.terrain.get_elevation(site.lat, site.lon),
+                    lat, lon, 1.5 + await self.terrain.get_elevation(lat, lon),
+                    building
+                )
+                if intersection is not None:
+                    # Building penetration loss (~20dB for concrete)
+                    building_loss += 20.0
+                    has_los = False
+                    break  # One building is enough
+
+        # Calculate RSRP
+        # RSRP = Tx Power + Tx Gain - Path Loss - Antenna Loss - Terrain Loss - Building Loss
+        rsrp = site.power + site.gain - path_loss - antenna_loss - terrain_loss - building_loss
+
+        return CoveragePoint(
+            lat=lat,
+            lon=lon,
+            rsrp=rsrp,
+            distance=distance,
+            has_los=has_los,
+            terrain_loss=terrain_loss,
+            building_loss=building_loss
+        )
+
+    def _okumura_hata(
+        self,
+        distance: float,  # meters
+        frequency: float,  # MHz
+        tx_height: float,  # meters
+        rx_height: float   # meters
+    ) -> float:
+        """
+        Okumura-Hata path loss model (urban)
+
+        Returns path loss in dB
+        """
+        d_km = distance / 1000
+
+        if d_km < 0.1:
+            d_km = 0.1  # Minimum distance
+
+        # Mobile antenna height correction (urban)
+        a_hm = (1.1 * np.log10(frequency) - 0.7) * rx_height - (1.56 * np.log10(frequency) - 0.8)
+
+        # Path loss
+        L = (69.55 + 26.16 * np.log10(frequency) - 13.82 * np.log10(tx_height) - a_hm +
+             (44.9 - 6.55 * np.log10(tx_height)) * np.log10(d_km))
+
+        return L
+
+    def _antenna_pattern_loss(
+        self,
+        site_lat: float, site_lon: float,
+        point_lat: float, point_lon: float,
+        azimuth: float, beamwidth: float
+    ) -> float:
+        """Calculate antenna pattern attenuation"""
+        # Calculate bearing from site to point
+        bearing = self._calculate_bearing(site_lat, site_lon, point_lat, point_lon)
+
+        # Angle difference from main lobe
+        angle_diff = abs(bearing - azimuth)
+        if angle_diff > 180:
+            angle_diff = 360 - angle_diff
+
+        # Simple cosine pattern approximation
+        # 3dB beamwidth = angle where power drops to half
+        half_beamwidth = beamwidth / 2
+
+        if angle_diff <= half_beamwidth:
+            # Within main lobe - minimal loss
+            loss = 3 * (angle_diff / half_beamwidth) ** 2
+        else:
+            # Outside main lobe - significant loss
+            loss = 3 + 12 * ((angle_diff - half_beamwidth) / half_beamwidth) ** 2
+            loss = min(loss, 25)  # Cap at 25dB (back lobe)
+
+        return loss
+
+    def _calculate_bearing(
+        self,
+        lat1: float, lon1: float,
+        lat2: float, lon2: float
+    ) -> float:
+        """Calculate bearing from point 1 to point 2 (degrees)"""
+        lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
+
+        dlon = lon2 - lon1
+
+        x = np.sin(dlon) * np.cos(lat2)
+        y = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(dlon)
+
+        bearing = np.degrees(np.arctan2(x, y))
+
+        return (bearing + 360) % 360
+
+    def _diffraction_loss(self, clearance: float, frequency: float) -> float:
+        """
+        Knife-edge diffraction loss
+
+        Args:
+            clearance: Clearance in meters (negative = obstructed)
+            frequency: Frequency in MHz
+
+        Returns:
+            Additional loss in dB
+        """
+        if clearance >= 0:
+            return 0.0  # No obstruction
+
+        # Fresnel parameter approximation
+        # v ~ clearance * sqrt(2 / (lambda * d))
+        # Simplified: use clearance directly
+
+        v = abs(clearance) / 10  # Normalize
+
+        # Knife-edge loss approximation
+        if v <= 0:
+            loss = 0
+        elif v < 2.4:
+            loss = 6.02 + 9.11 * v - 1.27 * v**2
+        else:
+            loss = 13.0 + 20 * np.log10(v)
+
+        return min(loss, 40)  # Cap at 40dB
+
+
+# Singleton
+coverage_service = CoverageService()