Example: Runway Crosswind & Wind-Shift Operability
====================================================

An airport decides whether a runway is usable over the coming day. Two
things close it: the **crosswind** component of the wind (its speed times
the sine of the angle between the wind and the runway heading) exceeding an
aircraft limit, and abrupt **wind shifts** between hours. Both are driven
by a wind that evolves over the day — and wind direction is a *circular*
quantity that wanders across North and often switches between two
prevailing regimes. The full runnable source is
``examples/runway_crosswind/model.py``.

This example is the process side of prophys end-to-end: a coupled
stochastic wind, a circular direction with a bimodal step, and operability
read off the marginal of the day's worst crosswind — its expectation, its
tail, and the exceedance probability that decides downtime.

.. image:: /_static/generated/runway_crosswind_operability.png
   :width: 95%
   :align: center

Step 1 — The wind as a coupled process
------------------------------------------

Wind speed is a mean-reverting Ornstein–Uhlenbeck process. Wind direction
is a :func:`~prophys.circular_autoregressive` process reverting to the
prevailing heading, with a **two-component Gaussian mixture** innovation
for the backing/veering regimes — a step distribution that is bimodal, not
Gaussian. Sampling such a process is exactly what a mixture's derived
quantile makes possible (a bare mixture had no ``quantile`` before):

.. code-block:: python

   import prophys as prp

   speed = prp.ornstein_uhlenbeck(init=8.0, mean=8.0, theta=0.5, sigma=2.5, n_steps=24)

   direction = prp.circular_autoregressive(
       init=250.0, mean=250.0, phi=0.85,
       innovation=prp.Mixture([prp.Gaussian(-6.0, 4.0), prp.Gaussian(6.0, 4.0)],
                             logits=[0.0, 0.0]),
       n_steps=24, period=360.0,
   )

The circular geometry matters: reverting toward the prevailing 250° from a
direction near North follows the short arc across the 0/360 seam, and every
sampled hour stays a valid heading in ``[0, 360)`` — visible as the smooth
wander in the left panel above.

The two are coupled into one :class:`~prophys.JointProcess` so gusty hours
and shifting hours co-occur, then exposed as a single joint random variable
whose components share the underlying draw:

.. code-block:: python

   wind_process = prp.JointProcess([speed, direction], corr_raw=[0.35])
   wind = prp.JointRandomVariable("wind", wind_process, components=["speed", "direction"])

   speed_path = prp.clip(wind["speed"], 0.0, None)   # (..., 24) m/s
   direction_path = wind["direction"]                 # (..., 24) deg

Step 2 — Operability read off the wind
-----------------------------------------

The instantaneous crosswind on a runway of a given heading is
:math:`V\,\lvert\sin(\theta - \text{heading})\rvert`, and the operability
question is about the *worst* crosswind of the day — a maximum along the
time axis, observed through small anemometer noise so it is a proper
:class:`~prophys.domain.UncertainAttribute`:

.. code-block:: python

   def crosswind_component(heading_deg):
       return prp.abs_(speed_path * prp.sin(prp.deg2rad(direction_path - heading_deg)))

   max_crosswind = prp.UncertainAttribute(
       "max_crosswind_07_25",
       prp.Gaussian(mean=time_max(crosswind_component(250.0)), sigma=anemometer_sigma),
       unit="m/s",
   )

The largest hour-to-hour wind shift is a **circular** first difference, so
a swing across North is measured as the short arc rather than a spurious
near-360° jump:

.. code-block:: python

   from prophys.processes import Circular
   hourly_shift = prp.time_diff(direction_path, geometry=Circular(360.0))

Step 3 — Expectation, tail, and downtime probability
-------------------------------------------------------

Compiling marginalizes the wind process by Monte Carlo. For each runway the
model reports the expected worst crosswind, its **Conditional
Value-at-Risk** (the mean of the worst 5 % of days — the number that
actually governs operations), and the downtime probability that the day's
worst crosswind exceeds the limit. The exceedance is computed two ways that
agree: directly by Monte Carlo, and from a Gaussian fitted to the sampled
maxima via its first-class ``cdf`` (:math:`P(X > \ell) = 1 - F(\ell)`):

.. code-block:: python

   compiled = model.compile(mode="eval", n_samples=4000)

   e     = compiled.expectation("max_crosswind_07_25")
   cvar  = compiled.cvar("max_crosswind_07_25", alpha=0.95)
   p_mc  = 1.0 - compiled.prob("max_crosswind_07_25", CROSSWIND_LIMIT)

   samples = compiled.sample("max_crosswind_07_25", key, (4000,))
   fitted  = prp.Gaussian(float(samples.mean()), float(samples.std()))
   p_cdf   = 1.0 - float(fitted.cdf(CROSSWIND_LIMIT).evaluate({}))

For the modelled airport the runway aligned with the prevailing wind sees a
low expected crosswind and essentially no downtime, while the
near-perpendicular runway straddles the limit with ~50 % downtime (the two
distributions in the right panel above) — so the model picks the aligned
runway.

Primitives exercised
-----------------------

.. list-table::
   :header-rows: 1

   * - Category
     - Used for
   * - :func:`~prophys.ornstein_uhlenbeck`, :func:`~prophys.circular_autoregressive`
     - Wind speed (mean-reverting) and direction (circular, seam-correct)
   * - :class:`~prophys.distributions.Mixture`
     - Bimodal backing/veering hourly step — usable as an innovation via its derived quantile
   * - :class:`~prophys.JointProcess`, :class:`~prophys.domain.RandomVariable`
     - Coupling gusty and shifting hours; exposing a shared joint draw
   * - :func:`~prophys.time_diff` with :class:`~prophys.Circular`
     - Seam-safe hour-to-hour wind-shift magnitude
   * - ``expectation`` / ``cvar`` / ``prob`` and :class:`~prophys.distributions.Gaussian` ``cdf``/``quantile``
     - Expected & tail crosswind; downtime probability by Monte Carlo and by CDF

Running it
-----------

.. code-block:: bash

   python examples/runway_crosswind/model.py