@mytec: 1.4iter ready for testing
This commit is contained in:
394
backend/app/services/dominant_path_service.py
Normal file
394
backend/app/services/dominant_path_service.py
Normal file
@@ -0,0 +1,394 @@
|
||||
import numpy as np
|
||||
from typing import List, Tuple, Optional
|
||||
from dataclasses import dataclass
|
||||
from app.services.terrain_service import terrain_service
|
||||
from app.services.buildings_service import buildings_service, Building
|
||||
from app.services.materials_service import materials_service, BuildingMaterial
|
||||
|
||||
|
||||
@dataclass
|
||||
class RayPath:
|
||||
"""Single ray path from TX to RX"""
|
||||
path_type: str # "direct", "reflected", "diffracted", "street"
|
||||
total_distance: float # meters
|
||||
path_loss: float # dB
|
||||
reflection_points: List[Tuple[float, float]] # [(lat, lon), ...]
|
||||
materials_crossed: List[BuildingMaterial]
|
||||
is_valid: bool # Does this path exist?
|
||||
|
||||
|
||||
class DominantPathService:
|
||||
"""
|
||||
Find dominant propagation paths (2-3 strongest)
|
||||
|
||||
Path types:
|
||||
1. Direct (LoS if available)
|
||||
2. Single reflection off building
|
||||
3. Over-roof diffraction
|
||||
4. Around-corner diffraction
|
||||
"""
|
||||
|
||||
MAX_REFLECTIONS = 2
|
||||
MAX_PATHS = 3
|
||||
|
||||
async def find_dominant_paths(
|
||||
self,
|
||||
tx_lat: float, tx_lon: float, tx_height: float,
|
||||
rx_lat: float, rx_lon: float, rx_height: float,
|
||||
frequency_mhz: float,
|
||||
buildings: List[Building]
|
||||
) -> List[RayPath]:
|
||||
"""Find the dominant propagation paths"""
|
||||
|
||||
paths = []
|
||||
|
||||
# 1. Try direct path
|
||||
direct = await self._check_direct_path(
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz, buildings
|
||||
)
|
||||
if direct:
|
||||
paths.append(direct)
|
||||
|
||||
# 2. Try single-bounce reflections
|
||||
reflections = await self._find_reflection_paths(
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz, buildings
|
||||
)
|
||||
paths.extend(reflections[:2]) # Max 2 reflection paths
|
||||
|
||||
# 3. Try over-roof diffraction (if direct blocked)
|
||||
if not direct or not direct.is_valid:
|
||||
diffracted = await self._find_diffraction_path(
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz, buildings
|
||||
)
|
||||
if diffracted:
|
||||
paths.append(diffracted)
|
||||
|
||||
# Sort by path loss (best first) and return top N
|
||||
paths.sort(key=lambda p: p.path_loss)
|
||||
return paths[:self.MAX_PATHS]
|
||||
|
||||
async def _check_direct_path(
|
||||
self,
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz,
|
||||
buildings: List[Building]
|
||||
) -> Optional[RayPath]:
|
||||
"""Check if direct LoS path exists"""
|
||||
from app.services.los_service import los_service
|
||||
|
||||
# Check terrain LoS
|
||||
los_result = await los_service.check_line_of_sight(
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height
|
||||
)
|
||||
|
||||
if not los_result["has_los"]:
|
||||
distance = terrain_service.haversine_distance(tx_lat, tx_lon, rx_lat, rx_lon)
|
||||
return RayPath(
|
||||
path_type="direct",
|
||||
total_distance=distance,
|
||||
path_loss=float('inf'),
|
||||
reflection_points=[],
|
||||
materials_crossed=[],
|
||||
is_valid=False
|
||||
)
|
||||
|
||||
# Check building intersections
|
||||
materials_crossed = []
|
||||
for building in buildings:
|
||||
intersection = self._line_intersects_building_3d(
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
building
|
||||
)
|
||||
if intersection:
|
||||
material = materials_service.detect_material(building.tags)
|
||||
materials_crossed.append(material)
|
||||
|
||||
# Calculate path loss
|
||||
distance = terrain_service.haversine_distance(tx_lat, tx_lon, rx_lat, rx_lon)
|
||||
path_loss = self._calculate_path_loss(distance, frequency_mhz, tx_height, rx_height)
|
||||
|
||||
# Add material penetration losses
|
||||
for material in materials_crossed:
|
||||
path_loss += materials_service.get_penetration_loss(material, frequency_mhz)
|
||||
|
||||
return RayPath(
|
||||
path_type="direct",
|
||||
total_distance=distance,
|
||||
path_loss=path_loss,
|
||||
reflection_points=[],
|
||||
materials_crossed=materials_crossed,
|
||||
is_valid=len(materials_crossed) < 3 # Too many walls = not viable
|
||||
)
|
||||
|
||||
async def _find_reflection_paths(
|
||||
self,
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz,
|
||||
buildings: List[Building]
|
||||
) -> List[RayPath]:
|
||||
"""Find viable single-bounce reflection paths"""
|
||||
|
||||
reflection_paths = []
|
||||
|
||||
for building in buildings:
|
||||
# Find potential reflection points on building walls
|
||||
reflection_point = self._find_reflection_point(
|
||||
tx_lat, tx_lon, rx_lat, rx_lon, building
|
||||
)
|
||||
|
||||
if not reflection_point:
|
||||
continue
|
||||
|
||||
ref_lat, ref_lon = reflection_point
|
||||
|
||||
# Check if both segments are clear
|
||||
# TX -> Reflection point
|
||||
dist1 = terrain_service.haversine_distance(tx_lat, tx_lon, ref_lat, ref_lon)
|
||||
# Reflection point -> RX
|
||||
dist2 = terrain_service.haversine_distance(ref_lat, ref_lon, rx_lat, rx_lon)
|
||||
|
||||
total_distance = dist1 + dist2
|
||||
|
||||
# Don't consider if much longer than direct path
|
||||
direct_distance = terrain_service.haversine_distance(tx_lat, tx_lon, rx_lat, rx_lon)
|
||||
if total_distance > direct_distance * 2:
|
||||
continue
|
||||
|
||||
# Calculate path loss
|
||||
path_loss = self._calculate_path_loss(total_distance, frequency_mhz, tx_height, rx_height)
|
||||
|
||||
# Add reflection loss
|
||||
material = materials_service.detect_material(building.tags)
|
||||
path_loss += materials_service.get_reflection_loss(material)
|
||||
|
||||
reflection_paths.append(RayPath(
|
||||
path_type="reflected",
|
||||
total_distance=total_distance,
|
||||
path_loss=path_loss,
|
||||
reflection_points=[(ref_lat, ref_lon)],
|
||||
materials_crossed=[],
|
||||
is_valid=True
|
||||
))
|
||||
|
||||
return reflection_paths
|
||||
|
||||
async def _find_diffraction_path(
|
||||
self,
|
||||
tx_lat, tx_lon, tx_height,
|
||||
rx_lat, rx_lon, rx_height,
|
||||
frequency_mhz,
|
||||
buildings: List[Building]
|
||||
) -> Optional[RayPath]:
|
||||
"""Find over-roof diffraction path"""
|
||||
|
||||
# Find highest obstacle between TX and RX
|
||||
max_height = 0
|
||||
obstacle_lat, obstacle_lon = None, None
|
||||
|
||||
# Sample points along direct path
|
||||
num_samples = 20
|
||||
for i in range(1, num_samples - 1):
|
||||
t = i / num_samples
|
||||
lat = tx_lat + t * (rx_lat - tx_lat)
|
||||
lon = tx_lon + t * (rx_lon - tx_lon)
|
||||
|
||||
# Check terrain
|
||||
terrain_elev = await terrain_service.get_elevation(lat, lon)
|
||||
if terrain_elev > max_height:
|
||||
max_height = terrain_elev
|
||||
obstacle_lat, obstacle_lon = lat, lon
|
||||
|
||||
# Check buildings at this point
|
||||
for building in buildings:
|
||||
if buildings_service.point_in_building(lat, lon, building):
|
||||
if building.height > max_height:
|
||||
max_height = building.height
|
||||
obstacle_lat, obstacle_lon = lat, lon
|
||||
|
||||
if not obstacle_lat:
|
||||
return None
|
||||
|
||||
# Calculate diffraction loss (simplified knife-edge)
|
||||
distance = terrain_service.haversine_distance(tx_lat, tx_lon, rx_lat, rx_lon)
|
||||
|
||||
# Fresnel parameter
|
||||
tx_elev = await terrain_service.get_elevation(tx_lat, tx_lon)
|
||||
rx_elev = await terrain_service.get_elevation(rx_lat, rx_lon)
|
||||
|
||||
tx_total = tx_elev + tx_height
|
||||
rx_total = rx_elev + rx_height
|
||||
|
||||
# Height of LoS at obstacle point
|
||||
d1 = terrain_service.haversine_distance(tx_lat, tx_lon, obstacle_lat, obstacle_lon)
|
||||
los_height = tx_total + (rx_total - tx_total) * (d1 / distance) if distance > 0 else tx_total
|
||||
|
||||
clearance = los_height - max_height
|
||||
|
||||
# Knife-edge diffraction loss
|
||||
diffraction_loss = self._knife_edge_loss(clearance, frequency_mhz, distance, d1)
|
||||
|
||||
path_loss = self._calculate_path_loss(distance, frequency_mhz, tx_height, rx_height)
|
||||
path_loss += diffraction_loss
|
||||
|
||||
return RayPath(
|
||||
path_type="diffracted",
|
||||
total_distance=distance,
|
||||
path_loss=path_loss,
|
||||
reflection_points=[(obstacle_lat, obstacle_lon)],
|
||||
materials_crossed=[],
|
||||
is_valid=True
|
||||
)
|
||||
|
||||
def _find_reflection_point(
|
||||
self,
|
||||
tx_lat: float, tx_lon: float,
|
||||
rx_lat: float, rx_lon: float,
|
||||
building: Building
|
||||
) -> Optional[Tuple[float, float]]:
|
||||
"""Find specular reflection point on building wall"""
|
||||
|
||||
# Simplified: find closest wall segment and calculate reflection
|
||||
geometry = building.geometry
|
||||
|
||||
best_point = None
|
||||
best_score = float('inf')
|
||||
|
||||
for i in range(len(geometry) - 1):
|
||||
wall_start = geometry[i]
|
||||
wall_end = geometry[i + 1]
|
||||
|
||||
# Find reflection point on this wall segment
|
||||
ref_point = self._specular_reflection(
|
||||
tx_lon, tx_lat, rx_lon, rx_lat,
|
||||
wall_start[0], wall_start[1],
|
||||
wall_end[0], wall_end[1]
|
||||
)
|
||||
|
||||
if ref_point:
|
||||
# Score by total path length
|
||||
d1 = np.sqrt((ref_point[0] - tx_lon)**2 + (ref_point[1] - tx_lat)**2)
|
||||
d2 = np.sqrt((ref_point[0] - rx_lon)**2 + (ref_point[1] - rx_lat)**2)
|
||||
score = d1 + d2
|
||||
|
||||
if score < best_score:
|
||||
best_score = score
|
||||
best_point = (ref_point[1], ref_point[0]) # Return as (lat, lon)
|
||||
|
||||
return best_point
|
||||
|
||||
def _specular_reflection(
|
||||
self,
|
||||
tx_x, tx_y, rx_x, rx_y,
|
||||
wall_x1, wall_y1, wall_x2, wall_y2
|
||||
) -> Optional[Tuple[float, float]]:
|
||||
"""Calculate specular reflection point on wall segment"""
|
||||
|
||||
# Wall vector
|
||||
wall_dx = wall_x2 - wall_x1
|
||||
wall_dy = wall_y2 - wall_y1
|
||||
wall_len = np.sqrt(wall_dx**2 + wall_dy**2)
|
||||
|
||||
if wall_len < 1e-10:
|
||||
return None
|
||||
|
||||
# Wall normal
|
||||
normal_x = -wall_dy / wall_len
|
||||
normal_y = wall_dx / wall_len
|
||||
|
||||
# Mirror TX across wall
|
||||
# Project TX onto wall
|
||||
tx_rel_x = tx_x - wall_x1
|
||||
tx_rel_y = tx_y - wall_y1
|
||||
|
||||
dot = tx_rel_x * normal_x + tx_rel_y * normal_y
|
||||
|
||||
mirror_x = tx_x - 2 * dot * normal_x
|
||||
mirror_y = tx_y - 2 * dot * normal_y
|
||||
|
||||
# Find intersection of (mirror -> RX) with wall
|
||||
# Parametric line: mirror + t * (rx - mirror)
|
||||
dx = rx_x - mirror_x
|
||||
dy = rx_y - mirror_y
|
||||
|
||||
# Wall parametric: wall1 + s * (wall2 - wall1)
|
||||
denom = dx * wall_dy - dy * wall_dx
|
||||
|
||||
if abs(denom) < 1e-10:
|
||||
return None # Parallel
|
||||
|
||||
t = ((wall_x1 - mirror_x) * wall_dy - (wall_y1 - mirror_y) * wall_dx) / denom
|
||||
s = ((wall_x1 - mirror_x) * dy - (wall_y1 - mirror_y) * dx) / (-denom)
|
||||
|
||||
# Check if intersection is on wall segment and between mirror and RX
|
||||
if 0 <= s <= 1 and 0 <= t <= 1:
|
||||
ref_x = mirror_x + t * dx
|
||||
ref_y = mirror_y + t * dy
|
||||
return (ref_x, ref_y)
|
||||
|
||||
return None
|
||||
|
||||
def _line_intersects_building_3d(
|
||||
self,
|
||||
lat1, lon1, height1,
|
||||
lat2, lon2, height2,
|
||||
building: Building
|
||||
) -> bool:
|
||||
"""Check if 3D line intersects building volume"""
|
||||
# Sample along line
|
||||
for t in np.linspace(0, 1, 20):
|
||||
lat = lat1 + t * (lat2 - lat1)
|
||||
lon = lon1 + t * (lon2 - lon1)
|
||||
height = height1 + t * (height2 - height1)
|
||||
|
||||
if buildings_service.point_in_building(lat, lon, building):
|
||||
if height < building.height:
|
||||
return True
|
||||
return False
|
||||
|
||||
def _calculate_path_loss(self, distance, frequency_mhz, tx_height, rx_height) -> float:
|
||||
"""Okumura-Hata path loss"""
|
||||
d_km = max(distance / 1000, 0.1)
|
||||
|
||||
a_hm = (1.1 * np.log10(frequency_mhz) - 0.7) * rx_height - (1.56 * np.log10(frequency_mhz) - 0.8)
|
||||
|
||||
L = (69.55 + 26.16 * np.log10(frequency_mhz) - 13.82 * np.log10(tx_height) - a_hm +
|
||||
(44.9 - 6.55 * np.log10(tx_height)) * np.log10(d_km))
|
||||
|
||||
return L
|
||||
|
||||
def _knife_edge_loss(self, clearance, frequency_mhz, total_distance, d1) -> float:
|
||||
"""Knife-edge diffraction loss"""
|
||||
if clearance >= 0:
|
||||
return 0.0
|
||||
|
||||
wavelength = 300 / frequency_mhz
|
||||
d2 = total_distance - d1
|
||||
|
||||
if d1 <= 0 or d2 <= 0 or wavelength <= 0:
|
||||
return 0.0
|
||||
|
||||
# Fresnel parameter v
|
||||
v = abs(clearance) * np.sqrt(2 * (d1 + d2) / (wavelength * d1 * d2))
|
||||
|
||||
# Lee's approximation
|
||||
if v <= -0.78:
|
||||
return 0
|
||||
elif v < 0:
|
||||
return 6.02 + 9.11 * v - 1.27 * v**2
|
||||
elif v < 2.4:
|
||||
return 6.02 + 9.11 * v + 1.27 * v**2
|
||||
else:
|
||||
return 13 + 20 * np.log10(v)
|
||||
|
||||
|
||||
dominant_path_service = DominantPathService()
|
||||
Reference in New Issue
Block a user