Rename the simulator to Overwatch and organize the repository for presentation.
The operator interface now uses that name. The radar design, parts list, and simulation sit in separate directories, and the earlier simulator and project license are removed. Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
10
.github/workflows/ci-tests.yml
vendored
10
.github/workflows/ci-tests.yml
vendored
@@ -46,8 +46,8 @@ jobs:
|
||||
- name: Unit tests
|
||||
run: >
|
||||
uv run pytest
|
||||
9_Firmware/9_3_GUI/test_GUI_V65_Tk.py
|
||||
9_Firmware/9_3_GUI/test_v7.py
|
||||
firmware/9_3_GUI/test_GUI_V65_Tk.py
|
||||
firmware/9_3_GUI/test_v7.py
|
||||
-v --tb=short
|
||||
|
||||
# ===========================================================================
|
||||
@@ -66,7 +66,7 @@ jobs:
|
||||
|
||||
- name: Build and run MCU tests
|
||||
run: make test
|
||||
working-directory: 9_Firmware/9_1_Microcontroller/tests
|
||||
working-directory: firmware/9_1_Microcontroller/tests
|
||||
|
||||
# ===========================================================================
|
||||
# FPGA RTL Regression (25 testbenches + lint)
|
||||
@@ -83,7 +83,7 @@ jobs:
|
||||
|
||||
- name: Run full FPGA regression
|
||||
run: bash run_regression.sh
|
||||
working-directory: 9_Firmware/9_2_FPGA
|
||||
working-directory: firmware/9_2_FPGA
|
||||
|
||||
# ===========================================================================
|
||||
# Cross-Layer Contract Tests (Python ↔ Verilog ↔ C)
|
||||
@@ -112,5 +112,5 @@ jobs:
|
||||
- name: Run cross-layer contract tests
|
||||
run: >
|
||||
uv run pytest
|
||||
9_Firmware/tests/cross_layer/test_cross_layer_contract.py
|
||||
firmware/tests/cross_layer/test_cross_layer_contract.py
|
||||
-v --tb=short
|
||||
|
||||
60
.gitignore
vendored
60
.gitignore
vendored
@@ -3,35 +3,35 @@
|
||||
*.vcd
|
||||
|
||||
# Debug / scratch RTL (not part of the design)
|
||||
9_Firmware/9_2_FPGA/debug_*.v
|
||||
9_Firmware/9_2_FPGA/tb/tb_fft_debug*.v
|
||||
9_Firmware/9_2_FPGA/tb/tb_fft_mini*.v
|
||||
9_Firmware/9_2_FPGA/tb/tb_bram_debug.v
|
||||
firmware/9_2_FPGA/debug_*.v
|
||||
firmware/9_2_FPGA/tb/tb_fft_debug*.v
|
||||
firmware/9_2_FPGA/tb/tb_fft_mini*.v
|
||||
firmware/9_2_FPGA/tb/tb_bram_debug.v
|
||||
|
||||
# Local simulation artifacts and CSV outputs
|
||||
9_Firmware/9_2_FPGA/cic_*.csv
|
||||
9_Firmware/9_2_FPGA/fir_*.csv
|
||||
9_Firmware/9_2_FPGA/nco_*.csv
|
||||
9_Firmware/9_2_FPGA/ddc_*.csv
|
||||
9_Firmware/9_2_FPGA/mf_pipeline_output.csv
|
||||
9_Firmware/9_2_FPGA/mf_chain_autocorr.csv
|
||||
9_Firmware/9_2_FPGA/rbd_mode00_ramp.csv
|
||||
9_Firmware/9_2_FPGA/rbd_mode01_peak.csv
|
||||
9_Firmware/9_2_FPGA/rbd_mode10_avg.csv
|
||||
9_Firmware/9_2_FPGA/rbd_mode10_ramp.csv
|
||||
9_Firmware/9_2_FPGA/rmc_autoscan.csv
|
||||
9_Firmware/9_2_FPGA/tb/mf_chain_autocorr.csv
|
||||
9_Firmware/9_2_FPGA/tb/rbd_mode00_ramp.csv
|
||||
9_Firmware/9_2_FPGA/tb/rbd_mode01_peak.csv
|
||||
9_Firmware/9_2_FPGA/tb/rbd_mode10_avg.csv
|
||||
9_Firmware/9_2_FPGA/tb/rbd_mode10_ramp.csv
|
||||
9_Firmware/9_2_FPGA/tb/rmc_autoscan.csv
|
||||
9_Firmware/9_2_FPGA/tb_usb_data_interface.csv
|
||||
firmware/9_2_FPGA/cic_*.csv
|
||||
firmware/9_2_FPGA/fir_*.csv
|
||||
firmware/9_2_FPGA/nco_*.csv
|
||||
firmware/9_2_FPGA/ddc_*.csv
|
||||
firmware/9_2_FPGA/mf_pipeline_output.csv
|
||||
firmware/9_2_FPGA/mf_chain_autocorr.csv
|
||||
firmware/9_2_FPGA/rbd_mode00_ramp.csv
|
||||
firmware/9_2_FPGA/rbd_mode01_peak.csv
|
||||
firmware/9_2_FPGA/rbd_mode10_avg.csv
|
||||
firmware/9_2_FPGA/rbd_mode10_ramp.csv
|
||||
firmware/9_2_FPGA/rmc_autoscan.csv
|
||||
firmware/9_2_FPGA/tb/mf_chain_autocorr.csv
|
||||
firmware/9_2_FPGA/tb/rbd_mode00_ramp.csv
|
||||
firmware/9_2_FPGA/tb/rbd_mode01_peak.csv
|
||||
firmware/9_2_FPGA/tb/rbd_mode10_avg.csv
|
||||
firmware/9_2_FPGA/tb/rbd_mode10_ramp.csv
|
||||
firmware/9_2_FPGA/tb/rmc_autoscan.csv
|
||||
firmware/9_2_FPGA/tb_usb_data_interface.csv
|
||||
|
||||
# Co-sim intermediate CSVs (regenerated by scripts)
|
||||
9_Firmware/9_2_FPGA/tb/cosim/rtl_doppler_*.csv
|
||||
9_Firmware/9_2_FPGA/tb/cosim/compare_doppler_*.csv
|
||||
9_Firmware/9_2_FPGA/tb/cosim/rtl_multiseg_*.csv
|
||||
firmware/9_2_FPGA/tb/cosim/rtl_doppler_*.csv
|
||||
firmware/9_2_FPGA/tb/cosim/compare_doppler_*.csv
|
||||
firmware/9_2_FPGA/tb/cosim/rtl_multiseg_*.csv
|
||||
|
||||
# macOS
|
||||
.DS_Store
|
||||
@@ -44,13 +44,13 @@ __pycache__/
|
||||
10_docs/
|
||||
|
||||
# Local simulation workspaces and generated outputs
|
||||
5_Simulations/generated/
|
||||
5_Simulations/aeris10_antenna_sim.py
|
||||
5_Simulations/aeris10_radar_sim.py
|
||||
hardware/rf-simulations/generated/
|
||||
hardware/rf-simulations/aeris10_antenna_sim.py
|
||||
hardware/rf-simulations/aeris10_radar_sim.py
|
||||
|
||||
# Local FPGA report dumps and scratch constraints
|
||||
9_Firmware/9_2_FPGA/reports/
|
||||
9_Firmware/9_2_FPGA/synth_only.xdc
|
||||
firmware/9_2_FPGA/reports/
|
||||
firmware/9_2_FPGA/synth_only.xdc
|
||||
|
||||
# Local timing closure report snapshots
|
||||
build*_reports/
|
||||
|
||||
@@ -20,11 +20,13 @@ Never run bare package installation commands without the `sfw` prefix.
|
||||
|
||||
| Path | Contents |
|
||||
|------|----------|
|
||||
| `4_Schematics and Boards Layout/` | KiCad schematics, Gerbers, BOM/CPL |
|
||||
| `9_Firmware/9_1_Microcontroller/` | STM32 MCU C/C++ firmware and unit tests |
|
||||
| `9_Firmware/9_2_FPGA/` | Verilog RTL, constraints, testbenches, build scripts |
|
||||
| `9_Firmware/9_3_GUI/` | Python radar dashboard (Tkinter/PyQt6) and CLI tools |
|
||||
| `9_Firmware/tests/cross_layer/` | Python-based system invariant/contract tests |
|
||||
| `simulation/` | Overwatch |
|
||||
| `schematics/` | Schematics, Gerbers, BOM/CPL |
|
||||
| `parts/` | Consolidated parts list and datasheets |
|
||||
| `firmware/9_1_Microcontroller/` | STM32 MCU C/C++ firmware and unit tests |
|
||||
| `firmware/9_2_FPGA/` | Verilog RTL, constraints, testbenches, build scripts |
|
||||
| `firmware/9_3_GUI/` | Python radar dashboard (Tkinter/PyQt6) and CLI tools |
|
||||
| `firmware/tests/cross_layer/` | Python-based system invariant/contract tests |
|
||||
| `docs/` | GitHub Pages documentation site |
|
||||
|
||||
## Code Standards & Tooling
|
||||
@@ -57,26 +59,26 @@ We use GitHub Actions for CI, which runs four main jobs on every PR. Run these l
|
||||
### 1. Python & Linting
|
||||
```bash
|
||||
uv run ruff check .
|
||||
cd 9_Firmware/9_3_GUI
|
||||
cd firmware/9_3_GUI
|
||||
uv run pytest test_GUI_V65_Tk.py test_v7.py -v
|
||||
```
|
||||
|
||||
### 2. FPGA Regression
|
||||
```bash
|
||||
cd 9_Firmware/9_2_FPGA
|
||||
cd firmware/9_2_FPGA
|
||||
bash run_regression.sh
|
||||
```
|
||||
This runs five phases (Lint, Changed Modules, Integration, Signal Processing, Infrastructure, and **P0 Adversarial Tests**). All must pass.
|
||||
|
||||
### 3. MCU Unit Tests
|
||||
```bash
|
||||
cd 9_Firmware/9_1_Microcontroller/tests
|
||||
cd firmware/9_1_Microcontroller/tests
|
||||
make clean && make
|
||||
```
|
||||
|
||||
### 4. Cross-Layer Contract Tests
|
||||
```bash
|
||||
uv run pytest 9_Firmware/tests/cross_layer/test_cross_layer_contract.py -v
|
||||
uv run pytest firmware/tests/cross_layer/test_cross_layer_contract.py -v
|
||||
```
|
||||
|
||||
## Before merging: CI checklist
|
||||
@@ -102,9 +104,9 @@ All PRs must pass CI:
|
||||
|
||||
- [ ] `uv run ruff check .` — no lint errors
|
||||
- [ ] `uv run pytest test_GUI_V65_Tk.py test_v7.py -v` — all pass
|
||||
- [ ] `cd 9_Firmware/9_2_FPGA && bash run_regression.sh` — all 5 phases pass
|
||||
- [ ] `cd 9_Firmware/9_1_Microcontroller/tests && make clean && make` — pass
|
||||
- [ ] `uv run pytest 9_Firmware/tests/cross_layer/test_cross_layer_contract.py` — pass
|
||||
- [ ] `cd firmware/9_2_FPGA && bash run_regression.sh` — all 5 phases pass
|
||||
- [ ] `cd firmware/9_1_Microcontroller/tests && make clean && make` — pass
|
||||
- [ ] `uv run pytest firmware/tests/cross_layer/test_cross_layer_contract.py` — pass
|
||||
- [ ] `git diff --check` — no whitespace issues
|
||||
- [ ] PR targets `develop` branch
|
||||
|
||||
|
||||
@@ -1,206 +0,0 @@
|
||||
# Generalized Projectile Defense Simulation
|
||||
|
||||
## Overview
|
||||
|
||||
This simulation models a realistic defense scenario where an attack projectile (drone, missile, or ballistic re-entry vehicle) approaches a city defended by a ring of radars and interceptor drones. The simulation uses proper physics, realistic flight profiles, and validated error models.
|
||||
|
||||
It is **generalized** — the attack profile is selected via **presets** for different projectile types. Each preset defines the speed, cruise altitude band, descent geometry, and terminal altitude. The default preset is the jet-powered **Shahed-238 (Geran-3)**. All attack parameters are also editable manually for custom scenarios.
|
||||
|
||||
## Projectile Presets
|
||||
|
||||
| Preset | Type | Speed | Cruise Alt | Descent Start | Terminal Alt | Notes |
|
||||
|--------|------|-------|-----------|---------------|--------------|-------|
|
||||
| `shahed-136` | Prop drone (Geran-2) | 50 m/s | 100–300 m | 5 km | 0 m | Slow, low-altitude |
|
||||
| `shahed-238` | Jet drone (Geran-3) | 150 m/s | 1000–2000 m | 8 km | 200 m | Default; fast, mid-altitude |
|
||||
| `cruise-missile` | Cruise missile | 300 m/s | 50–150 m | 5 km | 0 m | Fast, sea-skimming |
|
||||
| `ballistic-reentry` | Ballistic re-entry | 800 m/s | 60000–80000 m | 40 km | 0 m | Very fast, high-altitude entry |
|
||||
| `prop-drone` | Small UAS | 30 m/s | 100–300 m | 3 km | 0 m | Slow, low-altitude |
|
||||
|
||||
### Shahed-238 Specifications (Jet-Powered Variant)
|
||||
|
||||
Based on research of open-source intelligence and technical analysis:
|
||||
|
||||
| Parameter | Value | Notes |
|
||||
|-----------|-------|-------|
|
||||
| Speed (cruise) | ~150 m/s (540 km/h) | 3-4x faster than Shahed-136 |
|
||||
| Cruise altitude | 1000-2000 m | Typical operational band |
|
||||
| Terminal altitude | 200 m | Descent target altitude |
|
||||
| Descent start | ~8 km from target | Gradual descent begins |
|
||||
| Warhead | 30-50 kg | High explosive |
|
||||
| Guidance | GPS/INS + terminal seeker | Autonomous |
|
||||
| Range | ~1500-2000 km | One-way attack |
|
||||
|
||||
Sources:
|
||||
- [CSIS Missile Threat - Shahed-238](https://missilethreat.csis.org/missile/shahed-238/)
|
||||
- [BattlePolicy - Shahed-238](https://www.battlepolicy.com/shahed-238/)
|
||||
- [Wikipedia - HESA Shahed 136](https://en.wikipedia.org/wiki/HESA_Shahed_136)
|
||||
- [drone-warfare.com - Shahed variants](https://drone-warfare.com/research/shahed-136/)
|
||||
|
||||
## Flight Profile
|
||||
|
||||
Each preset follows a realistic two-phase flight profile:
|
||||
|
||||
1. **Cruise phase**: Maintains constant altitude (within the preset's cruise band, randomized per run) at cruise speed
|
||||
2. **Descent phase**: Begins `descent_start` distance from the target, gradually descends to the preset's terminal altitude
|
||||
|
||||
The descent is gradual (not a nosedive), consistent with reported loitering-munition attack patterns. The projectile uses GPS/INS navigation to follow a pre-planned waypoint course, with a terminal seeker for the final approach. For the ballistic re-entry preset, the high cruise altitude and long descent start model a high-altitude entry with a terminal pull-down.
|
||||
|
||||
## Defense System
|
||||
|
||||
### Radar Network
|
||||
- **8 AERIS-10 radars** positioned in a ring around the city (2 km outside city radius)
|
||||
- Each radar has **20 km range**
|
||||
- Single-measurement accuracy: **±20 m (1σ)**
|
||||
- Detection requires target to be within radar range
|
||||
- Multiple radars improve position accuracy: σ_combined = σ_single / √n
|
||||
|
||||
### Interceptor Drones
|
||||
- **24 drones** positioned in a ring around the city (8 km outside city radius)
|
||||
- Drone speed: **100 m/s** (slower than most presets' cruise speed, e.g. Shahed-238 at 150 m/s)
|
||||
- IR homing range: **5 km**
|
||||
- IR tracking accuracy: **±3 m (1σ)**
|
||||
- Kill radius: **20 m** (3D)
|
||||
- **3 closest drones** launch per threat (selected by fastest viable intercept time)
|
||||
|
||||
## Physics Models
|
||||
|
||||
### 1. Radar Detection
|
||||
Each radar measures:
|
||||
- Range (with Gaussian noise, σ = 20 m)
|
||||
- Azimuth (with noise proportional to 1/distance)
|
||||
- Elevation (with noise proportional to 1/distance)
|
||||
|
||||
The noise model reflects real radar behavior: angular accuracy improves with range, but absolute position error grows.
|
||||
|
||||
### 2. Position Estimation
|
||||
The 3D position is estimated by averaging the individual 3D point estimates from each radar (bearing + range). This is a least-squares estimate that properly constrains all three coordinates (x, y, z), including altitude.
|
||||
|
||||
**Key insight**: With coplanar radars (all at z=0), range measurements alone cannot determine altitude. The bearing (azimuth + elevation) measurements are required to constrain the z-coordinate.
|
||||
|
||||
### 3. Trajectory Estimation
|
||||
The velocity is estimated by linear regression of position vs. time:
|
||||
```
|
||||
pos(t) = pos0 + vel * t
|
||||
```
|
||||
The slope (vel) is the estimated velocity. R² measures fit quality. The standard error of the slope decreases as more data points accumulate: σ_vel = σ_pos / √(Sxx), where Sxx = Σ(t - t̄)².
|
||||
|
||||
**This is the key mechanism by which error decreases with more data points**: the best-fit line converges as more observations are added.
|
||||
|
||||
### 4. Intercept Calculation
|
||||
The drone intercept solves the classic pursuer-target problem:
|
||||
```
|
||||
|drone_pos + v_d * t - target_pos - v_s * t| = 0
|
||||
```
|
||||
where v_d = drone_speed * d (d is a unit vector).
|
||||
|
||||
This reduces to a quadratic in t:
|
||||
```
|
||||
(|v_s|² - v_d²) * t² + 2*(rel_pos · v_s) * t + |rel_pos|² = 0
|
||||
```
|
||||
|
||||
The **smallest positive root** gives the earliest intercept time. The intercept point is where the target will be at that time.
|
||||
|
||||
**Critical insight**: A slower pursuer (drone at 100 m/s vs a faster target such as Shahed-238 at 150 m/s) can only intercept if:
|
||||
- The target is moving TOWARD the pursuer (rel_pos · v_s < 0)
|
||||
- The geometry is favorable (drones positioned in the approach corridor)
|
||||
- The intercept happens before the target reaches the city
|
||||
|
||||
### 5. Drone Guidance
|
||||
Drones fly in 3D straight lines toward the computed intercept point. The intercept point is recomputed every 5 seconds as new trajectory estimates become available (proportional navigation). When within 5 km, the drone switches to IR homing for final precision.
|
||||
|
||||
### 6. Error Propagation
|
||||
- Radar position error: σ_pos = 20/√n m (decreases with data points)
|
||||
- 95% confidence interval: 1.96 * σ_pos
|
||||
- IR error: σ_IR = 3/√n_drones m (decreases with tracking drones)
|
||||
- Combined error: σ_combined = √(σ_radar² + σ_IR²)
|
||||
|
||||
## Validation Tests
|
||||
|
||||
The simulation is validated against known analytical solutions:
|
||||
|
||||
| Test | Description | Result |
|
||||
|------|-------------|--------|
|
||||
| Trajectory continuity | Target speed is constant, altitude profile correct | PASS |
|
||||
| 1D intercept | Analytical solution: t = 48 s for specific geometry | PASS (t = 48.0 s) |
|
||||
| No-solution case | Target moving away, no intercept possible | PASS (graceful fallback) |
|
||||
| 3D intercept | Intercept altitude matches target altitude | PASS |
|
||||
| Radar detection | In-range detected, out-of-range not | PASS |
|
||||
| Position estimation | Error < 30 m with 3 radars | PASS (10.2 m) |
|
||||
| Trajectory estimation | Velocity error < 10 m/s | PASS (0.00 m/s) |
|
||||
| Simulation runs | Full simulation completes without errors | PASS |
|
||||
| Interception rate | Rate is within [0, 100]% for default preset | PASS (~67%) |
|
||||
| Presets vary rates | Slow/low target >= fast/high target intercept rate | PASS |
|
||||
| Drone selection | Only viable drones launched | PASS |
|
||||
|
||||
**All 11 tests pass.**
|
||||
|
||||
## Monte Carlo Results (100 Trials per Preset)
|
||||
|
||||
| Preset | Interception Rate | Avg Detect | Avg Collision | Avg Impact |
|
||||
|--------|-------------------|-----------|---------------|-----------|
|
||||
| `shahed-136` (prop, low) | **100%** | 7.0 s | 88.5 s | — |
|
||||
| `shahed-238` (jet, mid) | **68%** | 27.2 s | 85.8 s | 170.0 s |
|
||||
| `cruise-missile` (fast, low) | **36%** | 64.7 s | 98.9 s | 137.6 s |
|
||||
| `ballistic-reentry` (very fast, high) | **0%** | 110.9 s | — | 114.6 s |
|
||||
| `prop-drone` (slow, low) | **100%** | 1.0 s | 44.2 s | — |
|
||||
|
||||
(Exact values are regenerated each run via `montecarlo_stats.json`; the table shows representative values.)
|
||||
|
||||
## Key Findings
|
||||
|
||||
1. **Interception rate scales with target speed and altitude.** Slow, low-altitude projectiles (Shahed-136, small UAS) are nearly always intercepted, while fast, high-altitude ones (ballistic re-entry) are effectively unstoppable by this drone defense. The jet-powered Shahed-238 sits in the middle at ~68%.
|
||||
|
||||
2. **Detection happens before the target is in the danger zone.** The target enters radar range (20 km) and is detected after accumulating 3+ observations. Faster targets get detected later in absolute terms but with less margin.
|
||||
|
||||
3. **The interception window is the gap between detection and impact.** When interception succeeds, it occurs well before the impact time. When it fails, the target reaches the city. This gap is the decision time available to the defense.
|
||||
|
||||
4. **The interception rate depends strongly on geometry.** Approaches from directions where the drones are well-positioned in the approach corridor have higher interception rates. Approaches from directions where the drones are on the opposite side of the city fail.
|
||||
|
||||
5. **Altitude matters.** Higher cruise altitudes give the drones more time to climb, but also increase the distance they must cover. The optimal altitude depends on the approach geometry.
|
||||
|
||||
## Why the Rates Are Realistic
|
||||
|
||||
The per-preset interception rates reflect the fundamental physics of a slower defender (100 m/s drones) vs. faster or equal-speed attackers:
|
||||
|
||||
- A slow target (30–50 m/s) is easily caught — the drones have ample time to close the gap, so interception is near 100%.
|
||||
- A mid-speed jet target (150 m/s) outruns part of the pursuit; only favorable geometries allow a timely intercept, giving ~68%.
|
||||
- A fast sea-skimming cruise missile (300 m/s) is hard to catch, especially at low altitude where climb time hurts; ~36%.
|
||||
- A very fast high-altitude ballistic re-entry (800 m/s) is far beyond the drones' capability; 0%.
|
||||
|
||||
This is consistent with real-world observations that drone-vs-drone interception is challenging against faster targets, and that the interception rate depends heavily on the target's speed, altitude, approach direction, and the positioning of the defenders.
|
||||
|
||||
## Files
|
||||
|
||||
| File | Description |
|
||||
|------|-------------|
|
||||
| `defense_sim.py` | Python simulation with proper physics and projectile presets |
|
||||
| `test_defense_sim.py` | Validation test suite (11 tests) |
|
||||
| `defense_sim.html` | Interactive HTML visualization with preset dropdown |
|
||||
| `montecarlo_stats.json` | Monte Carlo results (100 trials per preset) |
|
||||
|
||||
## Running the Simulation
|
||||
|
||||
```bash
|
||||
# Run the Python simulation (runs all presets)
|
||||
python3 defense_sim.py
|
||||
|
||||
# Run the validation tests
|
||||
python3 test_defense_sim.py
|
||||
|
||||
# Open the HTML viewer in a browser (serve over HTTP for the stats table)
|
||||
python3 -m http.server 8000
|
||||
# then open http://localhost:8000/defense_sim.html
|
||||
```
|
||||
|
||||
### Selecting a Projectile Preset
|
||||
|
||||
- **In the HTML:** use the **Projectile Preset** dropdown at the top of the Attack Parameters panel. Selecting a preset prefills all attack parameters (start distance, speed, cruise altitude band, descent start, terminal altitude) and highlights that preset's row in the stats table. Choose **Custom** to set all parameters manually.
|
||||
- **In Python:** pass `preset='shahed-238'` (or any preset name) to `run_simulation()` or `monte_carlo()`. Omit it to use the default (`shahed-238`).
|
||||
|
||||
## References
|
||||
|
||||
1. CSIS Missile Threat Project - Shahed-238: https://missilethreat.csis.org/missile/shahed-238/
|
||||
2. BattlePolicy - Shahed-238: https://www.battlepolicy.com/shahed-238/
|
||||
3. Wikipedia - HESA Shahed 136: https://en.wikipedia.org/wiki/HESA_Shahed_136
|
||||
4. drone-warfare.com - Shahed-136 and Geran variants: https://drone-warfare.com/research/shahed-136/
|
||||
5. AERIS-10 radar: 10.5 GHz X-band FMCW phased array, 20 km range
|
||||
154
Licence
154
Licence
@@ -1,154 +0,0 @@
|
||||
CERN Open Hardware Licence Version 2 - Permissive
|
||||
|
||||
Preamble
|
||||
|
||||
CERN has developed this licence to promote collaboration among
|
||||
hardware designers and to provide a legal tool for the dissemination
|
||||
of hardware designs. The CERN Open Hardware Licence Version 2 is
|
||||
available in two variants: Weakly Reciprocal (CERN-OHL-W) and
|
||||
Permissive (CERN-OHL-P). The CERN-OHL-P is for you if you want to
|
||||
maximise the ability of people to use your hardware design
|
||||
including for use in commercial applications with a minimum of
|
||||
obligations. It is a continuation of the CERN Open Hardware Licence
|
||||
version 1 (but with a new structure) and is focused more narrowly
|
||||
on the design itself and what you must do to give proper credit and
|
||||
to allow others to build and modify the design.
|
||||
|
||||
Please see the CERN-OHL version 2 family page (www.ohwr.org/cernohl)
|
||||
for more information.
|
||||
|
||||
Terms and conditions
|
||||
|
||||
1. Definitions
|
||||
|
||||
In this Licence, the following terms have the following meanings:
|
||||
|
||||
“Additive” means something you add to the Design that was not part
|
||||
of it when it was first published by the Licensor.
|
||||
|
||||
“Design” means the design documentation for the hardware that is
|
||||
licensed under this Licence, in any format or medium. It may
|
||||
include both the Hardware Description and the Modifications.
|
||||
|
||||
“Documentation” means all the materials that are provided by the
|
||||
Licensor or any Contributor to describe the Design, including
|
||||
schematics, netlists, PCB layout files, bills of materials,
|
||||
drawings, specifications, and the like, and also includes the
|
||||
Modifications.
|
||||
|
||||
“Hardware” means any physical device that embodies the Design or
|
||||
any part thereof, in any form.
|
||||
|
||||
“Hardware Description” means the files that specify the Hardware in
|
||||
a form that can be read by a computer and that, when used
|
||||
appropriately, can be used to fabricate Hardware, for example in
|
||||
whole or in part.
|
||||
|
||||
“Licence” means this licence (the CERN-OHL-P).
|
||||
|
||||
“Licensor” means the person or entity that first publishes the
|
||||
Design under the terms of this Licence.
|
||||
|
||||
“Modification” means any change to the Design, including changes to
|
||||
the Hardware Description.
|
||||
|
||||
“Product” means either an unchanged Design or a Design with
|
||||
Modifications, that is offered for sale, or otherwise made
|
||||
available for commercial use, whether as Hardware or in other forms.
|
||||
|
||||
“Source” means the format of the Design that is preferred for
|
||||
making Modifications, including all the files that are provided for
|
||||
that purpose.
|
||||
|
||||
“You” means any person or entity that exercises any right under
|
||||
this Licence.
|
||||
|
||||
2. What this Licence covers
|
||||
|
||||
The Licensor grants you a worldwide, royalty-free, non-exclusive,
|
||||
perpetual licence to:
|
||||
|
||||
a) use the Design and any associated Documentation,
|
||||
b) make Modifications,
|
||||
c) communicate the Design and any Documentation, including by
|
||||
making it available for copying and modification,
|
||||
d) manufacture and distribute Hardware that implements the Design
|
||||
or any part thereof,
|
||||
e) sell Products, and
|
||||
f) otherwise deal in the Design or any part thereof,
|
||||
|
||||
all in each case subject to the conditions set out in this Licence.
|
||||
|
||||
3. Conditions
|
||||
|
||||
3.1. Availability of Source
|
||||
|
||||
The Design and any Documentation that is made available under this
|
||||
Licence must be made available in Source format. If the Design is
|
||||
available in a format that is not Source, you must also make a
|
||||
version of the Design in Source format available, and you must
|
||||
inform others of how to obtain it.
|
||||
|
||||
3.2. Retention of copyright and other notices
|
||||
|
||||
You must retain all copyright and other notices that appear in the
|
||||
Design or Documentation.
|
||||
|
||||
3.3. Notice of Modifications
|
||||
|
||||
If you make Modifications, you must include a notice in any
|
||||
Distribution of the Design or Documentation, stating that
|
||||
Modifications have been made, by whom, and in what they consist.
|
||||
|
||||
3.4. Distribution of Modifications
|
||||
|
||||
If you distribute any Design or Documentation, whether unchanged or
|
||||
with Modifications, you must distribute it under this Licence, and
|
||||
you must inform the recipient that the Design or Documentation is
|
||||
licensed under this Licence.
|
||||
|
||||
3.5. No implicit approval
|
||||
|
||||
The Licensor shall not be taken to have approved any Modifications
|
||||
unless the Licessor expressly approves them in writing.
|
||||
|
||||
4. Disclaimer of warranty
|
||||
|
||||
THE DESIGN AND ANY DOCUMENTATION IS PROVIDED "AS IS" AND THE
|
||||
LICENSOR MAKES NO WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT
|
||||
LIMITED TO WARRANTIES OF MERCHANTABILITY, SATISFACTORY QUALITY,
|
||||
NON-INFRINGEMENT, TITLE, AND FITNESS FOR A PARTICULAR PURPOSE. THE
|
||||
DESIGN MAY CONTAIN DESIGN DEFECTS OR ERRORS THAT MAKE IT
|
||||
UNSUITABLE FOR USE IN CERTAIN APPLICATIONS. THE ENTIRE RISK AS TO
|
||||
ITS QUALITY AND PERFORMANCE REMAINS WITH YOU.
|
||||
|
||||
5. Limitation of liability
|
||||
|
||||
TO THE EXTENT NOT PROHIBITED BY LAW, IN NO EVENT WILL THE LICENSOR
|
||||
BE LIABLE FOR ANY DAMAGES WHATSOEVER, INCLUDING WITHOUT LIMITATION
|
||||
LOST PROFITS, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES, OR ANY
|
||||
OTHER LOSSES ARISING OUT OF THE USE OF THE DESIGN OR ANY
|
||||
DOCUMENTATION, HOWEVER CAUSED, WHETHER UNDER THEORY OF CONTRACT,
|
||||
TORT (INCLUDING NEGLIGENCE), PRODUCTS LIABILITY, OR OTHERWISE,
|
||||
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
|
||||
|
||||
6. Termination
|
||||
|
||||
This Licence is perpetual except that if you fail to comply with
|
||||
its terms, your rights under this Licence will terminate immediately.
|
||||
If you bring or defend any litigation or other claim alleging that
|
||||
the Design or any part thereof infringes any intellectual property
|
||||
right, all your rights under this Licence granted by the Licensor
|
||||
will terminate forty-five days after you bring or defend such claim.
|
||||
|
||||
7. General
|
||||
|
||||
This Licence shall be governed by the laws of the country of the
|
||||
Licensor, without reference to its conflict of laws provisions.
|
||||
If any part of this Licence is held to be invalid or unenforceable,
|
||||
the remaining parts shall remain in full force and effect.
|
||||
|
||||
---------------------------------------------------------------------
|
||||
|
||||
CERN-OHL-P is a variant of the CERN-OHL family. For more information
|
||||
see www.ohwr.org/cernohl.
|
||||
283
README.md
283
README.md
@@ -1,260 +1,55 @@
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR/stargazers)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR/network/members)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR/watchers)
|
||||
[](https://opensource.org/licenses/MIT)
|
||||
[](https://ohwr.org/cern_ohl_p_v2.txt)
|
||||
# AERIS-10
|
||||
|
||||
# AERIS-10: Open Source Pulse Linear Frequency Modulated Phased Array Radar
|
||||
10.5 GHz pulse linear-FM phased-array radar, with Overwatch beside the hardware design.
|
||||
|
||||
[](https://ohwr.org/cern_ohl_p_v2.txt)
|
||||
[](https://opensource.org/licenses/MIT)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR/issues)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR)
|
||||
[](https://github.com/NawfalMotii79/PLFM_RADAR/pulls)
|
||||
| | |
|
||||
|---|---|
|
||||
| [simulation/](simulation/) | Overwatch. Physics, radar model, and the operator view. |
|
||||
| [schematics/](schematics/) | Board schematics, stack-up, and production Gerbers. |
|
||||
| [parts/parts-list.csv](parts/parts-list.csv) | Parts list taken from the production bills of materials. |
|
||||
| [parts/cost.md](parts/cost.md) | Cost: the bills of materials do not contain prices. |
|
||||
|
||||
## Radar
|
||||
|
||||
AERIS-10 is an open-source, low-cost 10.5 GHz phased array radar system featuring Pulse Linear Frequency Modulated (LFM) modulation. Available in two versions (3km and 20km range), it's designed for researchers, drone developers, and serious SDR enthusiasts who want to explore and experiment with phased array radar technology.
|
||||
Two builds share the signal chain.
|
||||
|
||||
|
||||
- **AERIS-10N**, about 3 km, 8×16 patch array.
|
||||
- **AERIS-10X**, about 20 km, 32×16 dielectric-filled slotted waveguide, plus a 10 W QPA2962 stage on each of 16 amplifier boards.
|
||||
|
||||
## 📡 Overview
|
||||
The array steers elevation electronically over ±45°. Azimuth on the prototype is a stepper. Carrier is 10.5 GHz. The FPGA does chirp generation, pulse compression, Doppler, MTI, and CFAR. An STM32 sequences power, clocks, phase shifters, and the GPS/IMU.
|
||||
|
||||
The AERIS-10 project aims to democratize radar technology by providing a fully open-source, modular, and hackable radar system. Whether you're a university researcher, a drone startup, or an advanced maker, AERIS-10 offers a platform for experimenting with beamforming, pulse compression, Doppler processing, and target tracking.
|
||||
|
||||
|
||||
## 🔬 Key Features
|
||||
|
||||
|
||||
- **Open Source Hardware & Software** - Complete schematics, PCB layouts, firmware, and software available
|
||||
- **Dual Version Availability**:
|
||||
- **AERIS-10N (Nexus)**: 3km range with 8x16 patch antenna array
|
||||
- **AERIS-10E (Extended)**: 20km range with 32x16 dielectric-filled slotted waveguide array
|
||||
- **Full Electronic Beam Steering** - ±45° electronic steering in elevation and azimuth
|
||||
- **Advanced Signal Processing** - On-board FPGA handles pulse compression, Doppler FFT, MTI, and CFAR
|
||||
- **Python GUI** - User-friendly interface with map integration
|
||||
- **GPS/IMU Integration** - Real-time position and attitude correction
|
||||
- **Modular Design** - Separate power management, frequency synthesis, and RF boards
|
||||
### Where the rest of the design lives
|
||||
|
||||
## 🏗️ System Architecture
|
||||
| Path | Contents |
|
||||
|---|---|
|
||||
| `schematics/4_6_Schematics/` | Main board, power board, frequency synthesizer, power amplifier, antennas |
|
||||
| `schematics/4_7_Production Files/` | Gerbers, pick-and-place, and the original board bills of materials |
|
||||
| `parts/datasheets/` | Component datasheets |
|
||||
| `hardware/diagrams/` | Functional diagram |
|
||||
| `hardware/power/` | Power-rail and sequencing sheet |
|
||||
| `hardware/description/` | Project description |
|
||||
| `hardware/rf-simulations/` | Antenna, filter, and waveguide studies. These are circuit simulations, not Overwatch. |
|
||||
| `hardware/mechanical/` | Enclosure and waveguide drawings |
|
||||
| `hardware/application-notes/` | ADAR1000 user guide |
|
||||
| `firmware/` | STM32, FPGA, and the radar GUI |
|
||||
| `docs/` | Architecture notes, bring-up, and reports |
|
||||
|
||||
|
||||
## Overwatch
|
||||
|
||||
### Hardware Components
|
||||
From the repository root:
|
||||
|
||||
The AERIS-10 main sub-systems are:
|
||||
```bash
|
||||
cd simulation
|
||||
python -m pip install -e '.[server,dev,flight,performance]'
|
||||
python -m defense_sim.cli --steps 500 --seed 7
|
||||
```
|
||||
|
||||
- **Power Management Board** - Supplies all necessary voltage levels to the electronics components with proper filtering and sequencing (sequencing ensured by the microcontroller)
|
||||
The earlier simulator, the standalone scripts and the single-page viewer that used to sit in the repository root, has been removed.
|
||||
|
||||
- **Frequency Synthesizer Board** - Uses a high-performance Low Jitter Clock Generator (AD9523-1) that supplies phase-aligned clock references for:
|
||||
- RX and TX Frequency Synthesizers (ADF4382)
|
||||
- DAC
|
||||
- ADC
|
||||
- FPGA
|
||||
## Contact
|
||||
|
||||
- **Main Board** containing:
|
||||
- **DAC** - Generates the RADAR Chirps
|
||||
- **2x Microwave Mixers (LTC5552)** - For up-conversion and IF-down-conversion
|
||||
- **4x 4-Channel Phase Shifters (ADAR1000)** - For RX and TX chain beamforming
|
||||
- **16x Front End Chips (ADTR1107)** - Used for both Low Noise Amplifying (RX) and Power Amplifying (TX) stages
|
||||
- **XC7A50T FPGA** - Handles RADAR Signal Processing on the upstream FTG256 board:
|
||||
- PLFM Chirps generation via the DAC
|
||||
- Raw ADC data read
|
||||
- Hybrid Automatic Gain Control (AGC) — cross-layer FPGA/STM32/GUI loop
|
||||
- I/Q Baseband Down-Conversion
|
||||
- Decimation
|
||||
- Filtering
|
||||
- Forward FFT
|
||||
- Pulse Compression
|
||||
- Doppler, MTI and CFAR processing
|
||||
- USB Interface
|
||||
- **STM32F746xx Microcontroller** - Used for:
|
||||
- Power-up and power-down sequencing (see Power Management Excel File)
|
||||
- FPGA communication
|
||||
- Setup and Interface with:
|
||||
- Clock Generator (AD9523-1)
|
||||
- 2x Frequency Synthesizers (ADF4382)
|
||||
- 4x 4-Channel Phase Shifters (ADAR1000) for RADAR pulse sequencing
|
||||
- 2x ADS7830 8-channel I²C ADCs (Main Board, U88 @ 0x48 / U89 @ 0x4A) for 16x Idq measurement, one per PA channel, each sensed through a 5 mΩ shunt on the PA board and an INA241A3 current-sense amplifier (x50) on the Main Board
|
||||
- 2x DAC5578 8-channel I²C DACs (Main Board, U7 @ 0x48 / U69 @ 0x49) for 16x Vg control, one per PA channel; closed-loop calibrated at boot to the target Idq
|
||||
- GPS module (UM982) for GUI map centering and per-detection position tagging
|
||||
- GY-85 IMU for pitch/roll correction of target coordinates
|
||||
- BMP180 Barometer
|
||||
- Stepper Motor
|
||||
- 1x ADS7830 8-channel I²C ADC (Main Board, U10) reading 8 thermistors for thermal monitoring; a single GPIO (EN_DIS_COOLING) switches the cooling fans on when any channel exceeds the threshold
|
||||
- RF switches
|
||||
|
||||
- **16x Power Amplifier Boards** - Used only for AERIS-10E version, featuring 10Watt QPA2962 GaN amplifier for extended range
|
||||
|
||||
- **Antenna Arrays**:
|
||||
- **AERIS-10N (Nexus)** - 8x16 patch antenna array
|
||||
- **AERIS-10X (Extended)** - 32x16 dielectric-filled slotted waveguide antenna array
|
||||
|
||||
- **Miscellaneous Components**:
|
||||
- Slip-Ring
|
||||
- Stepper Motor and drivers
|
||||
- Cooling Fans
|
||||
- Enclosure
|
||||
|
||||
### Processing Pipeline
|
||||
|
||||
1. **Waveform Generation** - DAC creates LFM chirps
|
||||
2. **Up/Down Conversion** - LTC5552 mixers handle frequency translation
|
||||
3. **Beam Steering** - ADAR1000 phase shifters control 16 elements
|
||||
4. **Signal Processing (FPGA)**:
|
||||
- Raw ADC data capture
|
||||
- I/Q baseband down-conversion
|
||||
- Decimation & filtering (CIC/FIR)
|
||||
- Pulse compression
|
||||
- Doppler FFT processing
|
||||
- MTI & CFAR detection
|
||||
5. **System Management (STM32)**:
|
||||
- Power sequencing
|
||||
- Peripheral configuration
|
||||
- GPS/IMU integration
|
||||
- Stepper motor control
|
||||
6. **Visualization (Python GUI)**:
|
||||
- Real-time target plotting
|
||||
- Map integration
|
||||
- Radar control interface
|
||||
|
||||
|
||||
<!-- V6 GIF removed — V6 is deprecated. V65 Tk and V7 PyQt6 are the active GUIs. -->
|
||||
|
||||
## 📊 Technical Specifications
|
||||
|
||||
| Parameter | AERIS-10N (Nexus) | AERIS-10X (Extended) |
|
||||
|-----------|-------------------|----------------------|
|
||||
| **Frequency** | 10.5 GHz | 10.5 GHz |
|
||||
| **Max Range** | 3 km | 20 km |
|
||||
| **Antenna** | 8x16 Patch Array | 32x16 Slotted Waveguide |
|
||||
| **Beam Steering** | Electronic (±45°) | Electronic (±45°) |
|
||||
| **Mechanical Scan** | 360° (stepper motor) | 360° (stepper motor) |
|
||||
| **Output Power** | ~1W×16 | 10W×16 (GaN amplifier) |
|
||||
| **Processing** | FPGA + STM32 | FPGA + STM32 |
|
||||
|
||||
## 🚀 Getting Started
|
||||
|
||||
### 🧹 Repository File Placement Policy
|
||||
|
||||
To keep the repository root clean and make artifacts easy to find, place generated files in the following locations:
|
||||
|
||||
- **Published reports (tracked, GitHub Pages):** `docs/`
|
||||
- Example: `docs/AERIS_Simulation_Report_v2.pdf`
|
||||
- **Simulation-generated outputs (local, gitignored):** `5_Simulations/generated/`
|
||||
- Plots, scenario outputs, temporary analysis directories
|
||||
- **FPGA/Vivado generated artifacts (local, gitignored):** `9_Firmware/9_2_FPGA/reports/`
|
||||
- VCD/VVP dumps, temporary CSVs, local report snapshots
|
||||
- **Reusable FPGA automation scripts (tracked):** `9_Firmware/9_2_FPGA/scripts/`
|
||||
- TCL flows, helper scripts used by build/bring-up
|
||||
|
||||
**Do not leave generated artifacts in the repository root.**
|
||||
|
||||
### Prerequisites
|
||||
|
||||
- Basic understanding of radar principles
|
||||
- Experience with PCB assembly (for hardware build)
|
||||
- Python 3.8+ for the GUI software
|
||||
- FPGA development tools (Vivado) for signal processing modifications
|
||||
|
||||
### Hardware Assembly
|
||||
|
||||
1. **Order PCBs**: Production outputs are under `/4_Schematics and Boards Layout/4_7_Production Files`
|
||||
2. **Source Components**: BOM/CPL files are co-located under `/4_Schematics and Boards Layout/4_7_Production Files`
|
||||
3. **Assembly**: Use the schematics in `/4_Schematics and Boards Layout/4_6_Schematics` together with the production outputs above; a standalone assembly guide is not currently tracked
|
||||
4. **Antenna**: Choose appropriate array files for your target variant
|
||||
5. **Enclosure**: Mechanical drawings currently live in `/8_Utils/Mechanical_Drawings`
|
||||
|
||||
## 📜 License
|
||||
|
||||
This project is open-source but uses **different licenses for hardware and software** to ensure proper legal coverage.
|
||||
|
||||
### Hardware Documentation
|
||||
The hardware design files—including:
|
||||
- Schematics and PCB layouts (in `/4_Schematics and Boards Layout`)
|
||||
- Bill of Materials (BOM) files
|
||||
- Gerber files and manufacturing outputs
|
||||
- Mechanical drawings and enclosure designs
|
||||
|
||||
are licensed under the **CERN Open Hardware Licence Version 2 – Permissive (CERN-OHL-P)** .
|
||||
|
||||
This is a hardware-specific license that:
|
||||
- ✅ Clearly defines "Hardware," "Documentation," and "Products"
|
||||
- ✅ Includes explicit patent protection for contributors and users
|
||||
- ✅ Provides stronger liability limitations (important for high-power RF)
|
||||
- ✅ Aligns with professional open-hardware standards (CERN, OSHWA)
|
||||
|
||||
You may use, modify, and sell products based on these designs, provided you:
|
||||
- Maintain the original copyright notices
|
||||
- Distribute any modified designs under the same license
|
||||
- Make your modifications available in Source format
|
||||
|
||||
### Software and Firmware
|
||||
The software components—including:
|
||||
- FPGA code (VHDL/Verilog in `/9_Firmware`)
|
||||
- Microcontroller firmware (STM32)
|
||||
- Python GUI and utilities
|
||||
|
||||
remain under the **MIT License** for maximum flexibility.
|
||||
|
||||
### Full License Texts
|
||||
- The complete CERN-OHL-P license text is in the `LICENSE` file
|
||||
- MIT license terms apply to software where not otherwise specified
|
||||
|
||||
### Why This Change?
|
||||
Originally, the entire project used the MIT license. The community (special thanks to gmaynez!) pointed out that MIT lacks legal protections needed for physical hardware. The switch to CERN-OHL-P ensures the project is properly protected while maintaining the same permissive spirit.
|
||||
|
||||
## 📚 Documentation
|
||||
|
||||
Comprehensive documentation is available in the `/docs` folder and served via GitHub Pages at [https://NawfalMotii79.github.io/PLFM_RADAR/docs/](https://NawfalMotii79.github.io/PLFM_RADAR/docs/):
|
||||
|
||||
- [System Architecture](/docs/architecture.html)
|
||||
- [Implementation Log](/docs/implementation-log.html)
|
||||
- [Hardware Bring-Up Guide](/docs/bring-up.html)
|
||||
- [Test Reports](/docs/reports.html)
|
||||
- [Release Notes](/docs/release-notes.html)
|
||||
|
||||
## 🤝 Contributing
|
||||
|
||||
We welcome contributions! Please see our [Contributing Guidelines](/CONTRIBUTING.md) for details on repo layout, branch workflow, and basic PR checks.
|
||||
|
||||
Areas where help is especially appreciated:
|
||||
- **RF Engineers**: Review designs, optimize antenna performance
|
||||
- **FPGA Developers**: Optimize signal processing pipeline
|
||||
- **Software Developers**: Enhance Python GUI and SDK
|
||||
- **Beta Testers**: University researchers, drone startups, advanced makers
|
||||
|
||||
## 📞 Contact & Collaboration
|
||||
|
||||
I welcome serious inquiries from researchers, engineers, and potential collaborators. However, due to the high volume of interest in this project, please understand that I cannot guarantee a response to every message.
|
||||
|
||||
- **Technical questions or bug reports**: Please [open a GitHub issue](https://github.com/NawfalMotii79/PLFM_RADAR/issues) so the whole community can benefit from the discussion.
|
||||
- **Collaboration, licensing, or business inquiries**: 📧 nawfal.motii.33 [at] gmail [dot] com
|
||||
|
||||
## 💰 Sponsors
|
||||
|
||||
|
||||
|
||||
---
|
||||
|
||||
**Star ⭐ this repository if you're interested in open-source radar technology!**
|
||||
|
||||
*Note: This is an active development project. Some features are still in progress. Check the issues page for known limitations and upcoming features.*
|
||||
|
||||
## 19,000 stars – Thank you
|
||||
|
||||
This project started in a small workshop in Morocco. Today, 19,000 engineers on GitHub have starred it.
|
||||
|
||||
I am genuinely humbled.
|
||||
|
||||
**What this tells me:**
|
||||
|
||||
- Open source radar matters
|
||||
- Affordable sensing is needed
|
||||
- Engineers want to build, not just buy
|
||||
|
||||
Thank you to everyone who starred, forked, opened issues, submitted PRs, and shared this project.
|
||||
|
||||
The work continues.
|
||||
|
||||
**Nawfal Motii**
|
||||
ABAC INDUSTRY
|
||||
(http://www.abacindustry.com)
|
||||
Nawfal Motii, ABAC INDUSTRY. Technical questions go through the project issue tracker. Collaboration and business inquiries: nawfal.motii.33 [at] gmail [dot] com.
|
||||
|
||||
57
debug2.js
57
debug2.js
@@ -1,57 +0,0 @@
|
||||
const fs = require('fs');
|
||||
const html = fs.readFileSync('defense_sim.html', 'utf8');
|
||||
const script = html.match(/<script>([\s\S]*?)<\/script>/)[1];
|
||||
function extractFn(src, name) {
|
||||
const re = new RegExp('function ' + name + '\\(', 'g');
|
||||
const m = re.exec(src);
|
||||
if (!m) return null;
|
||||
let i = m.index, depth = 0, start = src.indexOf('{', i);
|
||||
for (let j = start; j < src.length; j++) {
|
||||
if (src[j] === '{') depth++;
|
||||
else if (src[j] === '}') { depth--; if (depth === 0) return src.slice(i, j + 1); }
|
||||
}
|
||||
return null;
|
||||
}
|
||||
const fns = ['targetTrajectory','detectFromRadars','estimatePosition3D','estimateTrajectory','computeInterceptPoint','perpDistToLine','gaussianNoise','kmToPixels'];
|
||||
eval(fns.map(n => extractFn(script, n)).filter(Boolean).join('\n'));
|
||||
const MEASUREMENT_ERROR = 20.0, IR_ERROR = 3.0, SCALE = 10;
|
||||
|
||||
// Single prop-drone trial, fixed geometry for reproducibility
|
||||
const P = { startDist:15, speed:30, cruiseMin:200, cruiseMax:200, descentStart:3, terminalAlt:0 };
|
||||
const numRadars=8, numDrones=24, radarRange=20, droneSpeed=100, irRange=5, cityRadius=10, droneLine=8, dt=0.5;
|
||||
const startAngle = Math.PI/3; // fixed
|
||||
const startPos = [P.startDist*Math.cos(startAngle), P.startDist*Math.sin(startAngle)];
|
||||
const targetAngle = Math.PI/3; // straight radial approach
|
||||
const targetDist = 2;
|
||||
const targetPos = [targetDist*Math.cos(targetAngle), targetDist*Math.sin(targetAngle)];
|
||||
const cruiseAlt = 200;
|
||||
const trajectory = targetTrajectory(startPos, targetPos, cruiseAlt, P.speed, P.descentStart, P.terminalAlt);
|
||||
const radarPositions = [];
|
||||
for (let i=0;i<numRadars;i++){const a=2*Math.PI*i/numRadars,d=cityRadius+2;radarPositions.push([d*Math.cos(a),d*Math.sin(a),0]);}
|
||||
|
||||
let time=0, detectionTime=null, observations=[];
|
||||
function getTruePos(t){let b=trajectory[0],bd=1e18;for(const p of trajectory){const d=Math.abs(p.t-t);if(d<bd){bd=d;b=p;}}return b;}
|
||||
|
||||
for (let step=0; step<20 && detectionTime===null; step++) {
|
||||
time += dt;
|
||||
const tp = getTruePos(time);
|
||||
const pos = [tp.x, tp.y, tp.z/1000];
|
||||
const detections = detectFromRadars(pos, radarPositions, radarRange, gaussianNoise);
|
||||
if (detections.length>0) observations.push({x:pos[0],y:pos[1],z:pos[2]*1000,t:time});
|
||||
if (observations.length>=3) {
|
||||
detectionTime = time;
|
||||
const posEst = estimatePosition3D(detections, radarPositions);
|
||||
if (posEst) observations[observations.length-1] = {x:posEst[0]/1000,y:posEst[1]/1000,z:posEst[2],t:time};
|
||||
console.log(`=== detection at t=${time.toFixed(1)}s, ${observations.length} observations ===`);
|
||||
console.log('true target pos (km):', pos.map(v=>v.toFixed(3)));
|
||||
if (posEst) console.log('posEst (km):', posEst.map(v=>v.toFixed(3)));
|
||||
console.log('observations (x km, y km, z m, t s):');
|
||||
for (const o of observations) console.log(' ', o.x.toFixed(4), o.y.toFixed(4), o.z.toFixed(2), o.t.toFixed(1));
|
||||
const traj = estimateTrajectory(observations);
|
||||
const velMs = [traj.vel[0]*1000, traj.vel[1]*1000, traj.vel[2]];
|
||||
console.log('vel (raw):', traj.vel.map(v=>v.toExponential(3)));
|
||||
console.log('velMs:', velMs.map(v=>v.toFixed(2)));
|
||||
console.log('speedEst:', Math.sqrt(velMs[0]**2+velMs[1]**2+velMs[2]**2).toFixed(1), 'R2:', traj.r2.toFixed(3));
|
||||
break;
|
||||
}
|
||||
}
|
||||
162
debug_sim.js
162
debug_sim.js
@@ -1,162 +0,0 @@
|
||||
const fs = require('fs');
|
||||
const html = fs.readFileSync('defense_sim.html', 'utf8');
|
||||
const script = html.match(/<script>([\s\S]*?)<\/script>/)[1];
|
||||
|
||||
function extractFn(src, name) {
|
||||
const re = new RegExp('function ' + name + '\\(', 'g');
|
||||
const m = re.exec(src);
|
||||
if (!m) return null;
|
||||
let i = m.index, depth = 0, start = src.indexOf('{', i);
|
||||
for (let j = start; j < src.length; j++) {
|
||||
if (src[j] === '{') depth++;
|
||||
else if (src[j] === '}') { depth--; if (depth === 0) return src.slice(i, j + 1); }
|
||||
}
|
||||
return null;
|
||||
}
|
||||
function extractConst(src, name) {
|
||||
const m = new RegExp('const ' + name + ' = ([^;\\n]+);').exec(src);
|
||||
return m ? m[1].trim() : null;
|
||||
}
|
||||
|
||||
const fns = ['targetTrajectory', 'detectFromRadars', 'estimatePosition3D', 'estimateTrajectory',
|
||||
'computeInterceptPoint', 'perpDistToLine', 'gaussianNoise', 'kmToPixels'];
|
||||
eval(fns.map(n => extractFn(script, n)).filter(Boolean).join('\n'));
|
||||
const MEASUREMENT_ERROR = parseFloat(extractConst(script, 'MEASUREMENT_ERROR'));
|
||||
const IR_ERROR = parseFloat(extractConst(script, 'IR_ERROR'));
|
||||
const SCALE = parseFloat(extractConst(script, 'SCALE'));
|
||||
console.log('MEASUREMENT_ERROR', MEASUREMENT_ERROR, 'IR_ERROR', IR_ERROR, 'SCALE', SCALE);
|
||||
|
||||
function runTrial(preset, verbose) {
|
||||
const P = {
|
||||
'shahed-136': { startDist:30, speed:50, cruiseMin:100, cruiseMax:300, descentStart:5, terminalAlt:0 },
|
||||
'shahed-238': { startDist:35, speed:150, cruiseMin:1000, cruiseMax:2000, descentStart:8, terminalAlt:200 },
|
||||
'prop-drone': { startDist:15, speed:30, cruiseMin:100, cruiseMax:300, descentStart:3, terminalAlt:0 },
|
||||
}[preset];
|
||||
const numRadars = 8, numDrones = 24, radarRange = 20, droneSpeed = 100, irRange = 5, cityRadius = 10, droneLine = 8;
|
||||
const dt = 0.5;
|
||||
|
||||
const startAngle = Math.random() * 2 * Math.PI;
|
||||
const startPos = [P.startDist * Math.cos(startAngle), P.startDist * Math.sin(startAngle)];
|
||||
const targetAngle = Math.random() * 2 * Math.PI;
|
||||
const targetDist = Math.random() * cityRadius * 0.8;
|
||||
const targetPos = [targetDist * Math.cos(targetAngle), targetDist * Math.sin(targetAngle)];
|
||||
const cruiseAlt = P.cruiseMin + Math.random() * (P.cruiseMax - P.cruiseMin);
|
||||
|
||||
const trajectory = targetTrajectory(startPos, targetPos, cruiseAlt, P.speed, P.descentStart, P.terminalAlt);
|
||||
|
||||
const radarPositions = [];
|
||||
for (let i = 0; i < numRadars; i++) {
|
||||
const a = 2 * Math.PI * i / numRadars, d = cityRadius + 2;
|
||||
radarPositions.push([d * Math.cos(a), d * Math.sin(a), 0]);
|
||||
}
|
||||
const dronePositions = [];
|
||||
for (let i = 0; i < numDrones; i++) {
|
||||
const a = 2 * Math.PI * i / numDrones, d = cityRadius + droneLine;
|
||||
dronePositions.push({ x: d * Math.cos(a), y: d * Math.sin(a), z: 0, state: 'waiting', intercept: null });
|
||||
}
|
||||
|
||||
let time = 0, detectionTime = null, launched = [], collision = false, impact = false, observations = [];
|
||||
|
||||
function getTruePos(t) {
|
||||
let best = trajectory[0], bd = Infinity;
|
||||
for (const p of trajectory) { const d = Math.abs(p.t - t); if (d < bd) { bd = d; best = p; } }
|
||||
return best;
|
||||
}
|
||||
|
||||
while (time < 1200) {
|
||||
time += dt;
|
||||
const tp = getTruePos(time);
|
||||
const pos = [tp.x, tp.y, tp.z / 1000];
|
||||
|
||||
const detections = detectFromRadars(pos, radarPositions, radarRange, gaussianNoise);
|
||||
if (detections.length > 0) observations.push({ x: pos[0], y: pos[1], z: pos[2] * 1000, t: time });
|
||||
|
||||
if (detectionTime === null && observations.length >= 3) {
|
||||
detectionTime = time;
|
||||
const posEst = estimatePosition3D(detections, radarPositions);
|
||||
if (posEst) observations[observations.length - 1] = { x: posEst[0]/1000, y: posEst[1]/1000, z: posEst[2], t: time };
|
||||
const traj = estimateTrajectory(observations);
|
||||
if (traj) {
|
||||
const velMs = [traj.vel[0]*1000, traj.vel[1]*1000, traj.vel[2]];
|
||||
const speedEst = Math.sqrt(velMs[0]**2 + velMs[1]**2 + velMs[2]**2);
|
||||
const horizDist = Math.sqrt(pos[0]**2 + pos[1]**2);
|
||||
const tti = horizDist / (speedEst / 1000);
|
||||
const velKs = [velMs[0]/1000, velMs[1]/1000, velMs[2]/1000];
|
||||
const posM = [pos[0]*1000, pos[1]*1000, pos[2]*1000];
|
||||
const cands = [];
|
||||
for (let i = 0; i < numDrones; i++) {
|
||||
const dp = dronePositions[i];
|
||||
const dpM = [dp.x*1000, dp.y*1000, dp.z*1000];
|
||||
const r = computeInterceptPoint(dpM, posM, velMs, droneSpeed);
|
||||
const dist = perpDistToLine(pos, velKs, [dp.x, dp.y, dp.z]);
|
||||
const viable = (r.t > 0 && r.t < tti);
|
||||
cands.push({ i, t: r.t, intercept: [r.intercept[0]/1000, r.intercept[1]/1000, r.intercept[2]/1000], dist, viable });
|
||||
}
|
||||
const vi = cands.filter(c => c.viable).sort((a,b)=>a.dist-b.dist);
|
||||
const nv = cands.filter(c=>!c.viable).sort((a,b)=>a.dist-b.dist);
|
||||
const chosen = vi.slice(0,3);
|
||||
for (const c of nv) { if (chosen.length >= 3) break; chosen.push(c); }
|
||||
launched = chosen.map(c=>c.i);
|
||||
for (const c of chosen) { dronePositions[c.i].intercept = c.intercept; dronePositions[c.i].state = 'flying'; }
|
||||
if (verbose) console.log(`[detect t=${time.toFixed(1)}s] horiz=${horizDist.toFixed(1)}km speedEst=${speedEst.toFixed(0)} tti=${tti.toFixed(0)}s viable=${vi.length} chosen=${launched.join(',')} R2=${traj.r2.toFixed(3)}`);
|
||||
}
|
||||
}
|
||||
|
||||
if (detectionTime !== null && !collision && observations.length >= 2) {
|
||||
const traj = estimateTrajectory(observations);
|
||||
if (traj) {
|
||||
const velMs = [traj.vel[0]*1000, traj.vel[1]*1000, traj.vel[2]];
|
||||
const speedEst = Math.sqrt(velMs[0]**2+velMs[1]**2+velMs[2]**2);
|
||||
const horizDist = Math.sqrt(pos[0]**2+pos[1]**2);
|
||||
const tti = horizDist / (speedEst/1000);
|
||||
if ((time - detectionTime) % 5 < dt) {
|
||||
const posM = [pos[0]*1000, pos[1]*1000, pos[2]*1000];
|
||||
for (const i of launched) {
|
||||
const dr = dronePositions[i];
|
||||
if (dr.state === 'flying') {
|
||||
const dpM = [dr.x*1000, dr.y*1000, dr.z*1000];
|
||||
const r = computeInterceptPoint(dpM, posM, velMs, droneSpeed);
|
||||
dr.intercept = [r.intercept[0]/1000, r.intercept[1]/1000, r.intercept[2]/1000];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (launched.length > 0 && !collision) {
|
||||
for (const i of launched) {
|
||||
const dr = dronePositions[i];
|
||||
if (dr.state === 'flying' || dr.state === 'ir-tracking') {
|
||||
const it = dr.intercept;
|
||||
if (it) {
|
||||
const d = [it[0]-dr.x, it[1]-dr.y, it[2]-dr.z];
|
||||
const L = Math.sqrt(d[0]**2+d[1]**2+d[2]**2);
|
||||
if (L > 0.001) { const s = droneSpeed*dt/1000; dr.x += d[0]/L*s; dr.y += d[1]/L*s; dr.z += d[2]/L*s; }
|
||||
}
|
||||
const dx = pos[0]-dr.x, dy = pos[1]-dr.y, dz = pos[2]-dr.z;
|
||||
const dist3D = Math.sqrt(dx*dx+dy*dy+dz*dz);
|
||||
if (dr.state === 'flying' && dist3D < irRange) dr.state = 'ir-tracking';
|
||||
if (dist3D*1000 < 20) { collision = true; if (verbose) console.log(`[COLLISION t=${time.toFixed(1)}s drone ${i}]`); }
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
const horizDist = Math.sqrt(pos[0]**2 + pos[1]**2);
|
||||
if (horizDist < cityRadius && !collision) { impact = true; if (verbose) console.log(`[IMPACT t=${time.toFixed(1)}s]`); break; }
|
||||
if (collision) break;
|
||||
}
|
||||
return { collision, impact, detectionTime, observations: observations.length, launched: launched.length, time };
|
||||
}
|
||||
|
||||
console.log('\n=== prop-drone verbose trials ===');
|
||||
for (let k = 0; k < 3; k++) {
|
||||
console.log(`--- trial ${k} ---`);
|
||||
const r = runTrial('prop-drone', true);
|
||||
console.log(` => collision=${r.collision} impact=${r.impact} detect=${r.detectionTime ? r.detectionTime.toFixed(1)+'s' : 'none'} obs=${r.observations} launched=${r.launched}`);
|
||||
}
|
||||
console.log('\n=== summary (30 trials each) ===');
|
||||
for (const preset of ['shahed-136','shahed-238','prop-drone']) {
|
||||
let c = 0, im = 0;
|
||||
for (let k = 0; k < 30; k++) { const r = runTrial(preset, false); if (r.collision) c++; else if (r.impact) im++; }
|
||||
console.log(`${preset.padEnd(14)} collision=${c}/30 (${(c/30*100).toFixed(0)}%) impact=${im}/30`);
|
||||
}
|
||||
1368
defense_sim.html
1368
defense_sim.html
File diff suppressed because it is too large
Load Diff
701
defense_sim.py
701
defense_sim.py
@@ -1,701 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Generalized Projectile Defense Simulation
|
||||
==========================================
|
||||
Models a projectile/drone attack profile and defender radar/drone
|
||||
intercept with proper physics and error propagation.
|
||||
|
||||
Supports multiple projectile types via presets:
|
||||
- Shahed-136 (Geran-2): prop, ~50 m/s, low altitude
|
||||
- Shahed-238 (Geran-3): jet, ~150 m/s, 1-2 km cruise (default)
|
||||
- Cruise Missile: ~300 m/s, sea-skimming
|
||||
- Ballistic (Re-entry): ~800 m/s, high altitude entry
|
||||
- Small UAS / Prop Drone: ~30 m/s, low altitude
|
||||
|
||||
Defender:
|
||||
- 8 AERIS-10 radars (20 km range, ±20 m single-shot accuracy)
|
||||
- 24 interceptor drones (100 m/s, IR homing within 5 km)
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import json
|
||||
import sys
|
||||
|
||||
# =============================================================================
|
||||
# Projectile Presets
|
||||
# =============================================================================
|
||||
# Each preset defines the attack profile:
|
||||
# start_dist: km - distance from city center at launch
|
||||
# speed: m/s - cruise speed
|
||||
# cruise_alt_min / cruise_alt_max: m - cruise altitude range (randomized within)
|
||||
# descent_start: km - distance from target at which descent begins
|
||||
# terminal_alt: m - altitude at impact
|
||||
PRESETS = {
|
||||
'shahed-136': {'start_dist': 30, 'speed': 50, 'cruise_alt_min': 100, 'cruise_alt_max': 300, 'descent_start': 5, 'terminal_alt': 0},
|
||||
'shahed-238': {'start_dist': 35, 'speed': 150, 'cruise_alt_min': 1000, 'cruise_alt_max': 2000, 'descent_start': 8, 'terminal_alt': 200},
|
||||
'cruise-missile': {'start_dist': 50, 'speed': 300, 'cruise_alt_min': 50, 'cruise_alt_max': 150, 'descent_start': 5, 'terminal_alt': 0},
|
||||
'ballistic-reentry':{'start_dist': 100, 'speed': 800, 'cruise_alt_min': 60000, 'cruise_alt_max': 80000, 'descent_start': 40, 'terminal_alt': 0},
|
||||
'prop-drone': {'start_dist': 15, 'speed': 30, 'cruise_alt_min': 100, 'cruise_alt_max': 300, 'descent_start': 3, 'terminal_alt': 0},
|
||||
}
|
||||
DEFAULT_PRESET = 'shahed-238'
|
||||
|
||||
|
||||
def get_preset(name=None):
|
||||
"""Return the preset dict for a given name, or the default."""
|
||||
if name is None:
|
||||
name = DEFAULT_PRESET
|
||||
return PRESETS.get(name, PRESETS[DEFAULT_PRESET])
|
||||
|
||||
# =============================================================================
|
||||
# Defender constants
|
||||
# =============================================================================
|
||||
|
||||
# Radar parameters
|
||||
RADAR_RANGE = 20000.0 # m
|
||||
RADAR_ERROR_1SIGMA = 20.0 # m (single measurement, 1 sigma)
|
||||
NUM_RADARS = 8
|
||||
|
||||
# Drone parameters
|
||||
DRONE_SPEED = 100.0 # m/s
|
||||
DRONE_MAX_ALT = 4000.0 # m
|
||||
DRONE_DEFENSE_LINE = 8000.0 # m from city center
|
||||
IR_RANGE = 5000.0 # m
|
||||
IR_ERROR_1SIGMA = 3.0 # m
|
||||
NUM_DRONES = 24
|
||||
CITY_RADIUS = 10000.0 # m
|
||||
|
||||
# Simulation
|
||||
DT = 0.5 # s (time step, coarser for speed)
|
||||
SIM_DURATION = 400.0 # s
|
||||
|
||||
# =============================================================================
|
||||
# Flight Profile
|
||||
# =============================================================================
|
||||
|
||||
def target_trajectory(start_pos, target_pos, cruise_alt, speed, descent_start, terminal_alt):
|
||||
"""
|
||||
Generate a realistic projectile flight trajectory.
|
||||
- Cruise at cruise_alt until `descent_start` distance from target
|
||||
- Gradual descent to `terminal_alt` over the remaining distance
|
||||
- Straight horizontal path (simplified for 2D projection)
|
||||
|
||||
descent_start: distance (m) from target at which descent begins
|
||||
terminal_alt: altitude (m) at impact
|
||||
|
||||
Returns: list of (x, y, z, t) tuples
|
||||
"""
|
||||
start_pos = np.array(start_pos)
|
||||
target_pos = np.array(target_pos)
|
||||
|
||||
# Horizontal direction
|
||||
dx = target_pos[0] - start_pos[0]
|
||||
dy = target_pos[1] - start_pos[1]
|
||||
horiz_dist = np.sqrt(dx**2 + dy**2)
|
||||
|
||||
if horiz_dist < 1e-6:
|
||||
return []
|
||||
|
||||
dir_x = dx / horiz_dist
|
||||
dir_y = dy / horiz_dist
|
||||
|
||||
# Descent begins at `descent_start` distance from target
|
||||
# So descent begins at (horiz_dist - descent_start) from start
|
||||
descent_start_dist = max(0.0, horiz_dist - descent_start)
|
||||
|
||||
# Descent distance (clamped to available distance)
|
||||
descent_dist = max(0.0, min(descent_start, horiz_dist))
|
||||
|
||||
# Total duration
|
||||
total_duration = horiz_dist / speed
|
||||
|
||||
trajectory = []
|
||||
t = 0.0
|
||||
|
||||
while t <= total_duration + DT:
|
||||
# Distance traveled
|
||||
dist_traveled = speed * t
|
||||
|
||||
# Determine altitude
|
||||
if descent_dist <= 0 or dist_traveled < descent_start_dist:
|
||||
# Cruise phase
|
||||
z = cruise_alt
|
||||
else:
|
||||
# Descent phase: linear interpolation from cruise_alt to terminal_alt
|
||||
descent_progress = min(1.0, (dist_traveled - descent_start_dist) / descent_dist)
|
||||
z = cruise_alt + (terminal_alt - cruise_alt) * descent_progress
|
||||
|
||||
# Horizontal position
|
||||
x = start_pos[0] + dir_x * dist_traveled
|
||||
y = start_pos[1] + dir_y * dist_traveled
|
||||
|
||||
trajectory.append((x, y, z, t))
|
||||
t += DT
|
||||
|
||||
return trajectory
|
||||
|
||||
# =============================================================================
|
||||
# Radar Detection Model
|
||||
# =============================================================================
|
||||
|
||||
def radar_positions(num_radars, city_radius):
|
||||
"""Generate radar positions on a circle around the city."""
|
||||
positions = []
|
||||
for i in range(num_radars):
|
||||
angle = 2 * np.pi * i / num_radars
|
||||
dist = city_radius + 2000 # 2 km outside city
|
||||
x = dist * np.cos(angle)
|
||||
y = dist * np.sin(angle)
|
||||
positions.append(np.array([x, y, 0.0]))
|
||||
return positions
|
||||
|
||||
def drone_positions(num_drones, city_radius, defense_line):
|
||||
"""Generate drone positions on a circle around the city."""
|
||||
positions = []
|
||||
for i in range(num_drones):
|
||||
angle = 2 * np.pi * i / num_drones
|
||||
dist = city_radius + defense_line
|
||||
x = dist * np.cos(angle)
|
||||
y = dist * np.sin(angle)
|
||||
positions.append({
|
||||
'pos': np.array([x, y, 0.0]),
|
||||
'state': 'waiting',
|
||||
'index': i
|
||||
})
|
||||
return positions
|
||||
|
||||
def detect_from_radars(pos, radar_pos_list, radar_range, rng):
|
||||
"""
|
||||
Detect target position from all radars in range.
|
||||
Returns list of (range, azimuth, elevation, radar_index) and estimated 3D position.
|
||||
"""
|
||||
detections = []
|
||||
for i, rpos in enumerate(radar_pos_list):
|
||||
dist = np.linalg.norm(pos - rpos)
|
||||
if dist < radar_range:
|
||||
# Add measurement noise
|
||||
range_meas = dist + rng.normal(0, RADAR_ERROR_1SIGMA)
|
||||
azimuth_true = np.arctan2(pos[1] - rpos[1], pos[0] - rpos[0])
|
||||
elevation_true = np.arctan2(pos[2] - rpos[2], np.sqrt((pos[0]-rpos[0])**2 + (pos[1]-rpos[1])**2))
|
||||
|
||||
azimuth_meas = azimuth_true + rng.normal(0, RADAR_ERROR_1SIGMA / max(100, dist))
|
||||
elevation_meas = elevation_true + rng.normal(0, RADAR_ERROR_1SIGMA / max(100, dist))
|
||||
|
||||
detections.append({
|
||||
'range': range_meas,
|
||||
'azimuth': azimuth_meas,
|
||||
'elevation': elevation_meas,
|
||||
'radar_index': i
|
||||
})
|
||||
|
||||
return detections
|
||||
|
||||
def estimate_position_3d(detections, radar_pos_list, rng):
|
||||
"""
|
||||
Estimate 3D position from multiple radar detections.
|
||||
Each radar provides azimuth, elevation, and range, which together define a 3D point.
|
||||
The estimated position is the least-squares fit of these 3D point estimates.
|
||||
Returns (x, y, z) estimate and covariance.
|
||||
"""
|
||||
if len(detections) < 2:
|
||||
return None, None
|
||||
|
||||
# Each radar's bearing + range defines a 3D point estimate
|
||||
estimates = []
|
||||
for det in detections:
|
||||
ri = radar_pos_list[det['radar_index']]
|
||||
az = det['azimuth']
|
||||
el = det['elevation']
|
||||
r = det['range']
|
||||
|
||||
# Convert bearing + range to 3D offset from radar
|
||||
dx = r * np.cos(el) * np.cos(az)
|
||||
dy = r * np.cos(el) * np.sin(az)
|
||||
dz = r * np.sin(el)
|
||||
|
||||
p_est = ri + np.array([dx, dy, dz])
|
||||
estimates.append(p_est)
|
||||
|
||||
# Least-squares estimate is the mean of the point estimates
|
||||
p = np.mean(estimates, axis=0)
|
||||
|
||||
# Estimate covariance from the spread of estimates
|
||||
if len(estimates) > 1:
|
||||
spread = np.std(np.array(estimates), axis=0)
|
||||
cov = np.diag(spread**2)
|
||||
else:
|
||||
sigma_pos = RADAR_ERROR_1SIGMA
|
||||
cov = np.eye(3) * sigma_pos**2
|
||||
|
||||
return p, cov
|
||||
|
||||
# =============================================================================
|
||||
# Trajectory Estimation (Linear Regression)
|
||||
# =============================================================================
|
||||
|
||||
def estimate_trajectory(observations, rng=None):
|
||||
"""
|
||||
Estimate trajectory (velocity) from a series of position observations.
|
||||
Each observation is (x, y, z, t).
|
||||
Returns (vx, vy, vz) and covariance.
|
||||
"""
|
||||
if len(observations) < 2:
|
||||
return None, None
|
||||
|
||||
times = np.array([o[3] for o in observations])
|
||||
xs = np.array([o[0] for o in observations])
|
||||
ys = np.array([o[1] for o in observations])
|
||||
zs = np.array([o[2] for o in observations])
|
||||
|
||||
n = len(observations)
|
||||
|
||||
# Linear regression for each axis: pos = pos0 + vel * t
|
||||
vel = np.zeros(3)
|
||||
pos0 = np.zeros(3)
|
||||
for axis in range(3):
|
||||
pos = np.array([o[axis] for o in observations])
|
||||
if n >= 2:
|
||||
# Least squares fit: pos = pos0 + vel*t
|
||||
A = np.vstack([np.ones(n), times]).T
|
||||
coeffs, _, _, _ = np.linalg.lstsq(A, pos, rcond=None)
|
||||
pos0[axis] = coeffs[0]
|
||||
vel[axis] = coeffs[1]
|
||||
|
||||
# Covariance estimate
|
||||
# For linear regression, Var(v) = sigma^2 / Sxx where Sxx = sum((t - t_mean)^2)
|
||||
t_mean = np.mean(times)
|
||||
Sxx = np.sum((times - t_mean)**2)
|
||||
sigma_vel = RADAR_ERROR_1SIGMA / np.sqrt(Sxx) if Sxx > 0 else RADAR_ERROR_1SIGMA
|
||||
cov = np.eye(3) * sigma_vel**2
|
||||
|
||||
return vel, cov
|
||||
|
||||
# =============================================================================
|
||||
# Drone Intercept Logic
|
||||
# =============================================================================
|
||||
|
||||
def compute_intercept_point(drone_pos, target_pos, target_vel, drone_speed):
|
||||
"""
|
||||
Compute the intercept point for a drone chasing a moving target in 3D.
|
||||
Solves the classic pursuer-target problem:
|
||||
|drone_pos + v_d*t - target_pos - v_s*t| = 0
|
||||
where v_d = drone_speed * unit_vector(direction to intercept)
|
||||
|
||||
This works even when the drone is slower than the target,
|
||||
as long as the drone is positioned close enough to the intercept path.
|
||||
|
||||
Returns: intercept point (3D) and estimated time
|
||||
"""
|
||||
drone_pos = np.array(drone_pos, dtype=float)
|
||||
target_pos = np.array(target_pos, dtype=float)
|
||||
target_vel = np.array(target_vel, dtype=float)
|
||||
|
||||
# Relative position and velocity
|
||||
rel_pos = target_pos - drone_pos # where target is relative to drone
|
||||
rel_vel = target_vel # relative velocity (drone direction is variable)
|
||||
|
||||
# We need to find a direction d (unit vector) and time t such that:
|
||||
# drone_speed * t * d = rel_pos + target_vel * t
|
||||
# => d = (rel_pos + target_vel * t) / (drone_speed * t)
|
||||
# And |d| = 1, so:
|
||||
# |rel_pos + target_vel * t| = drone_speed * t
|
||||
|
||||
# Let's expand:
|
||||
# |rel_pos + v_s*t|^2 = v_d^2 * t^2
|
||||
# |rel_pos|^2 + 2*t*(rel_pos·v_s) + t^2*(v_s·v_s) = v_d^2 * t^2
|
||||
# |rel_pos|^2 + 2*t*(rel_pos·v_s) + t^2*(|v_s|^2 - v_d^2) = 0
|
||||
|
||||
a = np.dot(target_vel, target_vel) - drone_speed**2
|
||||
b = 2 * np.dot(rel_pos, target_vel)
|
||||
c = np.dot(rel_pos, rel_pos)
|
||||
|
||||
# Solve quadratic: a*t^2 + b*t + c = 0
|
||||
discriminant = b*b - 4*a*c
|
||||
|
||||
if a == 0:
|
||||
# Linear equation: b*t + c = 0
|
||||
if abs(b) < 1e-10:
|
||||
# No solution
|
||||
t = np.linalg.norm(rel_pos) / drone_speed
|
||||
intercept = drone_pos + rel_pos
|
||||
return intercept, t
|
||||
t = -c / b
|
||||
if t <= 0:
|
||||
t = np.linalg.norm(rel_pos) / drone_speed
|
||||
intercept = drone_pos + rel_pos
|
||||
return intercept, t
|
||||
target_future = target_pos + target_vel * t
|
||||
intercept = target_future
|
||||
return intercept, t
|
||||
|
||||
if discriminant < 0:
|
||||
# No intercept possible with this geometry
|
||||
# Fallback: aim at current position
|
||||
t = np.linalg.norm(rel_pos) / drone_speed
|
||||
intercept = drone_pos + rel_pos
|
||||
return intercept, t
|
||||
|
||||
# Compute both roots, pick the smallest positive one
|
||||
sqrt_disc = np.sqrt(discriminant)
|
||||
t_candidates = [
|
||||
(-b + sqrt_disc) / (2 * a),
|
||||
(-b - sqrt_disc) / (2 * a)
|
||||
]
|
||||
|
||||
positive_roots = [t for t in t_candidates if t > 0]
|
||||
|
||||
if not positive_roots:
|
||||
# No valid positive intercept
|
||||
t = np.linalg.norm(rel_pos) / drone_speed
|
||||
intercept = drone_pos + rel_pos
|
||||
return intercept, t
|
||||
|
||||
# The smallest positive root is the earliest intercept
|
||||
t = min(positive_roots)
|
||||
|
||||
# Intercept point: where target will be at time t
|
||||
target_future = target_pos + target_vel * t
|
||||
intercept = target_future
|
||||
|
||||
return intercept, t
|
||||
|
||||
def drone_trajectory(drone_start, intercept_point, drone_speed):
|
||||
"""
|
||||
Generate a smooth 3D trajectory for the drone from start to intercept.
|
||||
Simple: linear interpolation with altitude management.
|
||||
"""
|
||||
drone_start = np.array(drone_start)
|
||||
intercept = np.array(intercept_point)
|
||||
|
||||
total_dist = np.linalg.norm(intercept - drone_start)
|
||||
duration = total_dist / drone_speed
|
||||
|
||||
# Generate waypoints
|
||||
num_steps = max(2, int(duration / DT))
|
||||
trajectory = []
|
||||
for i in range(num_steps + 1):
|
||||
frac = i / num_steps
|
||||
pos = drone_start + frac * (intercept - drone_start)
|
||||
t = i * DT
|
||||
trajectory.append((pos[0], pos[1], pos[2], t))
|
||||
|
||||
return trajectory
|
||||
|
||||
# =============================================================================
|
||||
# Main Simulation
|
||||
# =============================================================================
|
||||
|
||||
def run_simulation(seed=None, debug=False, preset=None):
|
||||
"""Run one full simulation and return results.
|
||||
preset: name of a projectile preset (or None for default).
|
||||
"""
|
||||
rng = np.random.default_rng(seed)
|
||||
preset_name = preset if preset is not None else DEFAULT_PRESET
|
||||
preset_data = get_preset(preset)
|
||||
|
||||
# Attack parameters from preset (with some randomization)
|
||||
start_dist = preset_data['start_dist'] * 1000.0 * rng.uniform(0.9, 1.1) # m, ±10%
|
||||
start_angle = rng.uniform(0, 2 * np.pi)
|
||||
target_angle = rng.uniform(0, 2 * np.pi)
|
||||
target_dist = rng.uniform(0, CITY_RADIUS * 0.8)
|
||||
|
||||
start_pos = np.array([
|
||||
start_dist * np.cos(start_angle),
|
||||
start_dist * np.sin(start_angle),
|
||||
0
|
||||
])
|
||||
|
||||
terminal_alt = preset_data['terminal_alt']
|
||||
target_pos = np.array([
|
||||
target_dist * np.cos(target_angle),
|
||||
target_dist * np.sin(target_angle),
|
||||
terminal_alt # terminal altitude
|
||||
])
|
||||
|
||||
cruise_alt = rng.uniform(preset_data['cruise_alt_min'], preset_data['cruise_alt_max'])
|
||||
speed = preset_data['speed']
|
||||
descent_start = preset_data['descent_start'] * 1000.0 # m
|
||||
|
||||
if debug:
|
||||
print(f" [DEBUG] preset={preset_data} start_dist={start_dist:.0f}m, cruise_alt={cruise_alt:.0f}m")
|
||||
print(f" [DEBUG] start_pos={start_pos}, target_pos={target_pos}")
|
||||
|
||||
# Generate trajectory
|
||||
trajectory = target_trajectory(start_pos, target_pos, cruise_alt, speed, descent_start, terminal_alt)
|
||||
if not trajectory:
|
||||
return None
|
||||
|
||||
# Setup radar and drones
|
||||
radar_pos_list = radar_positions(NUM_RADARS, CITY_RADIUS)
|
||||
drone_pos_list = drone_positions(NUM_DRONES, CITY_RADIUS, DRONE_DEFENSE_LINE)
|
||||
|
||||
# Simulation state
|
||||
time = 0.0
|
||||
detection_time = None
|
||||
launched_drones = []
|
||||
collision = False
|
||||
observations = []
|
||||
|
||||
# Pre-compute true target position at each time step
|
||||
true_positions = {t: np.array([x, y, z]) for x, y, z, t in trajectory}
|
||||
|
||||
def get_true_pos(t):
|
||||
# Find closest trajectory point
|
||||
times = sorted(true_positions.keys())
|
||||
# Binary search for closest
|
||||
idx = np.searchsorted(times, t)
|
||||
if idx >= len(times):
|
||||
return true_positions[times[-1]]
|
||||
if idx == 0:
|
||||
return true_positions[times[0]]
|
||||
# Choose closer of two adjacent
|
||||
t1, t2 = times[idx-1], times[idx]
|
||||
if abs(t - t1) < abs(t - t2):
|
||||
return true_positions[t1]
|
||||
return true_positions[t2]
|
||||
|
||||
# Main loop
|
||||
while time < SIM_DURATION and not collision:
|
||||
true_pos = get_true_pos(time)
|
||||
|
||||
# Check radar detection
|
||||
detections = detect_from_radars(true_pos, radar_pos_list, RADAR_RANGE, rng)
|
||||
active_radars = [d['radar_index'] for d in detections]
|
||||
|
||||
if active_radars:
|
||||
observations.append((true_pos[0], true_pos[1], true_pos[2], time))
|
||||
|
||||
# First detection triggers alert
|
||||
if detection_time is None and len(observations) >= 3:
|
||||
detection_time = time
|
||||
|
||||
if debug:
|
||||
print(f" [DEBUG] Detection at t={time:.1f}s, true_pos={true_pos}")
|
||||
|
||||
# Estimate 3D position (only once at detection)
|
||||
pos_est, cov = estimate_position_3d(detections, radar_pos_list, rng)
|
||||
if pos_est is not None:
|
||||
observations[-1] = (pos_est[0], pos_est[1], pos_est[2], time)
|
||||
if debug:
|
||||
print(f" [DEBUG] pos_est={pos_est}, true={true_pos}, err={np.linalg.norm(pos_est-true_pos):.1f}m")
|
||||
|
||||
# Estimate trajectory (only once at detection)
|
||||
vel_est, vel_cov = estimate_trajectory(observations, rng)
|
||||
if vel_est is not None:
|
||||
if debug:
|
||||
print(f" [DEBUG] vel_est={vel_est}, true_vel approx={speed} m/s")
|
||||
|
||||
# Compute time to impact for this trajectory
|
||||
# Target distance to city center
|
||||
target_horiz_dist = np.sqrt(true_pos[0]**2 + true_pos[1]**2)
|
||||
speed_est = np.linalg.norm(vel_est)
|
||||
t_to_impact = target_horiz_dist / speed_est if speed_est > 0 else float('inf')
|
||||
|
||||
if debug:
|
||||
print(f" [DEBUG] t_to_impact={t_to_impact:.1f}s, target_horiz_dist={target_horiz_dist:.0f}m")
|
||||
|
||||
# Launch drones that can intercept before impact
|
||||
drone_intercepts = []
|
||||
for i in range(NUM_DRONES):
|
||||
drone_pos = drone_pos_list[i]['pos']
|
||||
intercept, t_int = compute_intercept_point(
|
||||
drone_pos, true_pos, vel_est, DRONE_SPEED
|
||||
)
|
||||
# Only viable if intercept happens before impact
|
||||
if t_int < t_to_impact:
|
||||
drone_intercepts.append((t_int, i, intercept))
|
||||
|
||||
drone_intercepts.sort(key=lambda x: x[0])
|
||||
launched_drones = [d[1] for d in drone_intercepts[:3]]
|
||||
|
||||
if debug:
|
||||
print(f" [DEBUG] viable drones: {len(drone_intercepts)}, launched: {launched_drones}")
|
||||
for d in drone_intercepts[:3]:
|
||||
print(f" [DEBUG] drone {d[1]} t_int={d[0]:.1f}s intercept=({d[2][0]:.0f},{d[2][1]:.0f},{d[2][2]:.0f})")
|
||||
|
||||
# Compute intercept points
|
||||
for d in drone_intercepts[:3]:
|
||||
i = d[1]
|
||||
drone_pos_list[i]['intercept'] = d[2]
|
||||
drone_pos_list[i]['state'] = 'flying'
|
||||
|
||||
# Update drones
|
||||
if launched_drones:
|
||||
# Re-target drones periodically (proportional navigation)
|
||||
if detection_time is not None and (time - detection_time) % 5.0 < DT:
|
||||
vel_est_new, _ = estimate_trajectory(observations, rng)
|
||||
if vel_est_new is not None:
|
||||
# Recompute intercept points for launched drones
|
||||
for i in launched_drones:
|
||||
drone = drone_pos_list[i]
|
||||
if drone['state'] in ('flying', 'ir_tracking'):
|
||||
intercept, t_int = compute_intercept_point(
|
||||
drone['pos'], true_pos, vel_est_new, DRONE_SPEED
|
||||
)
|
||||
drone['intercept'] = intercept
|
||||
|
||||
for i in launched_drones:
|
||||
drone = drone_pos_list[i]
|
||||
if drone['state'] in ('flying', 'ir_tracking'):
|
||||
intercept = drone.get('intercept')
|
||||
if intercept is not None:
|
||||
delta = intercept - drone['pos']
|
||||
dist = np.linalg.norm(delta)
|
||||
if dist > 0.01:
|
||||
step = DRONE_SPEED * DT
|
||||
drone['pos'] = drone['pos'] + (delta / dist) * step
|
||||
else:
|
||||
drone['pos'] = intercept
|
||||
|
||||
# Check IR range
|
||||
target_dist_local = np.linalg.norm(true_pos - drone['pos'])
|
||||
if drone['state'] == 'flying' and target_dist_local < IR_RANGE:
|
||||
drone['state'] = 'ir_tracking'
|
||||
if debug:
|
||||
print(f" [DEBUG] t={time:.1f}s drone {i} -> IR tracking, dist={target_dist_local:.0f}m")
|
||||
|
||||
# Check collision
|
||||
if target_dist_local < 20.0: # 20 m kill radius
|
||||
collision = True
|
||||
if debug:
|
||||
print(f" [DEBUG] COLLISION at t={time:.1f}s drone {i}")
|
||||
break
|
||||
elif debug and time % 10 < DT:
|
||||
print(f" [DEBUG] t={time:.1f}s drone {i} pos=({drone['pos'][0]:.0f},{drone['pos'][1]:.0f},{drone['pos'][2]:.0f}), target=({true_pos[0]:.0f},{true_pos[1]:.0f},{true_pos[2]:.0f}), dist={target_dist_local:.0f}m")
|
||||
|
||||
# Check if target reached city
|
||||
horiz_dist = np.sqrt(true_pos[0]**2 + true_pos[1]**2)
|
||||
if horiz_dist < CITY_RADIUS and not collision:
|
||||
break
|
||||
|
||||
time += DT
|
||||
|
||||
# Calculate results
|
||||
result = {
|
||||
'preset': preset_name,
|
||||
'start_dist': start_dist,
|
||||
'cruise_alt': cruise_alt,
|
||||
'speed': speed,
|
||||
'terminal_alt': terminal_alt,
|
||||
'detection_time': detection_time,
|
||||
'collision': collision,
|
||||
'time_of_collision': time if collision else None,
|
||||
'time_of_impact': time if not collision else None,
|
||||
'num_observations': len(observations),
|
||||
'launched_drones': launched_drones,
|
||||
'seed': seed
|
||||
}
|
||||
|
||||
return result
|
||||
|
||||
# =============================================================================
|
||||
# Monte Carlo Validation
|
||||
# =============================================================================
|
||||
|
||||
def monte_carlo(n_trials=100, seed_base=42, preset=None):
|
||||
"""Run Monte Carlo simulation and return statistics.
|
||||
preset: name of a projectile preset (or None for default).
|
||||
"""
|
||||
results = []
|
||||
for i in range(n_trials):
|
||||
result = run_simulation(seed=seed_base + i, preset=preset)
|
||||
if result is not None:
|
||||
results.append(result)
|
||||
|
||||
if not results:
|
||||
return None
|
||||
|
||||
# Statistics
|
||||
collisions = sum(1 for r in results if r['collision'])
|
||||
detection_times = [r['detection_time'] for r in results if r['detection_time'] is not None]
|
||||
collision_times = [r['time_of_collision'] for r in results if r['collision']]
|
||||
impact_times = [r['time_of_impact'] for r in results if not r['collision']]
|
||||
num_obs = [r['num_observations'] for r in results]
|
||||
cruise_alts = [r['cruise_alt'] for r in results]
|
||||
preset_name = results[0]['preset']
|
||||
|
||||
stats = {
|
||||
'preset': preset_name,
|
||||
'total_trials': len(results),
|
||||
'collisions': collisions,
|
||||
'interception_rate': collisions / len(results) * 100,
|
||||
'avg_detection_time': np.mean(detection_times) if detection_times else None,
|
||||
'avg_collision_time': np.mean(collision_times) if collision_times else None,
|
||||
'avg_impact_time': np.mean(impact_times) if impact_times else None,
|
||||
'avg_observations': np.mean(num_obs),
|
||||
'min_observations': np.min(num_obs),
|
||||
'max_observations': np.max(num_obs),
|
||||
'avg_cruise_alt': np.mean(cruise_alts),
|
||||
'min_cruise_alt': np.min(cruise_alts),
|
||||
'max_cruise_alt': np.max(cruise_alts),
|
||||
}
|
||||
|
||||
return stats, results
|
||||
|
||||
# =============================================================================
|
||||
# Main
|
||||
# =============================================================================
|
||||
|
||||
if __name__ == '__main__':
|
||||
print("=" * 80)
|
||||
print("Generalized Projectile Defense Simulation")
|
||||
print("=" * 80)
|
||||
print()
|
||||
|
||||
# Run single simulation for demo (default preset)
|
||||
print(f"Running single simulation (preset={DEFAULT_PRESET})...")
|
||||
result = run_simulation(seed=42)
|
||||
if result:
|
||||
print(f" Preset: {result['preset']}")
|
||||
print(f" Start distance: {result['start_dist']/1000:.1f} km")
|
||||
print(f" Speed: {result['speed']:.0f} m/s")
|
||||
print(f" Cruise altitude: {result['cruise_alt']:.0f} m")
|
||||
print(f" Detection time: {result['detection_time']:.1f} s" if result['detection_time'] else " Detection time: N/A")
|
||||
print(f" Collision: {'YES' if result['collision'] else 'NO'}")
|
||||
if result['collision']:
|
||||
print(f" Collision time: {result['time_of_collision']:.1f} s")
|
||||
else:
|
||||
print(f" Impact time: {result['time_of_impact']:.1f} s")
|
||||
print(f" Observations: {result['num_observations']}")
|
||||
print()
|
||||
|
||||
# Run Monte Carlo for each preset
|
||||
print("Running Monte Carlo (100 trials per preset)...")
|
||||
all_stats = {}
|
||||
for preset_name in PRESETS.keys():
|
||||
stats, results = monte_carlo(n_trials=100, seed_base=100, preset=preset_name)
|
||||
if stats:
|
||||
all_stats[preset_name] = stats
|
||||
print(f" [{preset_name}]")
|
||||
print(f" Interception rate: {stats['interception_rate']:.1f}%")
|
||||
if stats['avg_detection_time']:
|
||||
print(f" Avg detection time: {stats['avg_detection_time']:.1f} s")
|
||||
if stats['avg_collision_time']:
|
||||
print(f" Avg collision time: {stats['avg_collision_time']:.1f} s")
|
||||
if stats['avg_impact_time']:
|
||||
print(f" Avg impact time: {stats['avg_impact_time']:.1f} s")
|
||||
print(f" Avg cruise altitude: {stats['avg_cruise_alt']:.0f} m")
|
||||
print()
|
||||
|
||||
# Save stats to JSON for the HTML viewer (per preset)
|
||||
import json
|
||||
json_stats = {
|
||||
preset_name: {
|
||||
'total_trials': stats['total_trials'],
|
||||
'collisions': stats['collisions'],
|
||||
'interception_rate': stats['interception_rate'],
|
||||
'avg_detection_time': stats['avg_detection_time'],
|
||||
'avg_collision_time': stats['avg_collision_time'],
|
||||
'avg_impact_time': stats['avg_impact_time'],
|
||||
'avg_observations': stats['avg_observations'],
|
||||
}
|
||||
for preset_name, stats in all_stats.items()
|
||||
}
|
||||
with open('montecarlo_stats.json', 'w') as f:
|
||||
json.dump(json_stats, f, indent=2)
|
||||
print()
|
||||
print("Stats saved to montecarlo_stats.json")
|
||||
|
||||
print()
|
||||
print("=" * 80)
|
||||
print("Simulation complete")
|
||||
print("=" * 80)
|
||||
@@ -1,21 +0,0 @@
|
||||
MIT License
|
||||
|
||||
Copyright (c) 2020-PRESENT Anthony Fu
|
||||
|
||||
Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
of this software and associated documentation files (the "Software"), to deal
|
||||
in the Software without restriction, including without limitation the rights
|
||||
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
copies of the Software, and to permit persons to whom the Software is
|
||||
furnished to do so, subject to the following conditions:
|
||||
|
||||
The above copyright notice and this permission notice shall be included in all
|
||||
copies or substantial portions of the Software.
|
||||
|
||||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||
SOFTWARE.
|
||||
@@ -1,50 +0,0 @@
|
||||
export interface SimulationEntity {
|
||||
id: string
|
||||
type: string
|
||||
active: boolean
|
||||
position_ned_m: number[]
|
||||
velocity_ned_m_s: number[]
|
||||
acceleration_ned_m_s2: number[]
|
||||
orientation_quaternion_wxyz?: number[]
|
||||
angular_velocity_body_rad_s?: number[]
|
||||
available?: boolean
|
||||
activated_for?: string | null
|
||||
position_error_ned_m?: number[] | null
|
||||
velocity_error_ned_m_s?: number[] | null
|
||||
trajectory_ned_m?: number[][]
|
||||
intended_trajectory_ned_m?: number[][]
|
||||
intention_kind?: string
|
||||
standby_altitude_m?: number | null
|
||||
flight_phase?: string
|
||||
}
|
||||
|
||||
export interface SimulationEvent {
|
||||
type: string
|
||||
time_s: number
|
||||
threat_id?: string
|
||||
defender_ids?: string[]
|
||||
projectile_id?: string
|
||||
target_id?: string
|
||||
position_ned_m?: number[]
|
||||
projected_city_crossing_ned_m?: number[]
|
||||
miss_distance_m?: number
|
||||
[key: string]: unknown
|
||||
}
|
||||
|
||||
export interface SimulationState {
|
||||
schema_version: number
|
||||
run_id: string
|
||||
sequence: number
|
||||
simulation_time_s: number
|
||||
kind: 'simulation_state'
|
||||
entities: SimulationEntity[]
|
||||
events: SimulationEvent[]
|
||||
environment: Record<string, unknown>
|
||||
devices?: Record<string, unknown>
|
||||
scenario?: Record<string, unknown>
|
||||
observations?: Record<string, unknown>
|
||||
radar_detection_points?: Record<string, number[][]>
|
||||
radar_tracks?: Record<string, { position_ned_m: number[], velocity_ned_m_s: number[], timestamp_s: number, age_s: number, velocity_observable: boolean }>
|
||||
inferred_trajectories?: Record<string, number[][]>
|
||||
statistics: Record<string, unknown>
|
||||
}
|
||||
@@ -1,3 +0,0 @@
|
||||
"""Backend-authoritative defense simulation."""
|
||||
|
||||
__version__ = "0.1.0"
|
||||
@@ -1,325 +0,0 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
import numpy as np
|
||||
from ..physics.validation import finite_numeric_fields
|
||||
from ..physics.atmosphere import isa
|
||||
from ..physics.integrator import NUMBA_AVAILABLE, rk4_fast
|
||||
if NUMBA_AVAILABLE:
|
||||
from ..physics.integrator import drone_rk4_kernel
|
||||
else: # pragma: no cover - numba optional
|
||||
def drone_rk4_kernel(state, wind, held, drag_k, rho):
|
||||
rel = state[3:] - wind
|
||||
speed = float(np.linalg.norm(rel))
|
||||
return np.concatenate((state[3:], held + np.array([0.0, 0.0, 9.80665])
|
||||
- 0.5 * rho * drag_k * speed * rel))
|
||||
|
||||
|
||||
@dataclass
|
||||
class DroneConfig:
|
||||
mass_kg: float = 2.0
|
||||
max_acceleration_m_s2: float = 18.0
|
||||
max_speed_m_s: float = 25.0
|
||||
max_turn_rate_rad_s: float = 1.2
|
||||
max_climb_rate_m_s: float = 5.0
|
||||
max_descent_rate_m_s: float = 4.0
|
||||
acceleration_time_constant_s: float = 0.25
|
||||
drag_coefficient_area_m2: float = 0.08
|
||||
collision_radius_m: float = 0.5
|
||||
actuator_sigma_m_s2: float = 0.0
|
||||
# IR seeker model (used once the drone is tasked against a threat).
|
||||
ir_range_m: float = 8_000.0
|
||||
ir_bearing_sigma_rad: float = .005
|
||||
ir_fov_deg: float = 120.0
|
||||
ir_position_sigma_m: float = 5.0
|
||||
ir_velocity_sigma_m_s: float = 1.0
|
||||
ir_update_rate_hz: float = 10.0
|
||||
ir_processing_latency_s: float = 0.01
|
||||
ir_lock_time_s: float = 1.0
|
||||
|
||||
@property
|
||||
def ir_acquisition_time_s(self) -> float:
|
||||
# Lock-on margin: seeker needs the configured lock time plus one
|
||||
# processing/sample period before the estimate is usable.
|
||||
return self.ir_lock_time_s + 1.0 / self.ir_update_rate_hz + self.ir_processing_latency_s
|
||||
|
||||
@property
|
||||
def ir_acquisition_range_m(self) -> float:
|
||||
# Minimum distance at which the seeker can hold a lock long enough to
|
||||
# produce a usable estimate. This is the standoff the interceptor
|
||||
# geometry is built around.
|
||||
return self.ir_acquisition_time_s * self.max_speed_m_s
|
||||
|
||||
def __post_init__(self):
|
||||
finite_numeric_fields(self)
|
||||
if self.acceleration_time_constant_s <= 0 or self.max_turn_rate_rad_s < 0 or self.max_climb_rate_m_s <= 0 or self.max_descent_rate_m_s <= 0:
|
||||
raise ValueError("invalid drone response/flight limits")
|
||||
if self.drag_coefficient_area_m2 < 0 or self.collision_radius_m < 0 or self.actuator_sigma_m_s2 < 0:
|
||||
raise ValueError("invalid drone drag/contact/noise parameters")
|
||||
if self.mass_kg <= 0 or self.max_acceleration_m_s2 <= 0 or self.max_speed_m_s <= 0:
|
||||
raise ValueError("invalid drone physical parameters")
|
||||
if self.ir_range_m <= 0 or self.ir_position_sigma_m < 0 or self.ir_velocity_sigma_m_s < 0:
|
||||
raise ValueError("invalid IR seeker parameters")
|
||||
if self.ir_bearing_sigma_rad < 0 or not 0 < self.ir_fov_deg <= 360:
|
||||
raise ValueError("invalid IR bearing noise/FOV")
|
||||
if self.ir_update_rate_hz <= 0 or self.ir_processing_latency_s < 0 or self.ir_lock_time_s < 0:
|
||||
raise ValueError("invalid IR seeker timing")
|
||||
|
||||
|
||||
@dataclass
|
||||
class Drone:
|
||||
entity_id: str
|
||||
entity_type: str
|
||||
position_ned_m: np.ndarray
|
||||
velocity_ned_m_s: np.ndarray = field(default_factory=lambda: np.zeros(3))
|
||||
config: DroneConfig = field(default_factory=DroneConfig)
|
||||
acceleration_ned_m_s2: np.ndarray = field(default_factory=lambda: np.zeros(3))
|
||||
estimated_position_ned_m: np.ndarray | None = None
|
||||
estimated_velocity_ned_m_s: np.ndarray | None = None
|
||||
actuator_acceleration_ned_m_s2: np.ndarray = field(default_factory=lambda: np.zeros(3))
|
||||
requested_thrust_ned_m_s2: np.ndarray = field(default_factory=lambda: np.zeros(3))
|
||||
orientation_quaternion_wxyz: np.ndarray = field(default_factory=lambda: np.array([1., 0., 0., 0.]))
|
||||
angular_velocity_body_rad_s: np.ndarray = field(default_factory=lambda: np.zeros(3))
|
||||
target_position_ned_m: np.ndarray | None = None
|
||||
active: bool = True
|
||||
available: bool = True
|
||||
activated_for: str | None = None
|
||||
# IR seeker state (threat-side): last estimate, lock age, acquisition flag.
|
||||
ir_threat_estimate: np.ndarray | None = None # [position(3), velocity(3)]
|
||||
ir_lock_age_s: float = 0.0
|
||||
ir_acquired: bool = False
|
||||
ir_next_sample_s: float = 0.0
|
||||
ir_rng: np.random.Generator | None = None
|
||||
flight_dynamics: object | None = None
|
||||
flight_reference: dict | None = None
|
||||
standby_altitude_m: float | None = None
|
||||
ir_pending: list = field(default_factory=list)
|
||||
ir_observation: dict | None = None
|
||||
|
||||
def estimate_threat(self, target, time_s: float, rng: np.random.Generator,
|
||||
event_sink: list | None = None, radar_prior=None) -> bool:
|
||||
"""Sample a rear-facing bearing-only IR camera, then deliver after delay.
|
||||
|
||||
Visibility uses truth only in measurement generation. Range is never
|
||||
measured by IR: the EKF retains the radar prior's range uncertainty.
|
||||
Ideal attitude/gyro knowledge and no terrain occlusion are assumptions.
|
||||
"""
|
||||
if self.activated_for is None or self.ir_rng is None or target is None or not target.active:
|
||||
return False
|
||||
if time_s + 1e-9 >= self.ir_next_sample_s:
|
||||
period = 1/self.config.ir_update_rate_hz
|
||||
while self.ir_next_sample_s <= time_s+1e-9:
|
||||
self.ir_next_sample_s += period
|
||||
relative = target.position_ned_m-self.position_ned_m
|
||||
distance = float(np.linalg.norm(relative))
|
||||
from scipy.spatial.transform import Rotation
|
||||
rotation = Rotation.from_quat(self.orientation_quaternion_wxyz, scalar_first=True).as_matrix()
|
||||
camera_axis = -rotation[:, 0] # rear-looking mount on body x axis
|
||||
visible = (0 < distance <= self.config.ir_range_m and
|
||||
float(relative @ camera_axis/distance) >= np.cos(np.deg2rad(self.config.ir_fov_deg/2)))
|
||||
if not visible:
|
||||
if self.ir_acquired and event_sink is not None:
|
||||
event_sink.append({"type": "ir_track_lost", "time_s": time_s, "drone_id": self.entity_id})
|
||||
self.ir_lock_age_s = 0.
|
||||
self.ir_acquired = False
|
||||
self.ir_pending.clear()
|
||||
self.ir_threat_estimate = None
|
||||
self.ir_observation = None
|
||||
else:
|
||||
self.ir_lock_age_s += period
|
||||
if self.ir_lock_age_s >= self.config.ir_lock_time_s:
|
||||
# Elevation is mathematical +down to match Stone Soup's
|
||||
# NED Cartesian model; physical elevation is its negative.
|
||||
angles = np.array([np.arctan2(relative[2], np.linalg.norm(relative[:2])),
|
||||
np.arctan2(relative[1], relative[0])])
|
||||
angles += self.ir_rng.normal(0, self.config.ir_bearing_sigma_rad, 2)
|
||||
sensor_position = self.estimated_position_ned_m if self.estimated_position_ned_m is not None else self.position_ned_m
|
||||
self.ir_pending.append({"sample_time_s": time_s,
|
||||
"ready_time_s": time_s+self.config.ir_processing_latency_s,
|
||||
"angles": angles, "sensor_position": sensor_position.copy()})
|
||||
due = [m for m in self.ir_pending if m["ready_time_s"] <= time_s+1e-9]
|
||||
self.ir_pending = [m for m in self.ir_pending if m["ready_time_s"] > time_s+1e-9]
|
||||
if not due or radar_prior is None:
|
||||
return False
|
||||
from ..sensors.tracking import refine_bearing
|
||||
measurement = due[-1]
|
||||
estimate = refine_bearing(radar_prior, measurement["sensor_position"], measurement["angles"],
|
||||
measurement["sample_time_s"], self.config.ir_bearing_sigma_rad)
|
||||
if estimate is None:
|
||||
return False
|
||||
acquired = not self.ir_acquired
|
||||
self.ir_acquired = True
|
||||
self.ir_threat_estimate = np.r_[estimate.position, estimate.velocity]
|
||||
self.ir_observation = {"sample_time_s": estimate.timestamp_s, "receive_time_s": time_s,
|
||||
"measurement_type": "bearing_only", "elevation_down_rad": float(measurement["angles"][0]),
|
||||
"bearing_rad": float(measurement["angles"][1]), "covariance": estimate.covariance.tolist()}
|
||||
if self.flight_reference is not None:
|
||||
# Bearing-informed corrections retain radar range uncertainty.
|
||||
# Keep along-path progress; smooth the lateral reference instead
|
||||
# of overwriting the drone's actual position or velocity.
|
||||
self.flight_reference["target_altitude_m"] = max(0., -estimate.position[2])
|
||||
direction = self.flight_reference["direction"][:2]
|
||||
relative_estimate = estimate.position[:2]-measurement["sensor_position"][:2]
|
||||
lateral = relative_estimate-direction*float(relative_estimate @ direction)
|
||||
bound = self.config.max_speed_m_s/self.config.ir_update_rate_hz
|
||||
lateral *= min(1., bound/max(float(np.linalg.norm(lateral)), 1e-9))
|
||||
self.flight_reference["anchor"][:2] += lateral
|
||||
if event_sink is not None:
|
||||
event_sink.append({"type": "ir_acquisition" if acquired else "ir_track_update",
|
||||
"time_s": time_s, "drone_id": self.entity_id, "threat_id": self.activated_for,
|
||||
**self.ir_observation})
|
||||
return True
|
||||
|
||||
def control(self, dt: float, rng: np.random.Generator) -> None:
|
||||
if self.target_position_ned_m is None:
|
||||
desired_net = np.zeros(3)
|
||||
else:
|
||||
est = self.estimated_position_ned_m if self.estimated_position_ned_m is not None else self.position_ned_m
|
||||
vel = self.estimated_velocity_ned_m_s if self.estimated_velocity_ned_m_s is not None else self.velocity_ned_m_s
|
||||
desired_net = 1.4 * (self.target_position_ned_m - est) - 1.1 * vel
|
||||
# Hover requires upward thrust to cancel NED gravity even at zero commanded net acceleration.
|
||||
desired = desired_net + np.array([0.0, 0.0, -9.80665])
|
||||
# Preserve vertical support before allocating remaining thrust to
|
||||
# planar motion. A distant horizontal reference must not silently
|
||||
# remove hover thrust and make a grounded drone slide indefinitely.
|
||||
limit = self.config.max_acceleration_m_s2
|
||||
desired[2] = np.clip(desired[2], -limit, limit)
|
||||
horizontal = float(np.linalg.norm(desired[:2]))
|
||||
budget = np.sqrt(max(0., limit*limit-desired[2]*desired[2]))
|
||||
if horizontal > budget:
|
||||
desired[:2] *= budget/max(horizontal, 1e-9)
|
||||
self.requested_thrust_ned_m_s2 = desired + rng.normal(0, self.config.actuator_sigma_m_s2, 3)
|
||||
n = np.linalg.norm(self.requested_thrust_ned_m_s2)
|
||||
if n > self.config.max_acceleration_m_s2:
|
||||
self.requested_thrust_ned_m_s2 *= self.config.max_acceleration_m_s2 / n
|
||||
|
||||
def update_flight_reference(self, dt):
|
||||
"""A bounded gameplay ascent reference, never a position overwrite."""
|
||||
ref = self.flight_reference
|
||||
if ref is None:
|
||||
if self.standby_altitude_m is not None and self.target_position_ned_m is not None:
|
||||
desired_down = -self.standby_altitude_m
|
||||
self.target_position_ned_m[2] += float(np.clip(desired_down-self.target_position_ned_m[2],
|
||||
-self.config.max_climb_rate_m_s*dt, self.config.max_descent_rate_m_s*dt))
|
||||
return
|
||||
ref["elapsed_s"] += dt
|
||||
previous = ref["previous_position"]
|
||||
ref["travel_m"] += float(np.linalg.norm(self.position_ned_m[:2]-previous[:2]))
|
||||
ref["previous_position"] = self.position_ned_m.copy()
|
||||
initial_h = ref["initial_altitude_m"]
|
||||
from ..guidance.ascent import altitude_reference
|
||||
height = altitude_reference(ref["travel_m"], initial_h, ref["target_altitude_m"],
|
||||
curve=ref.get("ascent_curve", "quintic"), length_m=ref.get("ascent_distance_m", 1000.),
|
||||
k_per_m=ref["ascent_k_per_m"])
|
||||
# Limit the vertical reference rate; the actual vehicle has its own
|
||||
# actuator/force dynamics and may fail to track this requested curve.
|
||||
old_height = -self.target_position_ned_m[2]
|
||||
height = old_height + float(np.clip(height-old_height,
|
||||
-self.config.max_descent_rate_m_s*dt, self.config.max_climb_rate_m_s*dt))
|
||||
target = ref["anchor"].copy()
|
||||
target[:2] += ref["direction"][:2]*ref["speed_m_s"]*ref["elapsed_s"]
|
||||
target[2] = -height
|
||||
# Smooth radar replans and lateral joins at the configured speed.
|
||||
delta = target[:2]-self.target_position_ned_m[:2]
|
||||
delta *= min(1., self.config.max_speed_m_s*dt/max(float(np.linalg.norm(delta)), 1e-9))
|
||||
self.target_position_ned_m[:2] += delta
|
||||
self.target_position_ned_m[2] = target[2]
|
||||
|
||||
def intended_trajectory(self, horizon_s=60., interval_s=1.):
|
||||
"""Preview the delivered command reference, not predicted physical flight.
|
||||
|
||||
Uses onboard estimates, rate limits and the selected altitude curve.
|
||||
Estimated horizontal travel advances the curve; actual flight may lag.
|
||||
No controller, actuator, wind or collision prediction is implied.
|
||||
"""
|
||||
if not self.active or self.target_position_ned_m is None:
|
||||
return []
|
||||
from ..guidance.ascent import altitude_reference
|
||||
origin = (self.estimated_position_ned_m if self.estimated_position_ned_m is not None
|
||||
else self.position_ned_m).copy()
|
||||
point = origin.copy()
|
||||
rows = [point.tolist()]
|
||||
ref = self.flight_reference
|
||||
travel = ref["travel_m"] if ref else 0.
|
||||
for elapsed in np.arange(interval_s, horizon_s+1e-9, interval_s):
|
||||
target = self.target_position_ned_m.copy()
|
||||
if ref is not None:
|
||||
target[:2] = ref["anchor"][:2] + ref["direction"][:2]*ref["speed_m_s"]*(ref["elapsed_s"]+elapsed)
|
||||
target[2] = -altitude_reference(travel, ref["initial_altitude_m"], ref["target_altitude_m"],
|
||||
curve=ref.get("ascent_curve", "quintic"), length_m=ref.get("ascent_distance_m", 1000.),
|
||||
k_per_m=ref["ascent_k_per_m"])
|
||||
elif self.standby_altitude_m is not None:
|
||||
target[2] = -self.standby_altitude_m
|
||||
delta = target[:2]-point[:2]
|
||||
delta *= min(1., self.config.max_speed_m_s*interval_s/max(float(np.linalg.norm(delta)), 1e-9))
|
||||
travel += float(np.linalg.norm(delta))
|
||||
point[:2] += delta
|
||||
point[2] += float(np.clip(target[2]-point[2], -self.config.max_climb_rate_m_s*interval_s,
|
||||
self.config.max_descent_rate_m_s*interval_s))
|
||||
rows.append(point.tolist())
|
||||
return rows
|
||||
|
||||
def step(self, dt: float, wind: np.ndarray, rng: np.random.Generator,
|
||||
control_due: bool = True, control_dt: float | None = None) -> None:
|
||||
if not self.active:
|
||||
return
|
||||
self.update_flight_reference(dt)
|
||||
if self.flight_dynamics is not None:
|
||||
if self.activated_for is None and self.target_position_ned_m is None and self.position_ned_m[2] >= 0:
|
||||
return
|
||||
self.flight_dynamics.step(self, dt, wind, control_due)
|
||||
return
|
||||
# A defender that has not been tasked simply holds station on the
|
||||
# loiter ring: it still responds to its actuator time constant (so a
|
||||
# gust or a fresh command produces a realistic ramp) but it does not
|
||||
# run the full guidance/controller/RK4 pipeline. This keeps idle
|
||||
# units cheap -- the engine should not pay flight-control cost for
|
||||
# drones that are not fighting -- and matches the doctrine that only
|
||||
# the tasked interceptors fly.
|
||||
if self.activated_for is None and self.target_position_ned_m is None:
|
||||
if self.position_ned_m[2] >= 0 and np.linalg.norm(self.velocity_ned_m_s) < 1e-9:
|
||||
self.position_ned_m[2] = 0.0
|
||||
self.acceleration_ned_m_s2 = np.zeros(3)
|
||||
return # supported by the ground, not an unpowered hover
|
||||
if not hasattr(self, "_hold_position_ned_m"):
|
||||
self._hold_position_ned_m = self.position_ned_m.copy()
|
||||
self.target_position_ned_m = self._hold_position_ned_m
|
||||
# Tasked flight: controller output sampled at its device cadence; the
|
||||
# held actuator state and vehicle dynamics continue at physics Hz.
|
||||
if control_due:
|
||||
self.control(control_dt if control_dt is not None else dt, rng)
|
||||
# Motor surrogate evolves at physics cadence, independently of the
|
||||
# command sample rate. Exact first-order decay avoids Euler clipping.
|
||||
old_actuator = self.actuator_acceleration_ned_m_s2.copy()
|
||||
tau = self.config.acceleration_time_constant_s
|
||||
decay = np.exp(-dt / tau)
|
||||
self.actuator_acceleration_ned_m_s2 = self.requested_thrust_ned_m_s2 + (old_actuator-self.requested_thrust_ned_m_s2)*decay
|
||||
wind = np.asarray(wind, dtype=np.float64)
|
||||
rho = isa(-self.position_ned_m[2]).density_kg_m3
|
||||
drag_k = self.config.drag_coefficient_area_m2 / self.config.mass_kg
|
||||
held = self.requested_thrust_ned_m_s2 + (old_actuator-self.requested_thrust_ned_m_s2)*tau*(1-decay)/dt
|
||||
old_velocity = self.velocity_ned_m_s.copy()
|
||||
y = rk4_fast(np.concatenate((self.position_ned_m, self.velocity_ned_m_s)), dt,
|
||||
drone_rk4_kernel, wind, held, drag_k, rho)
|
||||
self.velocity_ned_m_s = y[3:]
|
||||
speed = float(np.linalg.norm(self.velocity_ned_m_s))
|
||||
if speed > self.config.max_speed_m_s:
|
||||
self.velocity_ned_m_s *= self.config.max_speed_m_s / speed
|
||||
self.velocity_ned_m_s[2] = float(np.clip(self.velocity_ned_m_s[2], -self.config.max_climb_rate_m_s, self.config.max_descent_rate_m_s))
|
||||
# The reduced model treats speed bounds as kinematic constraints.
|
||||
# Integrate realized bounded velocity rather than unclipped endpoint
|
||||
# velocity; RotorPy is the force/attitude model for higher fidelity.
|
||||
self.position_ned_m += .5*dt*(old_velocity+self.velocity_ned_m_s)
|
||||
self.acceleration_ned_m_s2 = (self.velocity_ned_m_s-old_velocity)/dt
|
||||
if self.position_ned_m[2] > 0:
|
||||
self.position_ned_m[2] = 0.0
|
||||
self.velocity_ned_m_s[2] = min(0.0, self.velocity_ned_m_s[2])
|
||||
horizontal_speed = np.linalg.norm(self.velocity_ned_m_s[:2])
|
||||
if horizontal_speed > 1e-6:
|
||||
yaw = float(np.arctan2(self.velocity_ned_m_s[1], self.velocity_ned_m_s[0]))
|
||||
current_yaw = float(2 * np.arctan2(self.orientation_quaternion_wxyz[3], self.orientation_quaternion_wxyz[0]))
|
||||
error = (yaw - current_yaw + np.pi) % (2 * np.pi) - np.pi
|
||||
yaw_step = float(np.clip(error, -self.config.max_turn_rate_rad_s * dt, self.config.max_turn_rate_rad_s * dt))
|
||||
new_yaw = current_yaw + yaw_step
|
||||
self.angular_velocity_body_rad_s = np.array([0.0, 0.0, yaw_step / dt])
|
||||
self.orientation_quaternion_wxyz = np.array([np.cos(new_yaw/2), 0.0, 0.0, np.sin(new_yaw/2)])
|
||||
@@ -1,30 +0,0 @@
|
||||
"""Distance-based altitude references; these never prescribe actual motion."""
|
||||
import numpy as np
|
||||
|
||||
|
||||
def altitude_reference(distance_m, initial_altitude_m, target_altitude_m, *,
|
||||
curve="quintic", length_m=1000., k_per_m=.002):
|
||||
values = (distance_m, initial_altitude_m, target_altitude_m, length_m, k_per_m)
|
||||
if not all(np.isfinite(value) for value in values) or length_m <= 0 or k_per_m <= 0:
|
||||
raise ValueError("ascent reference parameters must be finite with positive L and k")
|
||||
x = max(0., distance_m)
|
||||
u = float(np.clip(x/length_m, 0., 1.))
|
||||
if curve == "quintic":
|
||||
fraction = u*u*u*(10+u*(-15+6*u))
|
||||
elif curve == "cubic":
|
||||
fraction = u*u*(3-2*u)
|
||||
elif curve == "exponential":
|
||||
fraction = -np.expm1(-k_per_m*x)
|
||||
else:
|
||||
raise ValueError("ascent curve must be quintic, cubic or exponential")
|
||||
return float(initial_altitude_m+(target_altitude_m-initial_altitude_m)*fraction)
|
||||
|
||||
|
||||
def automatic_ascent_length(rise_m, reference_speed_m_s, max_climb_rate_m_s, curve="quintic"):
|
||||
"""Heuristic L with 20% peak climb-rate margin at constant reference speed.
|
||||
|
||||
Quintic peak dh/dx is 1.875*rise/L; cubic peak is 1.5*rise/L.
|
||||
This does not guarantee acceleration feasibility or physical tracking.
|
||||
"""
|
||||
peak = 1.5 if curve == "cubic" else 1.875
|
||||
return float(max(500., peak*abs(rise_m)*reference_speed_m_s/(.8*max_climb_rate_m_s)))
|
||||
@@ -1,221 +0,0 @@
|
||||
"""Estimated constant-velocity forecasts and reduced-model reference rollouts.
|
||||
|
||||
The live controller uses ``predict_track`` with measured states only. Legacy
|
||||
point-mass rollout helpers remain available for numerical/reference checks;
|
||||
they must not be fed hidden attacker parameters by the defense controller.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
import numpy as np
|
||||
|
||||
from ..entities.projectile import Projectile
|
||||
from ..physics.integrator import NUMBA_AVAILABLE, rk4_fast
|
||||
from .activation import CityCircle
|
||||
from ..physics.events import ground_fraction, circle_entry_fraction
|
||||
|
||||
if NUMBA_AVAILABLE:
|
||||
from ..physics.integrator import projectile_rk4_kernel
|
||||
else: # pragma: no cover - numba optional
|
||||
from ..entities.projectile import _drag_scalar
|
||||
from ..physics.atmosphere import isa
|
||||
|
||||
def projectile_rk4_kernel(state, wind, drag_k, thrust):
|
||||
rel = state[3:] - wind
|
||||
return np.concatenate((state[3:], np.array([0.0, 0.0, 9.80665])
|
||||
+ isa(-state[2]).density_kg_m3 * _drag_scalar(rel, drag_k) + thrust))
|
||||
|
||||
PREDICTION_DT_S = 0.25
|
||||
PREDICTION_HORIZON_S = 600.0
|
||||
MIN_LINE_SAMPLES = 6
|
||||
|
||||
|
||||
@dataclass
|
||||
class TrajectoryPrediction:
|
||||
"""Best-fit threat trajectory and intercept geometry for one threat."""
|
||||
|
||||
track_time_s: float
|
||||
track_position_ned_m: np.ndarray
|
||||
track_velocity_ned_m_s: np.ndarray
|
||||
line_origin_ned_m: np.ndarray
|
||||
line_direction_ned: np.ndarray # unit vector
|
||||
crossing_ned_m: np.ndarray | None # primary intercept aim (first feasible crossing)
|
||||
crossing_time_s: float | None
|
||||
impact_ned_m: np.ndarray | None # predicted ground contact, if before the crossing
|
||||
impact_time_s: float | None
|
||||
samples: np.ndarray # (n, 4) -> time_s, x, y, z
|
||||
time_to_crossing_s: float | None
|
||||
|
||||
def point_at(self, delta_s: float) -> np.ndarray:
|
||||
"""Point on the threat line delta_s seconds from the crossing aim."""
|
||||
if self.crossing_ned_m is None:
|
||||
return self.line_origin_ned_m
|
||||
return self.crossing_ned_m + self.track_velocity_ned_m_s * delta_s
|
||||
|
||||
|
||||
def predict_track(position, velocity, city, time_s, ceiling_m=3000., horizon_s=600.):
|
||||
"""Declared constant-velocity estimate; no aircraft truth or true wind.
|
||||
|
||||
Forecast uncertainty grows with time. This extrapolation is useful for
|
||||
displaying estimated motion and is not a certified flight prediction.
|
||||
"""
|
||||
position, velocity = np.asarray(position, dtype=float), np.asarray(velocity, dtype=float)
|
||||
end = position + horizon_s*velocity
|
||||
ground = ground_fraction(position, end)
|
||||
entry = circle_entry_fraction(position, end, city.center_ned_m, city.radius_m)
|
||||
crossing = crossing_t = impact = impact_t = None
|
||||
if ground is not None:
|
||||
impact = position + ground*(end-position)
|
||||
impact_t = time_s+ground*horizon_s
|
||||
if entry is not None and (ground is None or entry < ground):
|
||||
point = position+entry*(end-position)
|
||||
if -ceiling_m <= point[2] < 0:
|
||||
crossing, crossing_t = point, time_s+entry*horizon_s
|
||||
speed = float(np.linalg.norm(velocity))
|
||||
if speed < 1e-6:
|
||||
return None
|
||||
duration = min(horizon_s, ground*horizon_s) if ground is not None else horizon_s
|
||||
times = np.linspace(0, duration, 61)
|
||||
samples = np.column_stack((times, position + times[:, None]*velocity))
|
||||
return TrajectoryPrediction(time_s, position.copy(), velocity.copy(), position.copy(), velocity/speed,
|
||||
crossing, crossing_t, impact, impact_t, samples,
|
||||
crossing_t-time_s if crossing_t is not None else None)
|
||||
|
||||
|
||||
def roll_out_threat(projectile: Projectile, position_ned_m: np.ndarray,
|
||||
velocity_ned_m_s: np.ndarray, wind_ned_m_s: np.ndarray,
|
||||
dt: float = PREDICTION_DT_S, horizon_s: float = PREDICTION_HORIZON_S) -> np.ndarray:
|
||||
"""Integrate the observed state forward; returns (n, 4) [t, x, y, z] rows.
|
||||
|
||||
The whole horizon is integrated step-at-a-time with the same RK4 kernel
|
||||
the real projectile uses (single-step fast path; the batch variant cannot
|
||||
hold fixed args across a rollout whose altitude/drag changes each step),
|
||||
and the rollout truncates at the first ground contact, exactly like the
|
||||
previous loop.
|
||||
"""
|
||||
wind = np.asarray(wind_ned_m_s, dtype=np.float64)
|
||||
drag_k = projectile.config.drag_coefficient * projectile.config.area_m2 / projectile.config.mass_kg
|
||||
thrust = projectile._thrust_cmd if projectile.guidance is not None and projectile.guidance_active else np.zeros(3)
|
||||
count = max(1, int(horizon_s / dt))
|
||||
state = np.ascontiguousarray(np.concatenate((position_ned_m, velocity_ned_m_s)), dtype=np.float64)
|
||||
rows = [np.concatenate(([0.0], position_ned_m))]
|
||||
for i in range(count):
|
||||
previous = state.copy()
|
||||
state = rk4_fast(state, dt, projectile_rk4_kernel, wind, drag_k, thrust)
|
||||
p = state[:3]
|
||||
ground = ground_fraction(previous[:3], p)
|
||||
if ground is not None:
|
||||
p = previous[:3] + ground*(p-previous[:3])
|
||||
p[2] = 0.0
|
||||
rows.append(np.concatenate(([(i+ground)*dt], p)))
|
||||
break
|
||||
rows.append(np.concatenate(([(i+1)*dt], p)))
|
||||
return np.asarray(rows, dtype=float)
|
||||
|
||||
|
||||
def fit_line(samples: np.ndarray) -> tuple[np.ndarray, np.ndarray] | None:
|
||||
"""Least-squares line through the rollout samples.
|
||||
|
||||
Model: x(t) = origin + direction * t with |direction| = 1, origin at t=0.
|
||||
Linear in origin and (unnormalized) direction: solved with lstsq.
|
||||
Returns (origin, unit_direction) or None if the fit is degenerate.
|
||||
"""
|
||||
if samples.shape[0] < MIN_LINE_SAMPLES:
|
||||
return None
|
||||
t = samples[:, 0]
|
||||
x = samples[:, 1:]
|
||||
t = t - t[0]
|
||||
# center the times for conditioning
|
||||
tc = t - float(t.mean())
|
||||
matrix = np.column_stack((np.ones_like(tc), tc[:, None]))
|
||||
solution, residuals, rank, _ = np.linalg.lstsq(matrix, x, rcond=None)
|
||||
origin = solution[0] - solution[1] * float(t.mean())
|
||||
direction = solution[1]
|
||||
norm = float(np.linalg.norm(direction))
|
||||
if norm < 1e-9:
|
||||
return None
|
||||
return origin, direction / norm
|
||||
|
||||
|
||||
def predict_trajectory(projectile: Projectile, position_ned_m: np.ndarray,
|
||||
velocity_ned_m_s: np.ndarray, wind_ned_m_s: np.ndarray,
|
||||
city: CityCircle, loiter_ceiling_m: float,
|
||||
track_time_s: float,
|
||||
dt: float = PREDICTION_DT_S, horizon_s: float = PREDICTION_HORIZON_S
|
||||
) -> TrajectoryPrediction | None:
|
||||
"""Full prediction: rollout + line fit + 3D protected-circle crossing."""
|
||||
samples = roll_out_threat(projectile, position_ned_m, velocity_ned_m_s, wind_ned_m_s, dt, horizon_s)
|
||||
fit = fit_line(samples)
|
||||
if fit is None:
|
||||
return None
|
||||
origin, direction = fit
|
||||
|
||||
center = np.asarray(city.center_ned_m, dtype=float)
|
||||
radius = city.radius_m
|
||||
crossing: np.ndarray | None = None
|
||||
crossing_t: float | None = None
|
||||
impact: np.ndarray | None = None
|
||||
impact_t: float | None = None
|
||||
|
||||
previous = samples[0, 1:]
|
||||
previous_t = 0.0
|
||||
for row in samples[1:]:
|
||||
t, current = float(row[0]), row[1:]
|
||||
entry = circle_entry_fraction(previous, current, center, radius)
|
||||
ground = ground_fraction(previous, current)
|
||||
if entry is not None and (ground is None or entry < ground):
|
||||
point = previous + entry*(current-previous)
|
||||
if -loiter_ceiling_m <= point[2] < 0:
|
||||
crossing = point
|
||||
crossing_t = track_time_s + previous_t + entry*(t-previous_t)
|
||||
break
|
||||
if ground is not None:
|
||||
impact = previous + ground*(current-previous)
|
||||
impact_t = track_time_s + previous_t + ground*(t-previous_t)
|
||||
break
|
||||
previous, previous_t = current, t
|
||||
time_to_crossing = crossing_t-track_time_s if crossing_t is not None else None
|
||||
|
||||
return TrajectoryPrediction(
|
||||
track_time_s=track_time_s,
|
||||
track_position_ned_m=np.asarray(position_ned_m, dtype=float),
|
||||
track_velocity_ned_m_s=np.asarray(velocity_ned_m_s, dtype=float),
|
||||
line_origin_ned_m=origin,
|
||||
line_direction_ned=direction,
|
||||
crossing_ned_m=crossing,
|
||||
crossing_time_s=crossing_t,
|
||||
impact_ned_m=impact,
|
||||
impact_time_s=impact_t,
|
||||
samples=samples,
|
||||
time_to_crossing_s=time_to_crossing,
|
||||
)
|
||||
|
||||
|
||||
def intercept_waypoints(prediction: TrajectoryPrediction, count: int,
|
||||
standoff_m: float, order: list[str]) -> dict[str, np.ndarray]:
|
||||
"""Waypoint per drone: lead ahead of the aim, middle on it, trail behind.
|
||||
|
||||
`order` lists drone IDs from nearest-to-farthest of the aim point; the
|
||||
nearest (lead) flies to the point standoff_m *before* the crossing along
|
||||
the threat line, the middle one walks the crossing point itself, and the
|
||||
trail one covers standoff_m *after* it so the threat cannot slip through.
|
||||
"""
|
||||
if prediction.crossing_ned_m is None:
|
||||
origin = prediction.line_origin_ned_m
|
||||
return {drone_id: origin.copy() for drone_id in order}
|
||||
direction = prediction.line_direction_ned
|
||||
# orient direction so that positive offset = toward the crossing
|
||||
# (the line was fit over the whole rollout; the crossing is usually the
|
||||
# latest feasible sample, so the direction points toward it already when
|
||||
# the origin is the track position). Guard against the reversed case.
|
||||
vector_to_crossing = prediction.crossing_ned_m - prediction.line_origin_ned_m
|
||||
if float(np.dot(vector_to_crossing, direction)) < 0.0:
|
||||
direction = -direction
|
||||
offsets = []
|
||||
n = len(order)
|
||||
for index in range(n):
|
||||
# nearest drone -> largest standoff before the aim, farthest -> closest
|
||||
remaining = n - 1 - index
|
||||
offset = -standoff_m * remaining if n > 1 else 0.0
|
||||
offsets.append(offset)
|
||||
return {drone_id: prediction.crossing_ned_m + direction * offset for drone_id, offset in zip(order, offsets)}
|
||||
@@ -1,165 +0,0 @@
|
||||
"""Device profiles and timing taken from this repository's AERIS-10 design.
|
||||
|
||||
The mechanical/electronic scan cadence is derived from the STM32 firmware's
|
||||
runRadarPulseSequence() constants. Measurement uncertainty and a 7.1 ms fixed
|
||||
processing delay follow simulate_trajectory_precision.py and remain explicit
|
||||
R&D assumptions, not bench-qualified radar performance.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
import numpy as np
|
||||
from ..physics.validation import finite_numeric_fields
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AerisRadarConfig:
|
||||
model: str = "aeris_10x"
|
||||
station_count: int = 8
|
||||
ring_radius_km: float = 12.0
|
||||
range_sigma_m: float = 2.0
|
||||
azimuth_sigma_deg: float = 0.5
|
||||
elevation_sigma_deg: float = 0.5
|
||||
doppler_sigma_m_s: float = 1.0
|
||||
wired_processing_latency_s: float = 0.0071
|
||||
half_sector_deg: float = 3.6
|
||||
elevation_limit_deg: float = 45.0
|
||||
elevation_step_deg: float = 3.0
|
||||
elevation_beam_half_width_deg: float = 1.5
|
||||
|
||||
@property
|
||||
def max_range_m(self) -> float:
|
||||
return 20_000.0 if self.model == "aeris_10x" else 3_000.0
|
||||
|
||||
@property
|
||||
def azimuth_step_period_s(self) -> float:
|
||||
# Firmware constants: 16 chirps at each PRI, 175.4 us guard,
|
||||
# three phase patterns x 15 elevation sweeps, four 1 ms motor steps.
|
||||
pulse_train_s = 16 * (167e-6 + 175e-6) + 175.4e-6
|
||||
return 15 * (3 * pulse_train_s) + 4e-3
|
||||
|
||||
@property
|
||||
def full_azimuth_revisit_s(self) -> float:
|
||||
return 50 * self.azimuth_step_period_s
|
||||
|
||||
@property
|
||||
def full_scan_update_hz(self) -> float:
|
||||
return 1.0 / self.full_azimuth_revisit_s
|
||||
|
||||
def __post_init__(self):
|
||||
finite_numeric_fields(self)
|
||||
if any(value < 0 for value in (self.range_sigma_m, self.azimuth_sigma_deg, self.elevation_sigma_deg, self.doppler_sigma_m_s, self.wired_processing_latency_s)):
|
||||
raise ValueError("radar noise and latency must be nonnegative")
|
||||
if self.elevation_step_deg <= 0 or self.half_sector_deg <= 0 or self.elevation_beam_half_width_deg <= 0:
|
||||
raise ValueError("radar beam geometry must be positive")
|
||||
if self.model not in {"aeris_10n", "aeris_10x"}:
|
||||
raise ValueError("radar model must be aeris_10n or aeris_10x")
|
||||
if not 1 <= self.station_count <= 16 or self.ring_radius_km <= 0:
|
||||
raise ValueError("radar station count/ring radius are outside supported bounds")
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class DroneLinkConfig:
|
||||
update_rate_hz: float = 20.0
|
||||
one_way_latency_s: float = 0.050
|
||||
latency_jitter_sigma_s: float = 0.005
|
||||
command_compute_latency_s: float = 0.010
|
||||
packet_loss_probability: float = 0.0
|
||||
bitrate_bps: float = 250_000.0
|
||||
packet_overhead_bytes: int = 24
|
||||
telemetry_payload_bytes: int = 64
|
||||
command_payload_bytes: int = 48
|
||||
|
||||
def __post_init__(self):
|
||||
finite_numeric_fields(self)
|
||||
if self.update_rate_hz <= 0 or self.one_way_latency_s < 0 or self.latency_jitter_sigma_s < 0 or self.bitrate_bps <= 0:
|
||||
raise ValueError("invalid drone datalink timing")
|
||||
if self.command_compute_latency_s < 0 or min(self.packet_overhead_bytes, self.telemetry_payload_bytes, self.command_payload_bytes) < 0:
|
||||
raise ValueError("invalid packet sizes/processing delay")
|
||||
if not 0 <= self.packet_loss_probability <= 1:
|
||||
raise ValueError("packet loss probability must be in [0, 1]")
|
||||
|
||||
|
||||
class Aeris10Radar:
|
||||
"""One scanning AERIS-10 site with a quantized 50-position azimuth scan."""
|
||||
|
||||
def __init__(self, radar_id: str, position_ned_m: np.ndarray, inward_heading_rad: float,
|
||||
config: AerisRadarConfig, rng: np.random.Generator, initial_slot: int = 0):
|
||||
self.radar_id = radar_id
|
||||
self.position_ned_m = np.asarray(position_ned_m, dtype=float)
|
||||
self.inward_heading_rad = float(inward_heading_rad)
|
||||
self.config = config
|
||||
self.rng = rng
|
||||
self.azimuth_slot = int(initial_slot) % 50
|
||||
self.next_sweep_time_s = 0.0
|
||||
self.last_detection_time_s: float | None = None
|
||||
|
||||
def scan(self, target, time_s: float) -> dict | None:
|
||||
measurements = self.scan_dwell([target], time_s)
|
||||
return measurements[0] if measurements else None
|
||||
|
||||
def is_dwell_due(self, time_s: float) -> bool:
|
||||
return time_s + 1e-9 >= self.next_sweep_time_s
|
||||
|
||||
def scan_dwell(self, targets, time_s: float) -> list[dict]:
|
||||
"""One dwell pass over every candidate target.
|
||||
|
||||
The slot advances once per dwell and each active target is evaluated
|
||||
with the same beam/time state; a first target that is out of range or
|
||||
out of beam can no longer consume the opportunity for later ones (F03).
|
||||
Detection remains deterministic within geometry/range constraints: no
|
||||
cross-section/SNR/clutter model is claimed here (R&D assumption).
|
||||
"""
|
||||
if not self.is_dwell_due(time_s):
|
||||
return []
|
||||
period = self.config.azimuth_step_period_s
|
||||
missed = max(0, int(np.floor((time_s - self.next_sweep_time_s + 1e-9) / period)))
|
||||
self.azimuth_slot = (self.azimuth_slot + missed) % 50
|
||||
self.next_sweep_time_s += (missed + 1) * period
|
||||
slot_heading = self.inward_heading_rad + self.azimuth_slot * (2 * np.pi / 50)
|
||||
self.azimuth_slot = (self.azimuth_slot + 1) % 50
|
||||
measurements = []
|
||||
for target in targets:
|
||||
measurement = self._evaluate_at_dwell(target, time_s, slot_heading)
|
||||
if measurement is not None:
|
||||
measurement["target_id"] = target.entity_id
|
||||
measurements.append(measurement)
|
||||
return measurements
|
||||
|
||||
def _evaluate_at_dwell(self, target, time_s: float, slot_heading: float) -> dict | None:
|
||||
relative = target.position_ned_m - self.position_ned_m
|
||||
distance = float(np.linalg.norm(relative))
|
||||
if distance <= 1e-6 or distance > self.config.max_range_m:
|
||||
return None
|
||||
north, east, down = relative
|
||||
azimuth = float(np.arctan2(east, north))
|
||||
azimuth_error = (azimuth - slot_heading + np.pi) % (2 * np.pi) - np.pi
|
||||
if abs(azimuth_error) > np.deg2rad(self.config.half_sector_deg):
|
||||
return None
|
||||
elevation = float(np.arctan2(-down, np.hypot(north, east)))
|
||||
if abs(elevation) > np.deg2rad(self.config.elevation_limit_deg):
|
||||
return None
|
||||
nearest_elevation_center = round(np.rad2deg(elevation) / self.config.elevation_step_deg) * self.config.elevation_step_deg
|
||||
if abs(np.rad2deg(elevation) - nearest_elevation_center) > self.config.elevation_beam_half_width_deg:
|
||||
return None
|
||||
range_measured = distance + self.rng.normal(0, self.config.range_sigma_m)
|
||||
azimuth_measured = azimuth + self.rng.normal(0, np.deg2rad(self.config.azimuth_sigma_deg))
|
||||
elevation_measured = elevation + self.rng.normal(0, np.deg2rad(self.config.elevation_sigma_deg))
|
||||
direction = np.array([np.cos(elevation_measured) * np.cos(azimuth_measured),
|
||||
np.cos(elevation_measured) * np.sin(azimuth_measured),
|
||||
-np.sin(elevation_measured)])
|
||||
position = self.position_ned_m + max(0.0, range_measured) * direction
|
||||
line_of_sight = direction
|
||||
radial_speed = float(np.dot(target.velocity_ned_m_s, relative / distance)
|
||||
+ self.rng.normal(0, self.config.doppler_sigma_m_s))
|
||||
self.last_detection_time_s = time_s
|
||||
return {"radar_id": self.radar_id, "measurement_time_s": time_s,
|
||||
"arrival_time_s": time_s + self.config.wired_processing_latency_s,
|
||||
"position_ned_m": position, "line_of_sight_ned": line_of_sight,
|
||||
"radial_velocity_m_s": radial_speed,
|
||||
"sensor_position_ned_m": self.position_ned_m.copy(),
|
||||
"range_sigma_m": self.config.range_sigma_m,
|
||||
"azimuth_sigma_rad": np.deg2rad(self.config.azimuth_sigma_deg),
|
||||
"elevation_sigma_rad": np.deg2rad(self.config.elevation_sigma_deg),
|
||||
"range_m": range_measured, "azimuth_rad": azimuth_measured,
|
||||
"elevation_rad": elevation_measured}
|
||||
@@ -1,782 +0,0 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
import uuid
|
||||
from copy import deepcopy
|
||||
import math
|
||||
import numpy as np
|
||||
from ..entities.drone import Drone
|
||||
from ..entities.projectile import Projectile
|
||||
from ..physics.environment import Environment
|
||||
from ..sensors.tracking import RadarTracker
|
||||
from ..sensors.models import Observation, Sensor, SensorConfig
|
||||
from ..guidance.activation import CityCircle, select_defenders, reachable_within
|
||||
from ..guidance.ascent import automatic_ascent_length
|
||||
from ..guidance.trajectory import predict_track
|
||||
from ..physics.atmosphere import isa
|
||||
from ..physics.events import sphere_entry_fraction, ground_fraction, circle_entry_fraction
|
||||
from ..sensors.devices import Aeris10Radar, AerisRadarConfig, DroneLinkConfig
|
||||
|
||||
|
||||
@dataclass
|
||||
class SimulationEngine:
|
||||
dt: float = 0.01
|
||||
seed: int | None = None
|
||||
environment: Environment = field(default_factory=Environment)
|
||||
city_circle: CityCircle | None = None
|
||||
sensor_config: SensorConfig = field(default_factory=SensorConfig)
|
||||
radar_config: AerisRadarConfig = field(default_factory=AerisRadarConfig)
|
||||
drone_link_config: DroneLinkConfig = field(default_factory=DroneLinkConfig)
|
||||
drone_control_rate_hz: float = 50.0
|
||||
projectiles: list[Projectile] = field(default_factory=list)
|
||||
drones: list[Drone] = field(default_factory=list)
|
||||
time_s: float = 0.0
|
||||
run_id: str = field(default_factory=lambda: str(uuid.uuid4()))
|
||||
events: list[dict] = field(default_factory=list)
|
||||
trajectory: dict = field(default_factory=dict)
|
||||
closest_approach_m: dict = field(default_factory=dict)
|
||||
record_interval_s: float = 0.2
|
||||
|
||||
def __post_init__(self):
|
||||
if not np.isfinite(self.dt) or self.dt <= 0:
|
||||
raise ValueError("simulation dt must be positive")
|
||||
if not np.isfinite(self.record_interval_s) or self.record_interval_s <= 0:
|
||||
raise ValueError("trajectory record interval must be positive")
|
||||
ids = [e.entity_id for e in [*self.projectiles, *self.drones]]
|
||||
if len(set(ids)) != len(ids):
|
||||
raise ValueError("entity ids must be unique")
|
||||
for entity in [*self.projectiles, *self.drones]:
|
||||
for name in ("position_ned_m", "velocity_ned_m_s"):
|
||||
value = np.asarray(getattr(entity, name), dtype=float)
|
||||
if value.shape != (3,) or not np.isfinite(value).all():
|
||||
raise ValueError(f"{name} must be a finite 3-vector")
|
||||
setattr(entity, name, value.copy())
|
||||
self.rng = np.random.default_rng(self.seed)
|
||||
if not np.isfinite(self.drone_control_rate_hz) or not 0 < self.drone_control_rate_hz <= 1/self.dt:
|
||||
raise ValueError("drone control rate must be positive")
|
||||
self.sensors = {d.entity_id: Sensor(self.sensor_config, np.random.default_rng(self.rng.integers(0, 2**63))) for d in self.drones}
|
||||
self.latest_observations = {}
|
||||
self.radar_tracks = {}
|
||||
self.radar_tracker = RadarTracker()
|
||||
self.track_predictions = {}
|
||||
self.radar_arrivals = []
|
||||
self.latest_detections: dict[str, dict] = {}
|
||||
self.radar_detection_points = {}
|
||||
self.pending_drone_commands = []
|
||||
self.standby_launches = {}
|
||||
self.pending_drone_telemetry = []
|
||||
self.wireless_channel_free_time_s = 0.0
|
||||
self.wireless_airtime_s = 0.0
|
||||
self.remote_drone_estimates = {}
|
||||
self.remote_drone_velocities = {}
|
||||
self.remote_drone_estimate_times = {}
|
||||
self.remote_drone_receive_times = {}
|
||||
self.remote_drone_link_latencies = {}
|
||||
self.remote_drone_telemetry_age_s = {}
|
||||
self.next_control_time_s = {d.entity_id: 0.0 for d in self.drones}
|
||||
self.activated_threats: set[str] = set()
|
||||
self._last_activation_check_s: dict[str, float] = {}
|
||||
for e in [*self.projectiles, *self.drones]:
|
||||
self.trajectory[e.entity_id] = [np.asarray(e.position_ned_m, dtype=float).copy()]
|
||||
self.closest_approach_m = {p.entity_id: float("inf") for p in self.projectiles}
|
||||
for p in self.projectiles:
|
||||
p.estimated_position_ned_m = None
|
||||
p.estimated_velocity_ned_m_s = None
|
||||
# Per-projectile defense lifecycle bookkeeping (assignment state
|
||||
# machine: proposed -> sent -> acknowledged -> active -> released).
|
||||
self.protected_entries: set[str] = set()
|
||||
self.projectile_outcomes: dict[str, str] = {}
|
||||
self.defense_assignments: dict[str, dict] = {} # threat_id -> assignment record
|
||||
self.pending_assignment_commands_by_drone: dict[str, list[dict]] = {}
|
||||
self.threat_release_events: list[tuple[float, str]] = [] # (time_s, threat_id)
|
||||
self.radars = []
|
||||
if self.city_circle is not None:
|
||||
ring_m = self.radar_config.ring_radius_km * 1000.0
|
||||
for index in range(self.radar_config.station_count):
|
||||
angle = 2 * np.pi * index / self.radar_config.station_count
|
||||
position = self.city_circle.center_ned_m + np.array([ring_m * np.cos(angle), ring_m * np.sin(angle), 0.0])
|
||||
radar_rng = np.random.default_rng(self.rng.integers(0, 2**63))
|
||||
self.radars.append(Aeris10Radar(f"aeris10_{index+1:02d}", position, angle + np.pi,
|
||||
self.radar_config, radar_rng, initial_slot=int(radar_rng.integers(0, 50))))
|
||||
|
||||
def step(self) -> "SimulationEngine":
|
||||
"""Advance one physics step. Never serializes; use ``step_and_snapshot``
|
||||
when a JSON-able state is needed."""
|
||||
if not hasattr(self, "record_interval_s"):
|
||||
self.record_interval_s = self.dt
|
||||
# Observation -> controller -> actuator/physics; each controller only sees estimates.
|
||||
moving = [p for p in self.projectiles if p.active] + [d for d in self.drones if d.active and (d.position_ned_m[2] < 0 or d.target_position_ned_m is not None)]
|
||||
if moving:
|
||||
altitude = float(np.mean([-e.position_ned_m[2] for e in moving]))
|
||||
speed = float(np.mean([np.linalg.norm(e.velocity_ned_m_s-self.environment.wind_ned_m_s) for e in moving]))
|
||||
else:
|
||||
altitude, speed = 0., 0.
|
||||
self.environment.update(self.dt, altitude, speed)
|
||||
for d in self.drones:
|
||||
if not d.active:
|
||||
continue
|
||||
sensor = self.sensors[d.entity_id]
|
||||
# Nav samples run at their own cadence (default 10 Hz): the truth
|
||||
# history is cheap to keep every step, but observe() (and with it
|
||||
# telemetry queuing and RNG draws) only fires when due.
|
||||
if not sensor.is_update_due(self.time_s) and not sensor.pending:
|
||||
continue
|
||||
obs = sensor.observe(d, self.time_s)
|
||||
if obs is not None:
|
||||
self.latest_observations[d.entity_id] = obs
|
||||
self._queue_drone_telemetry(d.entity_id, obs)
|
||||
# Radar scan cadence: a radar advances its azimuth dwell at most once
|
||||
# per physics step (its 0.258 s slot period is far longer than dt), and
|
||||
# one dwell evaluates ALL active targets with the same beam/time state
|
||||
# before the slot advances -- target list order must not decide which
|
||||
# threats get a chance to be detected (F03).
|
||||
self._scan_radars()
|
||||
self._deliver_radar_tracks()
|
||||
self._deliver_drone_telemetry()
|
||||
self._service_standby_launches()
|
||||
projectile_tracks = dict(self.radar_tracks)
|
||||
self._service_assignments()
|
||||
self._activate_defenders(projectile_tracks)
|
||||
self._deliver_drone_commands()
|
||||
wind = self.environment.air_velocity_ned_m_s
|
||||
previous = {e.entity_id: e.position_ned_m.copy() for e in [*self.projectiles, *self.drones]}
|
||||
previous_velocity = {e.entity_id: e.velocity_ned_m_s.copy() for e in [*self.projectiles, *self.drones]}
|
||||
previously_active = {e.entity_id: e.active for e in [*self.projectiles, *self.drones]}
|
||||
for p in self.projectiles:
|
||||
if p.guidance is not None and not p.guidance_active and p.active:
|
||||
p.guidance_active = True
|
||||
self.events.append({"type": "terminal_guidance_active", "time_s": self.time_s,
|
||||
"projectile_id": p.entity_id,
|
||||
"strike_target_ned_m": p.guidance.target_ned_m.tolist(),
|
||||
"note": "terminal glide steering engaged on the strike target"})
|
||||
p.step(self.dt, wind, self.time_s)
|
||||
control_period = 1.0 / self.drone_control_rate_hz
|
||||
record_every = max(1, int(round(self.record_interval_s / self.dt))) if hasattr(self, "record_interval_s") else 1
|
||||
step_count = int(round(self.time_s / self.dt))
|
||||
for d in self.drones:
|
||||
control_due = self.time_s + 1e-9 >= self.next_control_time_s[d.entity_id]
|
||||
if control_due:
|
||||
while self.next_control_time_s[d.entity_id] <= self.time_s + 1e-9:
|
||||
self.next_control_time_s[d.entity_id] += control_period
|
||||
d.step(self.dt, wind, self.rng, control_due=control_due, control_dt=control_period)
|
||||
self.time_s += self.dt
|
||||
self._resolve_physical_events(previous, previous_velocity)
|
||||
for drone in self.drones:
|
||||
if drone.active and drone.activated_for is not None:
|
||||
threat = next((p for p in self.projectiles if p.entity_id == drone.activated_for), None)
|
||||
drone.estimate_threat(threat, self.time_s, self.rng, self.events,
|
||||
radar_prior=self.radar_tracks.get(drone.activated_for))
|
||||
self._release_finished_assignments()
|
||||
record_now = step_count % record_every == 0
|
||||
for e in [*self.projectiles, *self.drones]:
|
||||
if record_now or (previously_active[e.entity_id] and not e.active):
|
||||
self.trajectory[e.entity_id].append(e.position_ned_m.copy())
|
||||
del self.trajectory[e.entity_id][:-1800]
|
||||
return self
|
||||
|
||||
def step_and_snapshot(self) -> dict:
|
||||
"""Advance one physics step and return a serialized snapshot.
|
||||
|
||||
``step`` itself never serializes (that was the dominant per-step cost);
|
||||
callers that need a JSON-able state use this wrapper or ``snapshot``.
|
||||
"""
|
||||
self.step()
|
||||
return self.snapshot()
|
||||
|
||||
def _scan_radars(self):
|
||||
"""One dwell per radar per step, evaluating all active targets together."""
|
||||
targets = [p for p in self.projectiles if p.active]
|
||||
if not targets:
|
||||
return
|
||||
for radar in self.radars:
|
||||
for measurement in radar.scan_dwell(targets, self.time_s):
|
||||
target_id = measurement["target_id"]
|
||||
self.radar_arrivals.append(measurement)
|
||||
history = self.radar_detection_points.setdefault(target_id, [])
|
||||
history.append(measurement["position_ned_m"].tolist())
|
||||
del history[:-600]
|
||||
prev = self.latest_detections.get(target_id)
|
||||
if prev is None or measurement["measurement_time_s"] > prev["measurement_time_s"]:
|
||||
self.latest_detections[target_id] = measurement
|
||||
|
||||
def _deliver_radar_tracks(self):
|
||||
due = [m for m in self.radar_arrivals if m["arrival_time_s"] <= self.time_s + 1e-9]
|
||||
self.radar_arrivals = [m for m in self.radar_arrivals if m["arrival_time_s"] > self.time_s + 1e-9]
|
||||
targets = {p.entity_id: p for p in self.projectiles if p.active}
|
||||
for measurement in sorted(due, key=lambda m: (m["measurement_time_s"], m["radar_id"])):
|
||||
target = targets.get(measurement["target_id"])
|
||||
if target is None:
|
||||
continue
|
||||
# One delivered position detection is enough for standby launch;
|
||||
# velocity evidence and reachability gate only the later assignment.
|
||||
self._request_standby_launches(target.entity_id)
|
||||
estimate = self.radar_tracker.update(measurement)
|
||||
if estimate is None:
|
||||
continue
|
||||
# Truth comparison is diagnostic only and never enters the filter.
|
||||
obs = Observation(target.entity_id, estimate.timestamp_s, estimate.position,
|
||||
estimate.velocity, target.position_ned_m-estimate.position,
|
||||
target.velocity_ned_m_s-estimate.velocity, estimate.covariance,
|
||||
estimate.sample_count, estimate.velocity_observable,
|
||||
error_reference="truth_at_receive_minus_sample_time_estimate")
|
||||
if self.city_circle is not None and estimate.velocity_observable:
|
||||
prediction = predict_track(estimate.position, estimate.velocity, self.city_circle, estimate.timestamp_s)
|
||||
if prediction is not None:
|
||||
self.track_predictions[target.entity_id] = prediction.samples
|
||||
self.radar_tracks[target.entity_id] = obs
|
||||
self.latest_observations[target.entity_id] = obs
|
||||
target.estimated_position_ned_m = estimate.position.copy()
|
||||
target.estimated_velocity_ned_m_s = estimate.velocity.copy()
|
||||
self.events.append({"type": "radar_track_update", "time_s": self.time_s,
|
||||
"measurement_time_s": estimate.timestamp_s, "target_id": target.entity_id,
|
||||
"radar_ids": [measurement["radar_id"]], "position_ned_m": estimate.position.tolist(),
|
||||
"velocity_ned_m_s": estimate.velocity.tolist(), "covariance": estimate.covariance.tolist(),
|
||||
"sample_count": estimate.sample_count, "velocity_observable": estimate.velocity_observable,
|
||||
"velocity_estimation_rule": "stonesoup_constant_velocity_position_filter"})
|
||||
|
||||
def _deliver_drone_commands(self):
|
||||
pending = []
|
||||
for command in self.pending_drone_commands:
|
||||
if command["delivery_time_s"] > self.time_s + 1e-9:
|
||||
pending.append(command)
|
||||
continue
|
||||
drone = next((d for d in self.drones if d.entity_id == command["drone_id"]), None)
|
||||
if drone is None:
|
||||
continue
|
||||
target = next((p for p in self.projectiles if p.entity_id == command["threat_id"]), None)
|
||||
if command.get("kind") == "standby_takeoff":
|
||||
if not drone.active or drone.activated_for is not None or target is None or not target.active:
|
||||
continue
|
||||
drone.standby_altitude_m = max(-drone.position_ned_m[2], command["altitude_m"])
|
||||
drone.target_position_ned_m = drone.position_ned_m.copy()
|
||||
self.standby_launches[drone.entity_id]["status"] = "active"
|
||||
self.events.append({"type": "standby_takeoff_started", "time_s": self.time_s,
|
||||
"drone_id": drone.entity_id, "threat_id": target.entity_id,
|
||||
"standby_altitude_m": drone.standby_altitude_m,
|
||||
"end_to_end_command_delay_s": self.time_s-command["computed_time_s"]})
|
||||
continue
|
||||
record = self.defense_assignments.get(command["threat_id"])
|
||||
if not drone.active or target is None or not target.active or record is None:
|
||||
continue
|
||||
if command["computed_time_s"] <= getattr(drone, "_last_command_computed_time_s", -float("inf")):
|
||||
continue
|
||||
drone._last_command_computed_time_s = command["computed_time_s"]
|
||||
drone.activated_for = command["threat_id"]
|
||||
drone.standby_altitude_m = None
|
||||
if drone.flight_reference is None:
|
||||
drone.target_position_ned_m = drone.position_ned_m.copy()
|
||||
reference = command.get("flight_reference")
|
||||
if reference is not None:
|
||||
if drone.flight_reference is None:
|
||||
drone.flight_reference = {**reference,
|
||||
"initial_altitude_m": -drone.position_ned_m[2], "travel_m": 0., "elapsed_s": 0.,
|
||||
"previous_position": drone.position_ned_m.copy()}
|
||||
else:
|
||||
# Preserve ascent progress and elapsed phase during replans.
|
||||
drone.flight_reference["direction"] = reference["direction"].copy()
|
||||
drone.flight_reference["target_altitude_m"] = reference["target_altitude_m"]
|
||||
else:
|
||||
drone.target_position_ned_m = command["target_position_ned_m"].copy()
|
||||
record["status"][drone.entity_id] = "active"
|
||||
if all(value == "active" for value in record["status"].values()) and not record.get("activation_recorded"):
|
||||
record["activation_recorded"] = True
|
||||
self.events.append({"type": "defenders_activated", "time_s": self.time_s,
|
||||
"threat_id": target.entity_id, "defender_ids": record["defender_ids"]})
|
||||
self.events.append({"type": "drone_command_received", "time_s": self.time_s,
|
||||
"command_sent_time_s": command["sent_time_s"], "drone_id": drone.entity_id,
|
||||
"threat_id": command["threat_id"],
|
||||
"end_to_end_command_delay_s": self.time_s - command["computed_time_s"],
|
||||
"wireless_queue_delay_s": command["queue_delay_s"],
|
||||
"wireless_airtime_s": command["airtime_s"],
|
||||
"wireless_one_way_latency_s": self.time_s - command["sent_time_s"]})
|
||||
self.pending_drone_commands = pending
|
||||
|
||||
def _request_standby_launches(self, threat_id):
|
||||
scenario = getattr(self, "scenario_config", {})
|
||||
if not scenario.get("early_takeoff_enabled", True):
|
||||
return
|
||||
for drone in self.drones:
|
||||
if drone.entity_type != "defense_drone" or not drone.active or not drone.available or drone.activated_for is not None:
|
||||
continue
|
||||
if drone.entity_id in self.standby_launches:
|
||||
record = self.standby_launches[drone.entity_id]
|
||||
previous_threat_active = any(p.active and p.entity_id == record["threat_id"] for p in self.projectiles)
|
||||
if drone.standby_altitude_m is not None or record["threat_id"] == threat_id or previous_threat_active:
|
||||
continue
|
||||
# An alert for a terminated threat must not permanently block
|
||||
# a later incident if its command never reached this drone.
|
||||
self.pending_drone_commands = [c for c in self.pending_drone_commands
|
||||
if not (c["drone_id"] == drone.entity_id and c.get("kind") == "standby_takeoff")]
|
||||
self.standby_launches[drone.entity_id] = {"threat_id": threat_id,
|
||||
"altitude_m": scenario.get("standby_altitude_m", 150.), "attempts": 0,
|
||||
"last_send": -float("inf"), "status": "proposed"}
|
||||
self._send_standby_command(drone.entity_id)
|
||||
|
||||
def _send_standby_command(self, drone_id):
|
||||
record = self.standby_launches[drone_id]
|
||||
record["attempts"] += 1
|
||||
record["last_send"] = self.time_s
|
||||
computed = self.time_s+self.drone_link_config.command_compute_latency_s
|
||||
_, sent = self._wireless_transmit(computed, self.drone_link_config.command_payload_bytes)
|
||||
if self.rng.random() < self.drone_link_config.packet_loss_probability:
|
||||
record["status"] = "dropped"
|
||||
self.events.append({"type": "standby_command_dropped", "time_s": self.time_s,
|
||||
"drone_id": drone_id, "threat_id": record["threat_id"]})
|
||||
return
|
||||
latency = max(0., self.drone_link_config.one_way_latency_s+self.rng.normal(
|
||||
0., self.drone_link_config.latency_jitter_sigma_s))
|
||||
record["status"] = "sent"
|
||||
self.pending_drone_commands.append({"kind": "standby_takeoff", "drone_id": drone_id,
|
||||
"threat_id": record["threat_id"], "altitude_m": record["altitude_m"],
|
||||
"computed_time_s": computed, "delivery_time_s": sent+latency})
|
||||
self.events.append({"type": "standby_takeoff_requested", "time_s": self.time_s,
|
||||
"drone_id": drone_id, "threat_id": record["threat_id"], "altitude_m": record["altitude_m"]})
|
||||
|
||||
def _service_standby_launches(self):
|
||||
active_threats = {p.entity_id for p in self.projectiles if p.active}
|
||||
for did, record in self.standby_launches.items():
|
||||
drone = next(d for d in self.drones if d.entity_id == did)
|
||||
if not drone.active or drone.activated_for is not None or record["threat_id"] not in active_threats:
|
||||
if record["status"] != "active":
|
||||
record["status"] = "cancelled"
|
||||
continue
|
||||
if record["status"] in {"active", "failed", "cancelled"}:
|
||||
continue
|
||||
if any(c["drone_id"] == did and c.get("kind") == "standby_takeoff" for c in self.pending_drone_commands):
|
||||
continue
|
||||
if self.time_s-record["last_send"] < 1.:
|
||||
continue
|
||||
if record["attempts"] >= 3:
|
||||
record["status"] = "failed"
|
||||
self.events.append({"type": "standby_takeoff_failed", "time_s": self.time_s,
|
||||
"drone_id": did, "threat_id": record["threat_id"]})
|
||||
else:
|
||||
self._send_standby_command(did)
|
||||
|
||||
def _queue_drone_telemetry(self, drone_id, observation):
|
||||
if self.rng.random() < self.drone_link_config.packet_loss_probability:
|
||||
self.events.append({"type": "drone_telemetry_dropped", "time_s": self.time_s,
|
||||
"drone_id": drone_id})
|
||||
return
|
||||
computed_at = self.time_s
|
||||
request_time = computed_at
|
||||
tx_start, transmitted_at = self._wireless_transmit(request_time, self.drone_link_config.telemetry_payload_bytes)
|
||||
latency = max(0.0, self.drone_link_config.one_way_latency_s
|
||||
+ self.rng.normal(0.0, self.drone_link_config.latency_jitter_sigma_s))
|
||||
self.pending_drone_telemetry.append({"drone_id": drone_id,
|
||||
"measurement_time_s": observation.timestamp_s, "position_ned_m": observation.position_ned_m.copy(),
|
||||
"velocity_ned_m_s": observation.velocity_ned_m_s.copy(), "sent_time_s": float(transmitted_at),
|
||||
"delivery_time_s": float(transmitted_at + latency),
|
||||
"queue_delay_s": max(0.0, tx_start - request_time),
|
||||
"airtime_s": transmitted_at - tx_start, "wireless_latency_s": latency})
|
||||
|
||||
def _wireless_transmit(self, request_time_s: float, payload_bytes: int):
|
||||
"""Serialize packets through one shared R&D half-duplex radio channel."""
|
||||
period = 1.0 / self.drone_link_config.update_rate_hz
|
||||
slot = np.ceil((request_time_s - 1e-9) / period) * period
|
||||
start = max(float(slot), self.wireless_channel_free_time_s)
|
||||
airtime = (payload_bytes + self.drone_link_config.packet_overhead_bytes) * 8.0 / self.drone_link_config.bitrate_bps
|
||||
end = start + airtime
|
||||
self.wireless_channel_free_time_s = end
|
||||
self.wireless_airtime_s += airtime
|
||||
return start, end
|
||||
|
||||
def _deliver_drone_telemetry(self):
|
||||
pending = []
|
||||
for packet in self.pending_drone_telemetry:
|
||||
if packet["delivery_time_s"] > self.time_s + 1e-9:
|
||||
pending.append(packet)
|
||||
continue
|
||||
if packet["measurement_time_s"] <= self.remote_drone_estimate_times.get(packet["drone_id"], -float("inf")):
|
||||
continue
|
||||
self.remote_drone_estimates[packet["drone_id"]] = packet["position_ned_m"]
|
||||
self.remote_drone_velocities[packet["drone_id"]] = packet.get("velocity_ned_m_s", np.zeros(3))
|
||||
self.remote_drone_estimate_times[packet["drone_id"]] = packet["measurement_time_s"]
|
||||
self.remote_drone_receive_times[packet["drone_id"]] = self.time_s
|
||||
self.remote_drone_link_latencies[packet["drone_id"]] = packet["wireless_latency_s"]
|
||||
self.remote_drone_telemetry_age_s[packet["drone_id"]] = self.time_s - packet["measurement_time_s"]
|
||||
self.pending_drone_telemetry = pending
|
||||
|
||||
def _send_assignment_command(self, threat_id, drone_id):
|
||||
record = self.defense_assignments[threat_id]
|
||||
record["attempts"][drone_id] += 1
|
||||
record["last_send"][drone_id] = self.time_s
|
||||
computed = self.time_s+self.drone_link_config.command_compute_latency_s
|
||||
start, sent = self._wireless_transmit(computed, self.drone_link_config.command_payload_bytes)
|
||||
if self.rng.random() < self.drone_link_config.packet_loss_probability:
|
||||
if record["status"][drone_id] != "active":
|
||||
record["status"][drone_id] = "dropped"
|
||||
self.events.append({"type": "drone_command_dropped", "time_s": self.time_s,
|
||||
"drone_id": drone_id, "threat_id": threat_id})
|
||||
return
|
||||
latency = max(0., self.drone_link_config.one_way_latency_s + self.rng.normal(0., self.drone_link_config.latency_jitter_sigma_s))
|
||||
if record["status"][drone_id] != "active":
|
||||
record["status"][drone_id] = "sent"
|
||||
self.pending_drone_commands.append({"drone_id": drone_id, "threat_id": threat_id,
|
||||
"target_position_ned_m": record["waypoints"][drone_id].copy(),
|
||||
"flight_reference": deepcopy(record.get("references", {}).get(drone_id)),
|
||||
"sent_time_s": sent, "computed_time_s": computed,
|
||||
"queue_delay_s": max(0., start-computed), "airtime_s": sent-start,
|
||||
"delivery_time_s": sent+latency})
|
||||
|
||||
def _service_assignments(self):
|
||||
active_threats = {p.entity_id for p in self.projectiles if p.active}
|
||||
for threat_id, record in list(self.defense_assignments.items()):
|
||||
if threat_id not in active_threats:
|
||||
continue
|
||||
failed = False
|
||||
for did in record["defender_ids"]:
|
||||
if record["status"][did] == "active":
|
||||
continue
|
||||
# A queued packet may have long modeled serialization latency;
|
||||
# do not retransmit it while it is still awaiting delivery.
|
||||
if any(c["drone_id"] == did and c["threat_id"] == threat_id for c in self.pending_drone_commands):
|
||||
continue
|
||||
if self.time_s-record["last_send"][did] < 1.:
|
||||
continue
|
||||
if record["attempts"][did] >= 3:
|
||||
failed = True
|
||||
break
|
||||
self._send_assignment_command(threat_id, did)
|
||||
if failed:
|
||||
self._release_assignment(threat_id, "command_delivery_failed")
|
||||
continue
|
||||
track = self.radar_tracks.get(threat_id)
|
||||
if track is not None and track.velocity_observable and self.time_s-track.timestamp_s <= self.radar_config.full_azimuth_revisit_s*1.1 and self.time_s-record["last_replan_s"] >= 1.:
|
||||
velocity = track.velocity_ned_m_s.copy()
|
||||
velocity[2] = 0.
|
||||
norm = float(np.linalg.norm(velocity))
|
||||
if norm > 1e-6:
|
||||
for did, ref in record.get("references", {}).items():
|
||||
ref["direction"] = velocity/norm
|
||||
ref["target_altitude_m"] = max(0., -(track.position_ned_m[2]+track.velocity_ned_m_s[2]*(self.time_s-track.timestamp_s)))
|
||||
if record["status"][did] == "active" and not any(c["drone_id"] == did for c in self.pending_drone_commands):
|
||||
self._send_assignment_command(threat_id, did)
|
||||
record["last_replan_s"] = self.time_s
|
||||
|
||||
def _release_assignment(self, threat_id, reason):
|
||||
record = self.defense_assignments.pop(threat_id)
|
||||
self.pending_drone_commands = [c for c in self.pending_drone_commands if c["threat_id"] != threat_id]
|
||||
for did in record["defender_ids"]:
|
||||
drone = next(d for d in self.drones if d.entity_id == did)
|
||||
drone.activated_for = None
|
||||
drone.flight_reference = None
|
||||
drone.standby_altitude_m = None
|
||||
drone.available = drone.active
|
||||
drone.target_position_ned_m = drone.position_ned_m.copy() if drone.active and drone.position_ned_m[2] < 0 else None
|
||||
self.events.append({"type": "defenders_released", "time_s": self.time_s,
|
||||
"threat_id": threat_id, "defender_ids": record["defender_ids"], "reason": reason})
|
||||
|
||||
def _release_finished_assignments(self):
|
||||
still_active = {p.entity_id for p in self.projectiles if p.active}
|
||||
for threat_id in list(self.defense_assignments):
|
||||
if threat_id not in still_active:
|
||||
self._release_assignment(threat_id, "threat_terminated")
|
||||
|
||||
def _activate_defenders(self, tracks):
|
||||
"""Reserve exactly three nearest defenders after sufficient evidence.
|
||||
|
||||
Uses timestamped estimates and a declared constant-velocity forecast.
|
||||
Reachability is a necessary kinematic screen, not a flight guarantee.
|
||||
Commands, including replans, take effect only after link delivery.
|
||||
"""
|
||||
if self.city_circle is None:
|
||||
return
|
||||
self._release_finished_assignments()
|
||||
for p in self.projectiles:
|
||||
if not p.active or p.entity_id in self.activated_threats:
|
||||
continue
|
||||
if self.time_s - self._last_activation_check_s.get(p.entity_id, -1e9) < 1.0:
|
||||
continue
|
||||
self._last_activation_check_s[p.entity_id] = self.time_s
|
||||
track = tracks.get(p.entity_id)
|
||||
if track is None or not track.velocity_observable:
|
||||
continue
|
||||
if self.time_s - track.timestamp_s > self.radar_config.full_azimuth_revisit_s * 1.1:
|
||||
continue
|
||||
track_time = self.time_s
|
||||
track_velocity = track.velocity_ned_m_s.copy()
|
||||
track_position = track.position_ned_m + track_velocity * (self.time_s-track.timestamp_s)
|
||||
prediction = predict_track(track_position, track_velocity, self.city_circle,
|
||||
self.time_s, ceiling_m=3000.)
|
||||
if prediction is not None:
|
||||
self.track_predictions[p.entity_id] = prediction.samples
|
||||
if prediction is None or prediction.crossing_ned_m is None:
|
||||
continue
|
||||
aim = prediction.crossing_ned_m
|
||||
max_age_s = 2.0 / self.sensor_config.update_rate_hz + self.drone_link_config.one_way_latency_s + 0.05
|
||||
remote_positions = {drone_id: position+self.remote_drone_velocities.get(drone_id, np.zeros(3))*(self.time_s-self.remote_drone_estimate_times[drone_id])
|
||||
for drone_id, position in self.remote_drone_estimates.items()
|
||||
if self.time_s - self.remote_drone_estimate_times.get(drone_id, -1e9) <= max_age_s}
|
||||
if len(remote_positions) < 3:
|
||||
continue
|
||||
chosen = select_defenders(self.drones, aim, count=3, position_estimates=remote_positions,
|
||||
path_segment=(track_position, aim))
|
||||
if len(chosen) < 3:
|
||||
continue
|
||||
# Finite future path: current estimate through first protected entry.
|
||||
delta = aim-track_position
|
||||
length2 = float(delta @ delta)
|
||||
waypoints = {}
|
||||
references = {}
|
||||
direction = track_velocity.copy()
|
||||
direction[2] = 0.
|
||||
direction /= max(float(np.linalg.norm(direction)), 1e-9)
|
||||
for drone in chosen:
|
||||
position = remote_positions[drone.entity_id]
|
||||
fraction = float(np.clip((position-track_position) @ delta/max(length2, 1e-9), 0, 1))
|
||||
waypoint = track_position+fraction*delta
|
||||
waypoints[drone.entity_id] = waypoint
|
||||
scenario = getattr(self, "scenario_config", {})
|
||||
curve = scenario.get("ascent_curve", "quintic")
|
||||
speed = min(drone.config.max_speed_m_s*.8, float(np.linalg.norm(track_velocity[:2]))*.5)
|
||||
length = scenario.get("ascent_distance_m")
|
||||
if length is None:
|
||||
length = automatic_ascent_length(-waypoint[2]+position[2], speed, drone.config.max_climb_rate_m_s, curve)
|
||||
references[drone.entity_id] = {"anchor": waypoint.copy(), "direction": direction.copy(),
|
||||
"target_altitude_m": max(0., -waypoint[2]), "speed_m_s": speed,
|
||||
"ascent_curve": curve, "ascent_distance_m": length,
|
||||
"ascent_k_per_m": scenario.get("ascent_k_per_m", .002)}
|
||||
# Feasibility gate (F05): each defender must be able to reach its
|
||||
# waypoint before the threat is predicted to pass it, using only
|
||||
# estimates + declared vehicle limits. Infeasible assignments are
|
||||
# reported rather than silently committed; the outcome always comes
|
||||
# from simulated motion, never from this test declaring success.
|
||||
crossing_time = prediction.crossing_time_s
|
||||
velocity2 = float(track_velocity @ track_velocity)
|
||||
airtime = 3*(self.drone_link_config.command_payload_bytes+self.drone_link_config.packet_overhead_bytes)*8/self.drone_link_config.bitrate_bps
|
||||
delay = (self.drone_link_config.command_compute_latency_s + self.drone_link_config.one_way_latency_s
|
||||
+ 1/self.drone_link_config.update_rate_hz + airtime
|
||||
+ max(0., self.wireless_channel_free_time_s-self.time_s))
|
||||
available_times = {d.entity_id: float((waypoints[d.entity_id]-track_position) @ track_velocity/max(velocity2, 1e-9))-delay
|
||||
for d in chosen}
|
||||
feasible = [d for d in chosen if reachable_within(
|
||||
remote_positions[d.entity_id], waypoints[d.entity_id], available_times[d.entity_id],
|
||||
d.config.max_speed_m_s, d.config.max_climb_rate_m_s, d.config.max_descent_rate_m_s)]
|
||||
if len(feasible) < 3:
|
||||
self.events.append({"type": "defense_assignment_infeasible", "time_s": self.time_s,
|
||||
"threat_id": p.entity_id, "nearest_defender_ids": [d.entity_id for d in chosen],
|
||||
"reason": "reachability_gate_failed",
|
||||
"projected_city_crossing_ned_m": aim.tolist(),
|
||||
"crossing_time_s": crossing_time, "available_time_by_drone_s": available_times})
|
||||
continue
|
||||
for drone in chosen:
|
||||
drone.available = False
|
||||
# Reservation is not a delivered launch command.
|
||||
# Targets become executable only after command delivery.
|
||||
# Seed the IR seeker RNG so the tasked drone can acquire the
|
||||
# threat once inside seeker range.
|
||||
if drone.ir_rng is None:
|
||||
drone.ir_rng = np.random.default_rng(self.rng.integers(0, 2**63))
|
||||
self.defense_assignments[p.entity_id] = {"defender_ids": [d.entity_id for d in chosen],
|
||||
"waypoints": waypoints, "references": references, "last_replan_s": self.time_s, "status": {d.entity_id: "proposed" for d in chosen},
|
||||
"attempts": {d.entity_id: 0 for d in chosen}, "last_send": {d.entity_id: -1e9 for d in chosen}}
|
||||
for drone in chosen:
|
||||
self._send_assignment_command(p.entity_id, drone.entity_id)
|
||||
self.activated_threats.add(p.entity_id)
|
||||
self.events.append({"type": "defenders_assigned", "time_s": self.time_s,
|
||||
"threat_id": p.entity_id, "defender_ids": [d.entity_id for d in chosen],
|
||||
"projected_city_crossing_ned_m": aim.tolist(),
|
||||
"projected_threat_line_ned_m": prediction.line_origin_ned_m.tolist(),
|
||||
"projected_threat_direction_ned": prediction.line_direction_ned.tolist(),
|
||||
"waypoints_ned_m": {did: wp.tolist() for did, wp in waypoints.items()},
|
||||
"ascent_references": {did: {k: ref[k] for k in ("ascent_curve", "ascent_distance_m", "ascent_k_per_m")} for did, ref in references.items()},
|
||||
"ground_track_time_s": track.timestamp_s,
|
||||
"radar_to_control_latency_s": self.time_s - track.timestamp_s,
|
||||
"command_compute_latency_s": self.drone_link_config.command_compute_latency_s,
|
||||
"wireless_update_rate_hz": self.drone_link_config.update_rate_hz,
|
||||
"wireless_one_way_latency_s": self.drone_link_config.one_way_latency_s,
|
||||
"selection_rule": "three nearest available defenders by 3D distance to finite estimated future path, ties by entity id"})
|
||||
|
||||
def _resolve_physical_events(self, previous, previous_velocity):
|
||||
"""Apply earliest events globally; ground wins an exact contact tie.
|
||||
|
||||
Collision participants are consumed. Collision-as-neutralization is
|
||||
an explicit game assumption, without blast/debris modeling.
|
||||
"""
|
||||
candidates = []
|
||||
for p in self.projectiles:
|
||||
if not p.active:
|
||||
continue
|
||||
ground = ground_fraction(previous[p.entity_id], p.position_ned_m)
|
||||
limit = ground if ground is not None else 1.0
|
||||
if ground is not None:
|
||||
candidates.append((ground, 0, p.entity_id, "", "ground_impact"))
|
||||
if self.city_circle is not None and p.entity_id not in self.protected_entries:
|
||||
entry = circle_entry_fraction(previous[p.entity_id], p.position_ned_m,
|
||||
self.city_circle.center_ned_m, self.city_circle.radius_m)
|
||||
if entry is not None and (ground is None or entry < ground):
|
||||
candidates.append((entry, 2, p.entity_id, "", "protected_area_entry"))
|
||||
for d in self.drones:
|
||||
if not d.active:
|
||||
continue
|
||||
r0 = previous[d.entity_id] - previous[p.entity_id]
|
||||
r1 = d.position_ned_m - p.position_ned_m
|
||||
delta = r1 - r0
|
||||
denom = float(delta @ delta)
|
||||
nearest = float(np.clip(-float(r0 @ delta) / denom, 0, limit)) if denom else 0.0
|
||||
distance = float(np.linalg.norm(r0 + nearest * delta))
|
||||
self.closest_approach_m[p.entity_id] = min(self.closest_approach_m[p.entity_id], distance)
|
||||
contact = sphere_entry_fraction(r0, r1,
|
||||
p.config.collision_radius_m + d.config.collision_radius_m)
|
||||
if contact is not None and contact <= limit:
|
||||
candidates.append((contact, 1, p.entity_id, d.entity_id, "contact"))
|
||||
by_projectile = {p.entity_id: p for p in self.projectiles}
|
||||
by_drone = {d.entity_id: d for d in self.drones}
|
||||
endpoints = {e.entity_id: e.position_ned_m.copy() for e in [*self.projectiles, *self.drones]}
|
||||
endvel = {e.entity_id: e.velocity_ned_m_s.copy() for e in [*self.projectiles, *self.drones]}
|
||||
for fraction, _, pid, did, kind in sorted(candidates):
|
||||
p = by_projectile[pid]
|
||||
if not p.active:
|
||||
continue
|
||||
position = previous[pid] + fraction * (endpoints[pid] - previous[pid])
|
||||
event_time = self.time_s - self.dt + fraction * self.dt
|
||||
if kind == "protected_area_entry":
|
||||
self.protected_entries.add(pid)
|
||||
self.events.append({"type": kind, "time_s": event_time,
|
||||
"projectile_id": pid, "position_ned_m": position.tolist()})
|
||||
continue
|
||||
if kind == "contact":
|
||||
d = by_drone[did]
|
||||
if not d.active:
|
||||
continue
|
||||
d.position_ned_m = previous[did] + fraction * (endpoints[did] - previous[did])
|
||||
d.velocity_ned_m_s = previous_velocity[did] + fraction * (endvel[did] - previous_velocity[did])
|
||||
d.active, d.available = False, False
|
||||
d.target_position_ned_m = None
|
||||
kind = "interception" if d.entity_type == "defense_drone" else "collision"
|
||||
outcome = "intercepted" if kind == "interception" else "collision"
|
||||
event = {"type": kind, "time_s": event_time, "projectile_id": pid,
|
||||
"target_id": did, "position_ned_m": position.tolist(),
|
||||
"miss_distance_m": float(np.linalg.norm(d.position_ned_m - position)),
|
||||
"intercept_error_ned_m": (d.position_ned_m - position).tolist(),
|
||||
"neutralization_assumption": "contact_consumes_both"}
|
||||
else:
|
||||
position[2] = 0.0
|
||||
inside = self.city_circle is not None and float(np.linalg.norm(
|
||||
position[:2] - self.city_circle.center_ned_m[:2])) <= self.city_circle.radius_m
|
||||
outcome = "protected_impact" if inside else "safe_ground_impact"
|
||||
event = {"type": "ground_impact", "time_s": event_time,
|
||||
"projectile_id": pid, "position_ned_m": position.tolist(), "outcome": outcome}
|
||||
p.position_ned_m = position
|
||||
p.velocity_ned_m_s = previous_velocity[pid] + fraction * (endvel[pid] - previous_velocity[pid])
|
||||
p.active = False
|
||||
event["velocity_ned_m_s"] = p.velocity_ned_m_s.tolist()
|
||||
self.projectile_outcomes[pid] = outcome
|
||||
self.events.append(event)
|
||||
|
||||
def snapshot(self) -> dict:
|
||||
entities = []
|
||||
for e in [*self.projectiles, *self.drones]:
|
||||
kind = "projectile" if isinstance(e, Projectile) else e.entity_type
|
||||
row = {"id": e.entity_id, "type": kind, "active": e.active,
|
||||
"position_ned_m": e.position_ned_m.tolist(), "velocity_ned_m_s": e.velocity_ned_m_s.tolist(),
|
||||
"acceleration_ned_m_s2": e.acceleration_ned_m_s2.tolist()}
|
||||
if isinstance(e, Drone):
|
||||
row.update({"orientation_quaternion_wxyz": e.orientation_quaternion_wxyz.tolist(),
|
||||
"angular_velocity_body_rad_s": e.angular_velocity_body_rad_s.tolist(),
|
||||
"position_error_ned_m": (e.position_ned_m - e.estimated_position_ned_m).tolist() if e.estimated_position_ned_m is not None else None,
|
||||
"velocity_error_ned_m_s": (e.velocity_ned_m_s - e.estimated_velocity_ned_m_s).tolist() if e.estimated_velocity_ned_m_s is not None else None,
|
||||
"actuator_acceleration_ned_m_s2": e.actuator_acceleration_ned_m_s2.tolist()})
|
||||
row.update({"available": e.available, "activated_for": e.activated_for,
|
||||
"standby_altitude_m": e.standby_altitude_m,
|
||||
"flight_phase": ("inactive" if not e.active else "assigned" if e.activated_for is not None
|
||||
else ("standby_climb" if -e.position_ned_m[2] < e.standby_altitude_m-2 else "standby_hold")
|
||||
if e.standby_altitude_m is not None else "takeoff_pending"
|
||||
if self.standby_launches.get(e.entity_id, {}).get("status") == "sent" else "idle"),
|
||||
"intended_trajectory_ned_m": e.intended_trajectory(),
|
||||
"intention_kind": "delivered_command_reference",
|
||||
"flight_dynamics": e.flight_dynamics.model_name if e.flight_dynamics is not None else "reduced_kinematic_surrogate",
|
||||
"ir_acquired": e.ir_acquired, "ir_observation": e.ir_observation})
|
||||
elif e.estimated_position_ned_m is not None:
|
||||
row.update({"estimated_position_ned_m": e.estimated_position_ned_m.tolist(),
|
||||
"estimated_velocity_ned_m_s": e.estimated_velocity_ned_m_s.tolist(),
|
||||
"position_error_ned_m": (e.position_ned_m - e.estimated_position_ned_m).tolist(),
|
||||
"velocity_error_ned_m_s": (e.velocity_ned_m_s - e.estimated_velocity_ned_m_s).tolist()})
|
||||
if isinstance(e, Projectile) and e.guidance is not None:
|
||||
row.update({"strike_target_ned_m": e.guidance.target_ned_m.tolist(),
|
||||
"guidance_active": e.guidance_active,
|
||||
"terminal_guidance": {"max_thrust_m_s2": e.guidance.max_thrust_m_s2,
|
||||
"command_rate_hz": e.guidance.command_rate_hz}})
|
||||
row["trajectory_ned_m"] = np.asarray(self.trajectory[e.entity_id][-1800:]).tolist()
|
||||
entities.append(row)
|
||||
if entities:
|
||||
reference_altitude = -float(np.mean([e["position_ned_m"][2] for e in entities]))
|
||||
else:
|
||||
reference_altitude = 0.0
|
||||
air = isa(reference_altitude)
|
||||
observations = {key: {"timestamp_s": value.timestamp_s,
|
||||
"position_ned_m": value.position_ned_m.tolist(),
|
||||
"velocity_ned_m_s": value.velocity_ned_m_s.tolist(),
|
||||
"position_error_ned_m": value.position_error_ned_m.tolist(),
|
||||
"velocity_error_ned_m_s": value.velocity_error_ned_m_s.tolist(),
|
||||
"error_reference": value.error_reference,
|
||||
**({"control_station_receive_time_s": self.remote_drone_receive_times[key],
|
||||
"wireless_one_way_latency_s": self.remote_drone_link_latencies[key],
|
||||
"remote_measurement_time_s": self.remote_drone_estimate_times[key],
|
||||
"remote_position_ned_m": self.remote_drone_estimates[key].tolist(),
|
||||
"data_age_at_control_station_s": self.time_s - self.remote_drone_estimate_times[key]}
|
||||
if key in self.remote_drone_receive_times else {})}
|
||||
for key, value in self.latest_observations.items()}
|
||||
closest = min(self.closest_approach_m.values(), default=float("inf"))
|
||||
closest = closest if math.isfinite(closest) else None
|
||||
return {"schema_version": 1, "run_id": self.run_id, "sequence": int(round(self.time_s / self.dt)),
|
||||
"simulation_time_s": self.time_s, "kind": "simulation_state", "entities": entities,
|
||||
"events": list(self.events[-250:]), "environment": {"wind_ned_m_s": self.environment.wind_ned_m_s.tolist(),
|
||||
"gust_ned_m_s": self.environment.gust_ned_m_s.tolist(),
|
||||
"altitude_m": reference_altitude, "temperature_k": air.temperature_k, "pressure_pa": air.pressure_pa,
|
||||
"density_kg_m3": air.density_kg_m3, "speed_of_sound_m_s": air.speed_of_sound_m_s,
|
||||
"viscosity_pa_s": air.viscosity_pa_s,
|
||||
"city_center_ned_m": self.city_circle.center_ned_m.tolist() if self.city_circle else [0.0, 0.0, 0.0],
|
||||
"city_radius_m": self.city_circle.radius_m if self.city_circle else 0.0,
|
||||
"protected_radius_m": self.city_circle.radius_m if self.city_circle else 0.0,
|
||||
"radar_radius_m": self.radar_config.ring_radius_km * 1000.0,
|
||||
"drone_radius_m": (getattr(self, "scenario_config", {}).get("drone_radius_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_perimeter_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_ring_km", 0.0)) * 1000.0,
|
||||
"defense_perimeter_radius_m": (getattr(self, "scenario_config", {}).get("drone_radius_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_perimeter_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_ring_km", 0.0)) * 1000.0,
|
||||
"defense_ring_radius_m": (getattr(self, "scenario_config", {}).get("drone_radius_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_perimeter_km", 0.0)
|
||||
or getattr(self, "scenario_config", {}).get("defense_ring_km", 0.0)) * 1000.0},
|
||||
"devices": {"radars": [{"id": r.radar_id, "position_ned_m": r.position_ned_m.tolist(),
|
||||
"azimuth_slot": r.azimuth_slot, "last_detection_time_s": r.last_detection_time_s}
|
||||
for r in self.radars],
|
||||
"radar_profile": {"model": self.radar_config.model, "max_range_m": self.radar_config.max_range_m,
|
||||
"azimuth_step_period_s": self.radar_config.azimuth_step_period_s,
|
||||
"full_azimuth_revisit_s": self.radar_config.full_azimuth_revisit_s,
|
||||
"wired_processing_latency_s": self.radar_config.wired_processing_latency_s},
|
||||
"wireless_link_profile": {"update_rate_hz": self.drone_link_config.update_rate_hz,
|
||||
"one_way_latency_s": self.drone_link_config.one_way_latency_s,
|
||||
"latency_jitter_sigma_s": self.drone_link_config.latency_jitter_sigma_s,
|
||||
"command_compute_latency_s": self.drone_link_config.command_compute_latency_s,
|
||||
"bitrate_bps": self.drone_link_config.bitrate_bps,
|
||||
"packet_loss_probability": self.drone_link_config.packet_loss_probability,
|
||||
"channel_airtime_s": self.wireless_airtime_s},
|
||||
"drone_control_rate_hz": self.drone_control_rate_hz},
|
||||
"scenario": getattr(self, "scenario_config", {}), "observations": observations,
|
||||
"radar_tracks": {key: {"timestamp_s": value.timestamp_s,
|
||||
"position_ned_m": value.position_ned_m.tolist(),
|
||||
"velocity_ned_m_s": value.velocity_ned_m_s.tolist(),
|
||||
"covariance": value.covariance.tolist() if value.covariance is not None else None,
|
||||
"sample_count": value.sample_count, "velocity_observable": value.velocity_observable,
|
||||
"age_s": self.time_s - value.timestamp_s}
|
||||
for key, value in self.radar_tracks.items()},
|
||||
"radar_detection_points": self.radar_detection_points,
|
||||
"inferred_trajectories": {key: value[:, 1:].tolist() for key, value in self.track_predictions.items()},
|
||||
"latest_detections": {key: {"radar_id": value["radar_id"],
|
||||
"measurement_time_s": value["measurement_time_s"],
|
||||
"position_ned_m": value["position_ned_m"].tolist(),
|
||||
"line_of_sight_ned": value["line_of_sight_ned"].tolist(),
|
||||
"radial_velocity_m_s": value["radial_velocity_m_s"]}
|
||||
for key, value in self.latest_detections.items()},
|
||||
"statistics": {"elapsed_time_s": self.time_s,
|
||||
"projectile_outcomes": {p.entity_id: self.projectile_outcomes.get(p.entity_id, "unresolved") for p in self.projectiles},
|
||||
"neutralization_assumption": "contact_consumes_both",
|
||||
"protected_area_entry_count": len(self.protected_entries),
|
||||
"interception_count": sum(e["type"] == "interception" for e in self.events),
|
||||
"closest_approach_m": closest}}
|
||||
|
||||
def run(self, steps: int) -> list[dict]:
|
||||
return [self.step_and_snapshot() for _ in range(steps)]
|
||||
@@ -1,271 +0,0 @@
|
||||
import json
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from defense_sim.entities.drone import Drone, DroneConfig
|
||||
from defense_sim.entities.projectile import Projectile, ProjectileConfig
|
||||
from defense_sim.physics.events import sphere_entry_fraction, circle_entry_fraction
|
||||
from defense_sim.physics.reference import integrate_reference
|
||||
from defense_sim.sensors.tracking import RadarTracker
|
||||
from defense_sim.simulation.engine import SimulationEngine
|
||||
from defense_sim.statistics.monte_carlo import run_many
|
||||
from defense_sim.sensors.models import Observation, SensorConfig
|
||||
from defense_sim.guidance.activation import CityCircle
|
||||
|
||||
|
||||
def projectile(position, velocity, radius=0.1, name="p"):
|
||||
return Projectile(name, np.array(position, dtype=float), np.array(velocity, dtype=float),
|
||||
ProjectileConfig(area_m2=0, collision_radius_m=radius))
|
||||
|
||||
|
||||
def drone(position, name="d", kind="defense_drone"):
|
||||
return Drone(name, kind, np.array(position, dtype=float), config=DroneConfig(collision_radius_m=0.1))
|
||||
|
||||
|
||||
def test_contact_root_and_shallow_circle_chord():
|
||||
assert sphere_entry_fraction([-1, 0, 0], [1, 0, 0], 0.2) == pytest.approx(0.4)
|
||||
assert circle_entry_fraction([-2, .9, -1], [2, .9, -1], [0, 0, 0], 1) is not None
|
||||
|
||||
|
||||
def test_ground_precedes_underground_collision():
|
||||
p = projectile([0, 0, -1], [0, 0, 20], radius=.01)
|
||||
sim = SimulationEngine(dt=.1, drone_control_rate_hz=10, seed=1, projectiles=[p], drones=[drone([0, 0, 1])])
|
||||
sim.step()
|
||||
assert [e["type"] for e in sim.events] == ["ground_impact"]
|
||||
assert p.position_ned_m[2] == 0
|
||||
assert p.velocity_ned_m_s[2] == pytest.approx(sim.events[0]["velocity_ned_m_s"][2])
|
||||
|
||||
|
||||
@pytest.mark.parametrize("reverse", [False, True])
|
||||
def test_first_contact_is_order_independent_and_consumes_drone(reverse):
|
||||
p = projectile([-1, 0, -10], [200, 0, 0])
|
||||
drones = [drone([.5, 0, -10], "later"), drone([0, 0, -10], "earlier")]
|
||||
if reverse:
|
||||
drones.reverse()
|
||||
sim = SimulationEngine(dt=.01, seed=1, projectiles=[p], drones=drones)
|
||||
sim.step()
|
||||
event = sim.events[0]
|
||||
assert event["type"] == "interception" and event["target_id"] == "earlier"
|
||||
assert event["time_s"] == pytest.approx(.004, abs=1e-6)
|
||||
assert p.position_ned_m[0] == pytest.approx(-.2, abs=1e-5)
|
||||
assert not next(d for d in drones if d.entity_id == "earlier").active
|
||||
assert sim.trajectory["p"][-1][0] == pytest.approx(p.position_ned_m[0])
|
||||
|
||||
|
||||
def test_idle_ground_support_under_wind():
|
||||
d = drone([1, 2, 0])
|
||||
sim = SimulationEngine(seed=1, drones=[d])
|
||||
sim.environment.wind_ned_m_s = np.array([20., 3., 0.])
|
||||
for _ in range(500):
|
||||
sim.step()
|
||||
assert np.array_equal(d.position_ned_m, [1, 2, 0])
|
||||
assert np.array_equal(d.velocity_ned_m_s, np.zeros(3))
|
||||
|
||||
|
||||
def test_scipy_reference_matches_vacuum_ground_time_and_rk4():
|
||||
p = projectile([0, 0, -10], [20, 0, 0])
|
||||
state = np.r_[p.position_ned_m, p.velocity_ned_m_s]
|
||||
reference = integrate_reference(lambda y: p.derivative(y, np.zeros(3)), state, 2, ground_event=True)
|
||||
ground_time = np.sqrt(20/9.80665)
|
||||
assert reference.t_events[0][0] == pytest.approx(ground_time, abs=1e-9)
|
||||
p.step(.1, np.zeros(3))
|
||||
assert np.allclose(np.r_[p.position_ned_m, p.velocity_ned_m_s], reference.sol(.1), atol=1e-9)
|
||||
|
||||
|
||||
def measurement(t, radar="r"):
|
||||
position = np.array([1000.+10*t, 0, -10.])
|
||||
return {"target_id": "p", "radar_id": radar, "measurement_time_s": t,
|
||||
"position_ned_m": position, "range_m": np.linalg.norm(position),
|
||||
"azimuth_rad": 0., "elevation_rad": np.arctan2(10., position[0]),
|
||||
"range_sigma_m": 1., "azimuth_sigma_rad": .001, "elevation_sigma_rad": .001}
|
||||
|
||||
|
||||
def test_tracking_timestamp_covariance_and_velocity_evidence():
|
||||
tracker = RadarTracker(acceleration_diffusion=.01)
|
||||
first = tracker.update(measurement(0))
|
||||
assert not first.velocity_observable
|
||||
for t in (1., 2., 3.):
|
||||
result = tracker.update(measurement(t))
|
||||
assert result.velocity_observable
|
||||
assert result.velocity[0] == pytest.approx(10., abs=.1)
|
||||
assert np.linalg.eigvalsh(result.covariance).min() > 0
|
||||
assert tracker.update(measurement(3)) is None
|
||||
assert tracker.update(measurement(1, "other")) is None
|
||||
|
||||
|
||||
def test_mature_tracks_and_event_covariance_serialize_as_strict_json():
|
||||
sim = SimulationEngine(seed=1, projectiles=[projectile([1000, 0, -10], [10, 0, 0])])
|
||||
for t in (0., 1., 2., 3.):
|
||||
sim.time_s = t
|
||||
sim.radar_arrivals.append({**measurement(t), 'arrival_time_s': t})
|
||||
sim._deliver_radar_tracks()
|
||||
assert sim.radar_tracks['p'].velocity_observable
|
||||
json.dumps(sim.snapshot(), allow_nan=False)
|
||||
|
||||
|
||||
def test_monte_carlo_outcomes_are_per_projectile_and_cutoffs_censored():
|
||||
def factory(seed):
|
||||
return SimulationEngine(seed=seed, projectiles=[
|
||||
projectile([-1, 0, -10], [200, 0, 0], name="hit"),
|
||||
projectile([100, 0, -100], [0, 0, 0], name="unfinished")],
|
||||
drones=[drone([0, 0, -10])])
|
||||
result = run_many(factory, runs=2, steps=1)
|
||||
groups = result.aggregate["interception_by_projectile"]
|
||||
assert groups["hit"]["interception_probability"] == 1
|
||||
assert groups["unfinished"]["interception_probability"] is None
|
||||
assert groups["unfinished"]["unresolved_count"] == 2
|
||||
json.dumps(result.aggregate, allow_nan=False)
|
||||
|
||||
|
||||
def launch_fixture(loss=0):
|
||||
from defense_sim.sensors.devices import DroneLinkConfig
|
||||
p = projectile([1000, 0, -20], [-50, 0, 0])
|
||||
from defense_sim.entities.projectile import ProjectileGuidanceConfig
|
||||
p.guidance = ProjectileGuidanceConfig(np.array([0., 0., 0.]), max_thrust_m_s2=20.)
|
||||
p.guidance_active = True
|
||||
drones = [drone([500, i*5, 0], f'd{i}') for i in range(4)]
|
||||
sim = SimulationEngine(seed=3, projectiles=[p], drones=drones,
|
||||
city_circle=CityCircle(np.zeros(3), 500),
|
||||
sensor_config=SensorConfig(position_sigma_m=0, velocity_sigma_m_s=0),
|
||||
drone_link_config=DroneLinkConfig(packet_loss_probability=loss))
|
||||
track = Observation('p', 0., p.position_ned_m.copy(), p.velocity_ned_m_s.copy(),
|
||||
np.zeros(3), np.zeros(3), np.eye(6), 5, True)
|
||||
sim.radars = [] # this fixture isolates explicitly injected tracking evidence
|
||||
sim.radar_tracks['p'] = track
|
||||
sim.remote_drone_estimates = {d.entity_id: d.position_ned_m.copy() for d in drones}
|
||||
sim.remote_drone_estimate_times = {d.entity_id: 0. for d in drones}
|
||||
return sim
|
||||
|
||||
|
||||
def test_launch_requires_evidence_and_delivered_commands():
|
||||
sim = launch_fixture()
|
||||
sim.radar_tracks['p'].velocity_observable = False
|
||||
sim._activate_defenders(sim.radar_tracks)
|
||||
assert all(d.available for d in sim.drones)
|
||||
sim.radar_tracks['p'].velocity_observable = True
|
||||
sim.time_s = 1.
|
||||
sim.remote_drone_estimate_times = {d.entity_id: 1. for d in sim.drones}
|
||||
sim._activate_defenders(sim.radar_tracks)
|
||||
assert sum(not d.available for d in sim.drones) == 3
|
||||
assert all(d.activated_for is None and d.target_position_ned_m is None for d in sim.drones)
|
||||
sim.time_s = max(c['delivery_time_s'] for c in sim.pending_drone_commands)+.01
|
||||
sim._deliver_drone_commands()
|
||||
assert sum(d.activated_for == 'p' for d in sim.drones) == 3
|
||||
for _ in range(500):
|
||||
sim.step()
|
||||
assert any(d.position_ned_m[2] < -.01 for d in sim.drones)
|
||||
assert all(d.position_ned_m[2] <= 0 for d in sim.drones)
|
||||
|
||||
|
||||
def test_total_packet_loss_retries_and_releases_reservations():
|
||||
sim = launch_fixture(loss=1.)
|
||||
sim._activate_defenders(sim.radar_tracks)
|
||||
for t in (1., 2., 3.):
|
||||
sim.time_s = t
|
||||
sim._service_assignments()
|
||||
assert not sim.defense_assignments
|
||||
assert all(d.available and d.activated_for is None for d in sim.drones)
|
||||
assert len([e for e in sim.events if e['type'] == 'drone_command_dropped']) == 9
|
||||
|
||||
|
||||
def test_ir_bearing_only_processing_delay_and_loss():
|
||||
d = drone([0, 0, -10])
|
||||
d.activated_for = 'p'
|
||||
d.ir_rng = np.random.default_rng(1)
|
||||
d.config.ir_lock_time_s = 0.
|
||||
d.config.ir_processing_latency_s = .05
|
||||
p = projectile([-100, 0, -10], [0, 0, 0])
|
||||
prior = Observation('p', 0., p.position_ned_m.copy(), np.zeros(3), np.zeros(3),
|
||||
np.zeros(3), np.eye(6)*100, 5, True)
|
||||
events = []
|
||||
assert not d.estimate_threat(p, 0., d.ir_rng, events, prior)
|
||||
assert d.estimate_threat(p, .05, d.ir_rng, events, prior)
|
||||
assert d.ir_observation['measurement_type'] == 'bearing_only'
|
||||
assert 'range_m' not in d.ir_observation
|
||||
p.position_ned_m = np.array([100., 0., -10.]) # outside rear-facing FOV
|
||||
assert not d.estimate_threat(p, .1, d.ir_rng, events, prior)
|
||||
assert not d.ir_acquired and d.ir_threat_estimate is None
|
||||
|
||||
|
||||
@pytest.mark.parametrize('dt', [.1, .05, .025])
|
||||
def test_ground_event_converges_to_scipy_reference(dt):
|
||||
p = projectile([0, 0, -10], [0, 0, 0])
|
||||
sim = SimulationEngine(dt=dt, drone_control_rate_hz=min(50, 1/dt), seed=1, projectiles=[p])
|
||||
while p.active:
|
||||
sim.step()
|
||||
analytic = np.sqrt(20/9.80665)
|
||||
assert abs(sim.events[-1]['time_s']-analytic) < dt*dt
|
||||
|
||||
|
||||
def test_atmosphere_high_altitude_and_nonfinite_rejected():
|
||||
from defense_sim.physics.atmosphere import isa
|
||||
assert 1e-5 < isa(80000).density_kg_m3 < 3e-5
|
||||
assert isa(60000).density_kg_m3 < isa(20000).density_kg_m3/100
|
||||
with pytest.raises(ValueError):
|
||||
DroneConfig(max_speed_m_s=float('nan'))
|
||||
with pytest.raises(ValueError):
|
||||
isa(90000)
|
||||
|
||||
|
||||
def test_radar_dwell_evaluates_all_targets_and_preserves_clock_phase():
|
||||
from defense_sim.sensors.devices import Aeris10Radar, AerisRadarConfig
|
||||
config = AerisRadarConfig(range_sigma_m=0, azimuth_sigma_deg=0, elevation_sigma_deg=0, doppler_sigma_m_s=0)
|
||||
visible = projectile([1000, 0, -10], [-10, 0, 0], name='visible')
|
||||
outside = projectile([30000, 0, -10], [-10, 0, 0], name='outside')
|
||||
for targets in ([outside, visible], [visible, outside]):
|
||||
radar = Aeris10Radar('r', np.zeros(3), 0, config, np.random.default_rng(3))
|
||||
assert [m['target_id'] for m in radar.scan_dwell(targets, 0)] == ['visible']
|
||||
for t in np.arange(.01, 13.01, .01):
|
||||
radar.scan_dwell([], float(t))
|
||||
assert radar.next_sweep_time_s == pytest.approx(51*config.azimuth_step_period_s)
|
||||
|
||||
|
||||
def test_telemetry_out_of_order_cannot_overwrite_newer_sample():
|
||||
sim = SimulationEngine(seed=1)
|
||||
sim.time_s = 2.
|
||||
base = {'drone_id': 'd', 'delivery_time_s': 1., 'wireless_latency_s': .1}
|
||||
sim.pending_drone_telemetry = [
|
||||
{**base, 'measurement_time_s': 1., 'position_ned_m': np.ones(3)},
|
||||
{**base, 'measurement_time_s': .5, 'position_ned_m': np.zeros(3)}]
|
||||
sim._deliver_drone_telemetry()
|
||||
assert np.array_equal(sim.remote_drone_estimates['d'], np.ones(3))
|
||||
assert sim.remote_drone_estimate_times['d'] == 1.
|
||||
|
||||
|
||||
def test_guidance_does_not_overwrite_velocity_or_create_zero_thrust():
|
||||
from defense_sim.entities.projectile import ProjectileGuidanceConfig
|
||||
p = projectile([1000, 0, -100], [-50, 0, 0])
|
||||
p.guidance = ProjectileGuidanceConfig(np.zeros(3), max_thrust_m_s2=5.)
|
||||
p.guidance_active = True
|
||||
initial = p.velocity_ned_m_s.copy()
|
||||
p.update_guidance(0., np.zeros(3))
|
||||
assert np.array_equal(p.velocity_ned_m_s, initial)
|
||||
assert np.linalg.norm(p._thrust_cmd) <= 5.+1e-12
|
||||
p.guidance.max_thrust_m_s2 = 0.
|
||||
p.update_guidance(1., np.zeros(3))
|
||||
assert np.array_equal(p._thrust_cmd, np.zeros(3))
|
||||
|
||||
|
||||
def test_ground_wins_tie_and_other_drone_contact_is_not_defense_success():
|
||||
p = projectile([0, 0, 0], [0, 0, 0])
|
||||
sim = SimulationEngine(seed=1, projectiles=[p], drones=[drone([0, 0, 0])])
|
||||
sim.step()
|
||||
assert sim.events[0]['type'] == 'ground_impact'
|
||||
q = projectile([-1, 0, -10], [200, 0, 0])
|
||||
sim = SimulationEngine(seed=1, projectiles=[q], drones=[drone([0, 0, -10], kind='attack_drone')])
|
||||
sim.step()
|
||||
assert sim.events[0]['type'] == 'collision'
|
||||
assert sim.projectile_outcomes['p'] == 'collision'
|
||||
|
||||
|
||||
@pytest.mark.parametrize('reverse', [False, True])
|
||||
def test_consumed_drone_cannot_contact_later_threat(reverse):
|
||||
first = projectile([-.8, 0, -10], [200, 0, 0], name='first')
|
||||
later = projectile([-1., 0, -10], [200, 0, 0], name='later')
|
||||
threats = [later, first]
|
||||
if reverse:
|
||||
threats.reverse()
|
||||
sim = SimulationEngine(seed=1, projectiles=threats, drones=[drone([0, 0, -10])])
|
||||
sim.step()
|
||||
assert not first.active and later.active
|
||||
assert [e['projectile_id'] for e in sim.events if e['type'] == 'interception'] == ['first']
|
||||
@@ -1,101 +0,0 @@
|
||||
import json
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from defense_sim.simulation.scenarios import basic_intercept
|
||||
from test_audit_regressions import launch_fixture, measurement
|
||||
|
||||
|
||||
def detect_once(sim):
|
||||
sim.radar_tracks.clear()
|
||||
sim.radar_arrivals.append({**measurement(sim.time_s), 'arrival_time_s': sim.time_s})
|
||||
sim._deliver_radar_tracks()
|
||||
|
||||
|
||||
def deliver_standby(sim):
|
||||
sim.time_s = max(c['delivery_time_s'] for c in sim.pending_drone_commands)+.01
|
||||
sim._deliver_drone_commands()
|
||||
|
||||
|
||||
def test_first_delivered_detection_launches_before_velocity_or_assignment():
|
||||
sim = launch_fixture()
|
||||
sim.scenario_config = {'standby_altitude_m': 42., 'early_takeoff_enabled': True}
|
||||
detect_once(sim)
|
||||
assert not sim.radar_tracks['p'].velocity_observable
|
||||
assert len(sim.pending_drone_commands) == 4
|
||||
assert all(d.target_position_ned_m is None for d in sim.drones)
|
||||
assert not sim.defense_assignments
|
||||
deliver_standby(sim)
|
||||
assert all(d.standby_altitude_m == 42. and d.available and d.activated_for is None for d in sim.drones)
|
||||
for _ in range(300):
|
||||
sim.step()
|
||||
assert all(d.position_ned_m[2] < -1 for d in sim.drones)
|
||||
snapshot = sim.snapshot()
|
||||
row = next(r for r in snapshot['entities'] if r['id'] == 'd0')
|
||||
assert row['flight_phase'] == 'standby_climb'
|
||||
intention = np.array(row['intended_trajectory_ned_m'])
|
||||
assert intention[-1, 2] == pytest.approx(-42.)
|
||||
assert np.allclose(intention[:, :2], intention[0, :2])
|
||||
assert not np.array_equal(intention, row['trajectory_ned_m'])
|
||||
json.dumps(snapshot, allow_nan=False)
|
||||
|
||||
|
||||
def test_early_takeoff_can_be_disabled_and_altitude_is_validated():
|
||||
sim = launch_fixture()
|
||||
sim.scenario_config = {'early_takeoff_enabled': False}
|
||||
detect_once(sim)
|
||||
assert not sim.pending_drone_commands
|
||||
for altitude in (0, -1, float('nan'), 3001):
|
||||
with pytest.raises(ValueError):
|
||||
basic_intercept(standby_altitude_m=altitude)
|
||||
assert basic_intercept(standby_altitude_m=250.).scenario_config['standby_altitude_m'] == 250.
|
||||
|
||||
|
||||
def test_standby_packets_retry_and_total_loss_does_not_launch():
|
||||
sim = launch_fixture(loss=1.)
|
||||
detect_once(sim)
|
||||
for t in (1., 2., 3.):
|
||||
sim.time_s = t
|
||||
sim._service_standby_launches()
|
||||
assert not sim.pending_drone_commands
|
||||
assert all(d.target_position_ned_m is None for d in sim.drones)
|
||||
assert all(r['status'] == 'failed' and r['attempts'] == 3 for r in sim.standby_launches.values())
|
||||
|
||||
|
||||
def test_late_standby_packet_cannot_override_assignment():
|
||||
sim = launch_fixture()
|
||||
detect_once(sim)
|
||||
d = sim.drones[0]
|
||||
d.activated_for = 'p'
|
||||
d.target_position_ned_m = np.array([123., 456., -80.])
|
||||
deliver_standby(sim)
|
||||
assert d.standby_altitude_m is None
|
||||
assert np.array_equal(d.target_position_ned_m, [123., 456., -80.])
|
||||
|
||||
|
||||
def test_traces_persist_longer_than_old_six_second_window_and_samples_persist():
|
||||
sim = launch_fixture()
|
||||
sim.scenario_config = {'early_takeoff_enabled': False}
|
||||
sim.radar_detection_points['p'] = [[1000., 0., -20.]]
|
||||
for _ in range(1000):
|
||||
sim.step()
|
||||
snapshot = sim.snapshot()
|
||||
path = next(e['trajectory_ned_m'] for e in snapshot['entities'] if e['id'] == 'p')
|
||||
assert path[0] == [1000., 0., -20.]
|
||||
assert len(path) >= 50
|
||||
assert snapshot['radar_detection_points']['p'] == [[1000., 0., -20.]]
|
||||
|
||||
|
||||
def test_intention_preview_does_not_advance_controller_or_truth():
|
||||
sim = launch_fixture()
|
||||
sim._activate_defenders(sim.radar_tracks)
|
||||
deliver_standby(sim) # helper delivers all queued packets, including assignment
|
||||
d = next(d for d in sim.drones if d.flight_reference is not None)
|
||||
position, target = d.position_ned_m.copy(), d.target_position_ned_m.copy()
|
||||
travel, elapsed = d.flight_reference['travel_m'], d.flight_reference['elapsed_s']
|
||||
path = np.array(d.intended_trajectory())
|
||||
assert len(path) == 61 and np.isfinite(path).all()
|
||||
assert np.linalg.norm(path[-1, :2]-path[0, :2]) > 0
|
||||
assert np.array_equal(position, d.position_ned_m)
|
||||
assert np.array_equal(target, d.target_position_ned_m)
|
||||
assert (travel, elapsed) == (d.flight_reference['travel_m'], d.flight_reference['elapsed_s'])
|
||||
@@ -5,7 +5,7 @@
|
||||
- Date: 2026-03-21
|
||||
- Target: Trenz `TE0713-03-82C46-A` on `TE0701-06`
|
||||
- Top module: `radar_system_top_te0713_dev`
|
||||
- Constraint file: `9_Firmware/9_2_FPGA/constraints/te0713_te0701_minimal.xdc`
|
||||
- Constraint file: `firmware/9_2_FPGA/constraints/te0713_te0701_minimal.xdc`
|
||||
- Bitstream: `docs/artifacts/te0713-te0701-heartbeat-2026-03-21.bit`
|
||||
|
||||
Build result:
|
||||
|
||||
@@ -72,7 +72,7 @@ Previous builds failed setup timing due to source-synchronous clock skew:
|
||||
|
||||
## Source Files
|
||||
|
||||
- `9_Firmware/9_2_FPGA/radar_system_top_te0713_umft601x_dev.v` — wrapper top
|
||||
- `9_Firmware/9_2_FPGA/usb_data_interface.v` — FT601 USB data interface
|
||||
- `9_Firmware/9_2_FPGA/constraints/te0713_te0701_umft601x.xdc` — FMC LPC constraints
|
||||
- `9_Firmware/9_2_FPGA/scripts/build_te0713_umft601x_dev.tcl` — build script
|
||||
- `firmware/9_2_FPGA/radar_system_top_te0713_umft601x_dev.v` — wrapper top
|
||||
- `firmware/9_2_FPGA/usb_data_interface.v` — FT601 USB data interface
|
||||
- `firmware/9_2_FPGA/constraints/te0713_te0701_umft601x.xdc` — FMC LPC constraints
|
||||
- `firmware/9_2_FPGA/scripts/build_te0713_umft601x_dev.tcl` — build script
|
||||
|
||||
@@ -3,7 +3,7 @@ Copyright 1986-2022 Xilinx, Inc. All Rights Reserved. Copyright 2022-2025 Advanc
|
||||
| Tool Version : Vivado v.2025.2 (lin64) Build 6299465 Fri Nov 14 12:34:56 MST 2025
|
||||
| Date : Sat Mar 21 20:41:16 2026
|
||||
| Host : jason-pc running 64-bit Ubuntu 24.04.4 LTS
|
||||
| Command : report_timing_summary -report_unconstrained -max_paths 100 -file /home/jason-stone/PLFM_RADAR_work/PLFM_RADAR/9_Firmware/9_2_FPGA/vivado_te0713_umft601x_dev/reports/timing_summary.rpt
|
||||
| Command : report_timing_summary -report_unconstrained -max_paths 100 -file /home/jason-stone/PLFM_RADAR_work/PLFM_RADAR/firmware/9_2_FPGA/vivado_te0713_umft601x_dev/reports/timing_summary.rpt
|
||||
| Design : radar_system_top_te0713_umft601x_dev
|
||||
| Device : 7a200t-fbg484
|
||||
| Speed File : -2 PRODUCTION 1.23 2018-06-13
|
||||
|
||||
@@ -493,12 +493,12 @@
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr><td>Production-target XDC</td><td><code>9_Firmware/9_2_FPGA/constraints/xc7a200t_fbg484.xdc</code></td><td>Constraint source of truth for the production FPGA target</td><td>Tracked and validated after Build 16 cleanup port</td></tr>
|
||||
<tr><td>FPGA programming flow</td><td><code>9_Firmware/9_2_FPGA/scripts/program_fpga.tcl</code></td><td>Programs the device and reports DONE / INIT_COMPLETE / probes presence</td><td>Primary operator-facing programming script</td></tr>
|
||||
<tr><td>Debug probe insertion flow</td><td><code>9_Firmware/9_2_FPGA/scripts/insert_ila_probes.tcl</code></td><td>Used when generating or refreshing debug-capable images</td><td>Keep matched with the selected debug bitstream</td></tr>
|
||||
<tr><td>FPGA regression runner</td><td><code>9_Firmware/9_2_FPGA/run_regression.sh</code></td><td>Pre-arrival regression evidence for the tracked FPGA baseline</td><td>23 / 23 passing on the current tracked branch (includes CFAR + MTI + E2E tests)</td></tr>
|
||||
<tr><td>MCU regression harness</td><td><code>9_Firmware/9_1_Microcontroller/tests/Makefile</code></td><td>Pre-arrival firmware regression evidence before flashing hardware</td><td>20 / 20 passing on the current tracked branch</td></tr>
|
||||
<tr><td>Bring-up logging macros</td><td><code>9_Firmware/9_1_Microcontroller/9_1_1_C_Cpp_Libraries/diag_log.h</code></td><td>Defines the main first-power-on log vocabulary used over USART3</td><td>Observation-only instrumentation layer</td></tr>
|
||||
<tr><td>Production-target XDC</td><td><code>firmware/9_2_FPGA/constraints/xc7a200t_fbg484.xdc</code></td><td>Constraint source of truth for the production FPGA target</td><td>Tracked and validated after Build 16 cleanup port</td></tr>
|
||||
<tr><td>FPGA programming flow</td><td><code>firmware/9_2_FPGA/scripts/program_fpga.tcl</code></td><td>Programs the device and reports DONE / INIT_COMPLETE / probes presence</td><td>Primary operator-facing programming script</td></tr>
|
||||
<tr><td>Debug probe insertion flow</td><td><code>firmware/9_2_FPGA/scripts/insert_ila_probes.tcl</code></td><td>Used when generating or refreshing debug-capable images</td><td>Keep matched with the selected debug bitstream</td></tr>
|
||||
<tr><td>FPGA regression runner</td><td><code>firmware/9_2_FPGA/run_regression.sh</code></td><td>Pre-arrival regression evidence for the tracked FPGA baseline</td><td>23 / 23 passing on the current tracked branch (includes CFAR + MTI + E2E tests)</td></tr>
|
||||
<tr><td>MCU regression harness</td><td><code>firmware/9_1_Microcontroller/tests/Makefile</code></td><td>Pre-arrival firmware regression evidence before flashing hardware</td><td>20 / 20 passing on the current tracked branch</td></tr>
|
||||
<tr><td>Bring-up logging macros</td><td><code>firmware/9_1_Microcontroller/9_1_1_C_Cpp_Libraries/diag_log.h</code></td><td>Defines the main first-power-on log vocabulary used over USART3</td><td>Observation-only instrumentation layer</td></tr>
|
||||
<tr><td>Board-day worksheet</td><td><code>docs/board-day-worksheet.html</code></td><td>Record pass/fail, measurements, and blockers during first sessions</td><td>Use with this page and the bring-up plan</td></tr>
|
||||
<tr><td>Bring-up execution plan</td><td><code>docs/bring-up.html</code></td><td>Operator checklist, abort criteria, observability targets, and open risks</td><td>Primary readiness document</td></tr>
|
||||
</tbody>
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 195 KiB |
@@ -10,15 +10,15 @@ extern SPI_HandleTypeDef hspi1;
|
||||
extern UART_HandleTypeDef huart3;
|
||||
|
||||
// Chip Select GPIO definitions
|
||||
static const struct {
|
||||
GPIO_TypeDef* port;
|
||||
uint16_t pin;
|
||||
} CHIP_SELECTS[4] = {
|
||||
{ADAR_1_CS_3V3_GPIO_Port, ADAR_1_CS_3V3_Pin}, // ADAR1000 #1
|
||||
{ADAR_2_CS_3V3_GPIO_Port, ADAR_2_CS_3V3_Pin}, // ADAR1000 #2
|
||||
{ADAR_3_CS_3V3_GPIO_Port, ADAR_3_CS_3V3_Pin}, // ADAR1000 #3
|
||||
{ADAR_4_CS_3V3_GPIO_Port, ADAR_4_CS_3V3_Pin} // ADAR1000 #4
|
||||
};
|
||||
static const struct {
|
||||
GPIO_TypeDef* port;
|
||||
uint16_t pin;
|
||||
} CHIP_SELECTS[4] = {
|
||||
{ADAR_1_CS_3V3_GPIO_Port, ADAR_1_CS_3V3_Pin}, // ADAR1000 #1
|
||||
{ADAR_2_CS_3V3_GPIO_Port, ADAR_2_CS_3V3_Pin}, // ADAR1000 #2
|
||||
{ADAR_3_CS_3V3_GPIO_Port, ADAR_3_CS_3V3_Pin}, // ADAR1000 #3
|
||||
{ADAR_4_CS_3V3_GPIO_Port, ADAR_4_CS_3V3_Pin} // ADAR1000 #4
|
||||
};
|
||||
|
||||
// ADAR1000 Vector Modulator lookup tables (128-state phase grid, 2.8125 deg step).
|
||||
//
|
||||
@@ -255,15 +255,15 @@ bool ADAR1000Manager::setBeamAngle(float angle_degrees, BeamDirection direction)
|
||||
|
||||
for (uint8_t dev = 0; dev < devices_.size(); ++dev) {
|
||||
for (uint8_t ch = 0; ch < 4; ++ch) {
|
||||
if (direction == BeamDirection::TX) {
|
||||
adarSetTxPhase(dev, ch + 1, phase_settings[ch], BROADCAST_OFF);
|
||||
adarSetTxVgaGain(dev, ch + 1, kDefaultTxVgaGain, BROADCAST_OFF);
|
||||
} else {
|
||||
adarSetRxPhase(dev, ch + 1, phase_settings[ch], BROADCAST_OFF);
|
||||
adarSetRxVgaGain(dev, ch + 1, kDefaultRxVgaGain, BROADCAST_OFF);
|
||||
}
|
||||
}
|
||||
}
|
||||
if (direction == BeamDirection::TX) {
|
||||
adarSetTxPhase(dev, ch + 1, phase_settings[ch], BROADCAST_OFF);
|
||||
adarSetTxVgaGain(dev, ch + 1, kDefaultTxVgaGain, BROADCAST_OFF);
|
||||
} else {
|
||||
adarSetRxPhase(dev, ch + 1, phase_settings[ch], BROADCAST_OFF);
|
||||
adarSetRxVgaGain(dev, ch + 1, kDefaultRxVgaGain, BROADCAST_OFF);
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -522,7 +522,7 @@ bool ADAR1000Manager::initializeADTR1107Sequence() {
|
||||
HAL_UART_Transmit(&huart3, success, sizeof(success) - 1, 1000);
|
||||
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
bool ADAR1000Manager::setAllDevicesTXMode() {
|
||||
DIAG("BF", "setAllDevicesTXMode(): ADTR1107 -> TX, then configure ADAR1000s");
|
||||
@@ -648,7 +648,7 @@ void ADAR1000Manager::setADTR1107Mode(BeamDirection direction) {
|
||||
}
|
||||
DIAG("BF", " ADTR1107 RX mode complete");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void ADAR1000Manager::setADTR1107Control(bool tx_mode) {
|
||||
DIAG("BF", "setADTR1107Control(%s): setting TR switch on all %u devices, settling %lu us",
|
||||
@@ -727,7 +727,7 @@ void ADAR1000Manager::setLNABias(bool enable) {
|
||||
for (uint8_t dev = 0; dev < devices_.size(); ++dev) {
|
||||
adarWrite(dev, REG_LNA_BIAS_ON, lna_bias, BROADCAST_OFF);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void ADAR1000Manager::delayUs(uint32_t microseconds) {
|
||||
// Simple implementation - for F7 @ 216MHz, each loop ~7 cycles ≈ 0.032us
|
||||
@@ -903,11 +903,11 @@ void ADAR1000Manager::adarSetTxVgaGain(uint8_t deviceIndex, uint8_t channel, uin
|
||||
adarWrite(deviceIndex, REG_LOAD_WORKING, LD_WRK_REGS_LDTX_OVERRIDE, broadcast);
|
||||
}
|
||||
|
||||
void ADAR1000Manager::adarSetTxBias(uint8_t deviceIndex, uint8_t broadcast) {
|
||||
adarWrite(deviceIndex, REG_BIAS_CURRENT_TX, kTxBiasCurrent, broadcast);
|
||||
adarWrite(deviceIndex, REG_BIAS_CURRENT_TX_DRV, kTxDriverBiasCurrent, broadcast);
|
||||
adarWrite(deviceIndex, REG_LOAD_WORKING, 0x2, broadcast);
|
||||
}
|
||||
void ADAR1000Manager::adarSetTxBias(uint8_t deviceIndex, uint8_t broadcast) {
|
||||
adarWrite(deviceIndex, REG_BIAS_CURRENT_TX, kTxBiasCurrent, broadcast);
|
||||
adarWrite(deviceIndex, REG_BIAS_CURRENT_TX_DRV, kTxDriverBiasCurrent, broadcast);
|
||||
adarWrite(deviceIndex, REG_LOAD_WORKING, 0x2, broadcast);
|
||||
}
|
||||
|
||||
uint8_t ADAR1000Manager::adarAdcRead(uint8_t deviceIndex, uint8_t broadcast) {
|
||||
adarWrite(deviceIndex, REG_ADC_CONTROL, ADAR1000_ADC_ST_CONV, broadcast);
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user