import math
import numpy as np
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.models import BuildingFeature, SimulationRequest
from openworld_radio_twin.simulation.radio_map import RadioMapData
BOLTZMANN_J_PER_K = 1.380649e-23
PREVIEW_BANDWIDTH_HZ = 1_000_000.0
PREVIEW_TEMPERATURE_K = 293.0
def _point_in_polygon(x: float, y: float, polygon: list[tuple[float, float]]) -> bool:
inside = False
j = len(polygon) - 1
for i, (xi, yi) in enumerate(polygon):
xj, yj = polygon[j]
if (yi > y) != (yj > y):
crossing_x = (xj - xi) * (y - yi) / ((yj - yi) or 1e-12) + xi
if x < crossing_x:
inside = not inside
j = i
return inside
def _segment_intersects(
a: tuple[float, float],
b: tuple[float, float],
c: tuple[float, float],
d: tuple[float, float],
) -> bool:
def orientation(
p: tuple[float, float], q: tuple[float, float], r: tuple[float, float]
) -> float:
return (q[1] - p[1]) * (r[0] - q[0]) - (q[0] - p[0]) * (r[1] - q[1])
return (
orientation(a, b, c) * orientation(a, b, d) < 0
and orientation(c, d, a) * orientation(c, d, b) < 0
)
def _blocked_by_buildings(
tx_x: float,
tx_y: float,
x: float,
y: float,
tx_height: float,
rx_height: float,
polygons: list[tuple[float, float, list[tuple[float, float]]]],
) -> int:
blocked = 0
distance = math.hypot(x - tx_x, y - tx_y)
if distance < 1:
return 0
for minimum_height, maximum_height, polygon in polygons:
if _point_in_polygon(x, y, polygon):
continue
for index, start in enumerate(polygon[:-1]):
end = polygon[index + 1]
if not _segment_intersects((tx_x, tx_y), (x, y), start, end):
continue
wall_distance = min(distance, math.hypot(start[0] - tx_x, start[1] - tx_y))
ray_height = tx_height + (rx_height - tx_height) * (wall_distance / distance)
if minimum_height < ray_height < maximum_height:
blocked += 1
break
if blocked >= 3:
break
return blocked
[docs]
def simulate_preview(request: SimulationRequest, buildings: list[BuildingFeature]) -> RadioMapData:
"""Return an analytical preview in the same numerical contract as Sionna RT."""
center = request.scene_origin
diameter = request.radius_m * 2
width = max(1, math.ceil(diameter / request.resolution_m))
height = width
step_x = diameter / width
step_y = diameter / height
local_frame = LocalFrame.at(center.longitude, center.latitude)
polygons = [
(
building.min_height_m,
building.roof_height_agl_m,
[local_frame.to_local(lon, lat)[:2] for lon, lat in building.coordinates],
)
for building in buildings
]
transmitters = [
(
transmitter,
local_frame.to_local(
transmitter.position.longitude,
transmitter.position.latitude,
transmitter.position.altitude_m,
),
)
for transmitter in request.transmitters
]
path_gain = np.empty((len(transmitters), height, width), dtype=np.float32)
rss = np.empty_like(path_gain)
centers = np.empty((height, width, 3), dtype=np.float64)
for row in range(height):
north = -request.radius_m + (row + 0.5) * step_y
for column in range(width):
east = -request.radius_m + (column + 0.5) * step_x
centers[row, column] = (east, north, request.receiver_height_m)
for tx_index, (tx, (tx_east, tx_north, tx_up)) in enumerate(transmitters):
delta_east = east - tx_east
delta_north = north - tx_north
distance_2d = max(1.0, math.hypot(delta_east, delta_north))
distance_3d = math.hypot(distance_2d, tx_up - request.receiver_height_m)
fspl = (
32.44
+ 20 * math.log10(tx.frequency_ghz * 1000)
+ 20 * math.log10(distance_3d / 1000)
)
obstruction_count = _blocked_by_buildings(
tx_east,
tx_north,
east,
north,
tx_up,
request.receiver_height_m,
polygons,
)
directional_penalty = 0.0
if tx.antenna_pattern == "sector":
bearing = (math.degrees(math.atan2(delta_east, delta_north)) + 360) % 360
delta = abs((bearing - tx.azimuth_deg + 180) % 360 - 180)
directional_penalty = min(28.0, (delta / 65.0) ** 2 * 12.0)
loss_db = fspl + obstruction_count * 13.0 + directional_penalty
path_gain[tx_index, row, column] = 10.0 ** (-loss_db / 10.0)
transmit_power_w = 10.0 ** ((tx.power_dbm - 30.0) / 10.0)
rss[tx_index, row, column] = transmit_power_w * path_gain[
tx_index, row, column
]
noise_w = BOLTZMANN_J_PER_K * PREVIEW_TEMPERATURE_K * PREVIEW_BANDWIDTH_HZ
total_rss = np.sum(rss.astype(np.float64), axis=0, keepdims=True)
interference = np.maximum(0.0, total_rss - rss)
sinr = (rss / (interference + noise_w)).astype(np.float32)
association_values = {
"path_gain": path_gain,
"rss": rss,
"sinr": sinr,
}[request.association_metric]
association = np.argmax(association_values, axis=0).astype(np.int32)
association[np.max(association_values, axis=0) <= 0] = -1
return RadioMapData(path_gain, rss, sinr, association, centers)
